Heterogeneous integrated circuit thermal mechanical stress collaborative analysis method and system

Abstract

The invention provides a heterogeneous integrated circuit thermomechanical stress collaborative analysis method and system, and relates to the technical field of integrated circuit thermomechanical analysis, and the method comprises the steps: building a finite element analysis grid through obtaining geometric, material and temperature stress data, building a causal reasoning network guided by physical knowledge to recognize a thermal stress evolution path, and obtaining a thermal stress evolution path; and forming a thermal stress evolution analysis map, and constructing a dynamic network to carry out system characteristic analysis so as to realize thermal stress instability early warning.

Description

Technical Field

[0001] The present invention relates to the technical field of integrated circuit thermomechanical analysis, and in particular to a method and system for collaborative analysis of thermomechanical stress of heterogeneous integrated circuits. Background Art

[0002] As electronic products develop toward high integration and multifunctionality, heterogeneous integrated circuit technology has become an important approach to addressing the performance bottlenecks of traditional integrated circuits. Heterogeneous integrated circuits achieve functional complementarity and performance improvements by integrating devices made of different materials and using different processes. However, due to differences in thermal expansion coefficients, heterogeneous material systems are prone to thermal stress during operation, leading to problems such as chip deformation, interface delamination, and reduced reliability. As the operating power density of integrated circuits continues to increase, thermomechanical stress issues are becoming increasingly serious, becoming a key factor affecting the reliability and lifespan of heterogeneous integrated circuits. Traditional thermomechanical stress analysis methods mainly rely on finite element analysis, which performs numerical simulations by constructing geometric models, setting boundary conditions, and setting material parameters. However, these methods have significant limitations under complex dynamic conditions. Existing thermomechanical stress analysis technologies fail to fully consider the bidirectional coupling effects of thermal and force fields, making it difficult to accurately simulate the thermal stress transfer mechanism at heterogeneous material interfaces. This leads to insufficient stress prediction accuracy under complex working conditions, a lack of in-depth research into the dynamic evolution characteristics of thermal stress systems, an inability to effectively identify key causal relationships and evolution paths in the system, and difficulty in accurately modeling the nonlinear behavior of thermal stress systems. Furthermore, there is a lack of an instability warning mechanism based on the system's dynamic characteristics, making it impossible to predict potential thermal stress instability risks in advance. Therefore, a solution is urgently needed to solve the problems existing in the prior art. Summary of the Invention

[0003] The embodiments of the present invention provide a method and system for collaborative analysis of thermal-mechanical stress in heterogeneous integrated circuits, which can at least solve some of the problems existing in the prior art.

[0004] A first aspect of an embodiment of the present invention provides a method for collaboratively analyzing thermal-mechanical stress in a heterogeneous integrated circuit, comprising: Obtaining the three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, constructing a finite element analysis mesh based on the three-dimensional geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; Constructing a causal reasoning network guided by physical knowledge, embedding the thermal-mechanical coupling coefficient matrix into the causal reasoning network, and identifying key causal transmission paths of thermal stress evolution through a structured causal discovery method, establishing a spatiotemporal causal mapping relationship of the thermal field, constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map; Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method and a thermal stress instability warning criterion is established. The key state variables are determined by embedding dimensionality analysis, and a dynamic prediction equation of the thermal stress system is established. The target working state is input into the dynamic prediction equation, and the system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified in combination with the thermal stress evolution analysis map, and a dynamic warning report is generated according to the thermal stress instability warning criterion.

[0005] In an optional embodiment, Constructing a finite element analysis mesh based on 3D geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on temperature and stress data include: Acquiring three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, wherein the three-dimensional geometric parameters include the dimensions of the interlayer interconnect structure, the distribution position of through-silicon vias, the chip layout, and the package boundary dimensions of the heterogeneous integrated circuit; the material parameters include the elastic modulus, Poisson's ratio, thermal expansion coefficient, and thermal conductivity of each layer of the medium; the temperature data includes the temperature distribution field and the temperature gradient field; and the stress data includes the thermal stress field distribution and the stress gradient field distribution; Dividing the heterogeneous integrated circuit into regular hexahedral grid units according to the three-dimensional geometric parameters and material parameters, setting grid nodes at interlayer interconnect structures and through-silicon via locations, and constructing a spatial mapping relationship between the nodes using an isoparametric unit interpolation function to obtain a finite element analysis grid; The temperature data and stress data are mapped to the grid nodes of the finite element analysis grid, the temperature-stress response relationship of the grid nodes is calculated, and a thermomechanical coupling coefficient matrix representing the degree of coupling between the temperature field and the stress field is generated.

[0006] In an optional embodiment, A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, and the key causal transmission path of thermal stress evolution is identified through a structured causal discovery method. The spatiotemporal causal mapping relationship of the thermal field is established, including: A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, an adjacency matrix is ​​constructed by setting a coupling strength threshold and performing a temporal causal test, and a directed graph containing thermal stress system nodes is constructed based on the adjacency matrix, wherein the edge set of the directed graph represents the thermal coupling relationship between the nodes; Establishing constraint conditions including the coupling relationship between the temperature field and the stress field based on the directed graph, constructing an optimization objective function including a temperature-stress data fitting error term and a physical constraint violation degree term, and obtaining the thermal stress evolution relationship between the nodes according to the optimization objective function; Based on the thermal stress evolution relationship between the nodes, the conditional probability of the node states at adjacent moments is calculated to obtain a transfer entropy that characterizes the degree of information transmission between the nodes. The product of the transfer entropy and the thermal coupling strength is used as the weight value of the corresponding edge in the directed graph. The maximum flow-minimum cut algorithm is combined with the weight value to identify the key causal transmission path of thermal stress evolution based on the weight value; The key causal transmission path is projected into the space-time domain to establish a space-time causal mapping function. The space-time causal mapping function maps the space-time coordinates of the source point to the space-time coordinates of the target point. The space-time causal mapping function is estimated by a nonlinear regression method to obtain a space-time causal mapping relationship that characterizes the space-time propagation characteristics of thermal stress.

[0007] In an optional embodiment, Constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map includes: Constructing an orthogonal variational encoder network structure including a mean encoder and a variance encoder, wherein the mean encoder encodes the thermal stress field state to obtain a mean parameter of the latent space variable, and the variance encoder encodes the thermal stress field state to obtain a variance parameter of the latent space variable, and constructing a probability distribution of the latent space variable based on the mean parameter and the variance parameter; A loss function is constructed based on the orthogonal variational encoder, wherein the loss function includes a reconstruction loss term for reconstructing the thermal stress field state using the latent space variables, a KL divergence loss term between the probability distribution of the latent space variables and the prior distribution, and an orthogonality constraint loss term between latent features, and the orthogonal decoupling feature of the thermal stress field state is obtained by minimizing the loss function; Constructing a thermal field embedding module, wherein the thermal field embedding module performs a first nonlinear transformation on the thermal field tensor using a preset first weight matrix and a first bias vector, performs a second nonlinear transformation on the output of the first nonlinear transformation using a second weight matrix and a second bias vector to obtain a thermal field embedding representation, combines the orthogonal decoupling feature with the thermal field embedding representation, and constructs a stress response function based on the spatiotemporal causal mapping relationship; The stress response function is used to calculate the stress response of each node in the finite element analysis grid, generate multiple sets of thermal stress evolution data including node position, time and stress value, and construct a thermal stress evolution analysis map, wherein the nodes in the thermal stress evolution analysis map represent stress states, the edges represent state transition relationships, and the weights represent transition probabilities.

[0008] In an optional embodiment, Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained, and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method, including: The thermal stress field sampling points in the thermal stress evolution data are constructed as a node set, the dynamic associations between the node sets are constructed as an edge set, and the weight matrix between the nodes is calculated by the Gaussian kernel function to obtain the dynamic network of the thermal stress system; Calculating the root mean square deviation of the node phase angle and the average phase angle in the dynamic network according to the inter-node weight matrix to obtain an order parameter, and analyzing the time evolution curve of the order parameter to obtain the evolution characteristics of the thermal stress system; Constructing a phase space reconstruction module for the thermal stress system according to the evolution characteristics, calculating an optimal time delay parameter of the time series by a mutual information method, determining a minimum embedding dimension of a phase space reconstruction vector by combining a phase space false neighbor method, and generating the phase space reconstruction vector based on the optimal time delay parameter and the minimum embedding dimension; The phase space reconstruction vector is input into a recursive quantitative analysis method to calculate the exponential growth rate of the disturbance in the orthogonal direction in the phase space to obtain the local Lyapunov exponent spectrum. The bifurcation point of the thermal stress system is identified by analyzing the sign change of the local Lyapunov exponent spectrum, and the attractor of the thermal stress system is calculated by constructing a correlation integral based on the Heisenberg step function.

[0009] In an optional embodiment, Establish thermal stress instability warning criteria, determine key state variables through embedded dimensional analysis, and establish the dynamic prediction equation of the thermal stress system including: Based on pre-acquired bifurcation points and attractors, pre-collected time-varying data of a thermal stress field are acquired, a simplicial complex sequence is constructed based on a distance metric between nodes in the thermal stress field, the rank of the dimensional homology group of the simplicial complex sequence is calculated, a Betti number sequence corresponding to each moment is obtained, and feature matching is performed with the Betti number sequence at the bifurcation point, a state space reconstruction dimension is determined through embedded dimensionality analysis, and the time-varying data is reconstructed to obtain a reconstructed state space, and a weighted summation of the Betti number sequences of different dimensions is performed to obtain key state variables; Recording the generation time and extinction time of the topological feature in the Betti number sequence, calculating the stability index of the persistent homology feature based on the time difference between the generation time and the extinction time in combination with the feature decay function, and constructing a persistent homology graph based on the bifurcation point and the attractor; Extracting a duration distribution of topological features based on the persistence homology graph, and fusing the duration distribution with the key state variables to obtain enhanced state features; Linearly combining the enhanced state feature, the first-order time derivative of the stability index, and the second-order time derivative of the stability index to generate an early warning index and construct a thermal stress instability early warning criterion; The state vector and its derivative in the reconstructed state space are obtained, and a dynamic prediction equation containing a basis function is constructed by combining the bifurcation point and the attractor with the thermal stress instability warning criterion as a constraint condition.

[0010] In an optional embodiment, The target working state is input into the dynamic prediction equation. The system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified by combining the thermal stress evolution analysis map. A dynamic warning report is generated according to the thermal stress instability warning criteria, including: Input the target working state into the dynamic prediction equation to calculate the evolution trajectory of the thermal stress system; Based on matching the topological features of the evolution trajectory with a thermal stress evolution analysis map, the thermal stress evolution analysis map includes the topological feature distribution of historical instability data, and identifying the potential instability region by calculating the similarity between the Betti number sequence of the evolution trajectory and the instability features in the thermal stress evolution analysis map; According to the distribution characteristics of the potential instability area and the thermal stress instability warning criteria, a dynamic warning report including the instability risk level and warning time is generated in combination with the spatiotemporal distribution characteristics of the potential instability area.

[0011] A second aspect of an embodiment of the present invention provides a system for collaboratively analyzing thermal mechanical stress of heterogeneous integrated circuits, comprising: The first unit is used to obtain the three-dimensional geometric parameters, material parameters, temperature data and stress data of the heterogeneous integrated circuit, construct a finite element analysis grid based on the three-dimensional geometric parameters and material parameters, and calculate the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; The second unit is used to construct a causal reasoning network guided by physical knowledge, embed the thermal-mechanical coupling coefficient matrix into the causal reasoning network, identify the key causal transmission paths of thermal stress evolution through a structured causal discovery method, establish a spatiotemporal causal mapping relationship of the thermal field, construct a stress response function based on the spatiotemporal causal mapping relationship, combine the stress response function with the finite element analysis grid to generate multiple sets of thermal stress evolution data, and form a thermal stress evolution analysis map; The third unit is used to construct the thermal stress system into a dynamic network based on the thermal stress evolution data, obtain the evolution characteristics of the thermal stress system and construct a phase space reconstruction module, identify the bifurcation points and attractors of the thermal stress system through a recursive quantitative analysis method and establish a thermal stress instability warning criterion, determine the key state variables through embedded dimensional analysis, and establish a dynamic prediction equation for the thermal stress system; The fourth unit is used to input the target working state into the dynamic prediction equation, determine the system state prediction result by calculating the evolution trajectory of the thermal stress system, identify potential unstable areas in combination with the thermal stress evolution analysis map, and generate a dynamic warning report based on the thermal stress instability warning criteria.

[0012] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0013] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0014] In the present invention, by constructing a causal reasoning network guided by physical knowledge, the key causal paths of thermal stress evolution can be accurately identified, the spatiotemporal causal mapping relationship of the thermal field can be established, and the complex thermo-mechanical coupling mechanism in heterogeneous integrated circuits can be effectively captured. Based on the phase space reconstruction and recursive quantitative analysis of the thermal stress system, the bifurcation points and attractors of the system can be scientifically identified, and a thermal stress instability warning criterion can be established. The dynamic prediction of the thermal stress state of the heterogeneous integrated circuit is realized, the accuracy and reliability of the thermal stress analysis are improved, and strong support is provided for the reliability design of the integrated circuit. By combining finite element analysis with dynamic system theory, a complete thermal stress evolution analysis map and dynamic warning mechanism are constructed, which can not only monitor the thermal stress state of the integrated circuit in real time, but also predict potential instability areas in advance, thereby effectively preventing circuit failure caused by thermal stress and greatly improving the reliability and service life of the heterogeneous integrated circuit. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1Schematic diagram of the process of the collaborative analysis method of thermal mechanical stress of heterogeneous integrated circuits according to an embodiment of the present invention; Figure 2 This is a simulation diagram of the thermal stress field distribution in the through-silicon via region of the heterogeneous integrated circuit thermal-mechanical stress collaborative analysis method according to an embodiment of the present invention; Figure 3 This is a comparison chart of the Betti number sequence feature extraction and bifurcation identification performance of the collaborative analysis method for thermal mechanical stress of heterogeneous integrated circuits according to an embodiment of the present invention. DETAILED DESCRIPTION

[0016] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0017] The following specific embodiments are used to describe the technical solution of the present invention in detail. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in detail in some embodiments.

[0018] Figure 1 FIG. 1 is a flow chart of a method for collaboratively analyzing thermal mechanical stress of heterogeneous integrated circuits according to an embodiment of the present invention. Figure 1 As shown, the method includes: Obtaining the three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, constructing a finite element analysis mesh based on the three-dimensional geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; Constructing a causal reasoning network guided by physical knowledge, embedding the thermal-mechanical coupling coefficient matrix into the causal reasoning network, and identifying key causal transmission paths of thermal stress evolution through a structured causal discovery method, establishing a spatiotemporal causal mapping relationship of the thermal field, constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map; Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method and a thermal stress instability warning criterion is established. The key state variables are determined by embedding dimensionality analysis, and a dynamic prediction equation of the thermal stress system is established. The target working state is input into the dynamic prediction equation, and the system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified in combination with the thermal stress evolution analysis map, and a dynamic warning report is generated according to the thermal stress instability warning criterion.

[0019] In an optional embodiment, Constructing a finite element analysis mesh based on 3D geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on temperature and stress data include: Acquiring three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, wherein the three-dimensional geometric parameters include the dimensions of the interlayer interconnect structure, the distribution position of through-silicon vias, the chip layout, and the package boundary dimensions of the heterogeneous integrated circuit; the material parameters include the elastic modulus, Poisson's ratio, thermal expansion coefficient, and thermal conductivity of each layer of the medium; the temperature data includes the temperature distribution field and the temperature gradient field; and the stress data includes the thermal stress field distribution and the stress gradient field distribution; Dividing the heterogeneous integrated circuit into regular hexahedral grid units according to the three-dimensional geometric parameters and material parameters, setting grid nodes at interlayer interconnect structures and through-silicon via locations, and constructing a spatial mapping relationship between the nodes using an isoparametric unit interpolation function to obtain a finite element analysis grid; The temperature data and stress data are mapped to the grid nodes of the finite element analysis grid, the temperature-stress response relationship of the grid nodes is calculated, and a thermomechanical coupling coefficient matrix representing the degree of coupling between the temperature field and the stress field is generated.

[0020] Obtain the three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit. The three-dimensional geometric parameters include the dimensions of the interlayer interconnect structure, the distribution of through-silicon vias (TSVs), the chip layout, and the package boundary dimensions. For example, a typical three-dimensional heterogeneous integrated circuit may include multiple silicon wafers, each with a thickness of 50 microns, an interlayer interconnect structure with a thickness of 10 microns, and TSVs with a diameter of 5 microns and a height of 60 microns, evenly distributed on the chip plane at a spacing of 100 microns by 100 microns. The package boundary dimensions are 10 mm by 10 mm by 2 mm. Material parameters include the elastic modulus, Poisson's ratio, thermal expansion coefficient, and thermal conductivity of each dielectric layer. For example, the elastic modulus of silicon material is 130 GPa, Poisson's ratio is 0.28, thermal expansion coefficient is 2.6 × 10^-6 / K, and thermal conductivity is 130 W / (m·K); the elastic modulus of the silicon dioxide dielectric layer is 70 GPa, Poisson's ratio is 0.17, thermal expansion coefficient is 0.5 × 10^-6 / K, and thermal conductivity is 1.4 W / (m·K); the elastic modulus of the copper interconnect is 110 GPa, Poisson's ratio is 0.34, thermal expansion coefficient is 17 × 10^-6 / K, and thermal conductivity is 400 W / (m·K). Temperature data includes temperature distribution field and temperature gradient field, which can be obtained through thermal conduction analysis or experimental measurement. Stress data includes thermal stress field distribution and stress gradient field distribution, which can be obtained through thermo-mechanical coupling analysis or experimental measurement.

[0021] Based on the acquired three-dimensional geometric and material parameters, the heterogeneous integrated circuit is divided into regular hexahedral mesh cells. During the meshing process, mesh nodes are set according to the locations of interlayer interconnects and through-silicon vias (TSVs) to ensure accurate representation of key structures. For example, for a heterogeneous integrated circuit consisting of three silicon wafers, a 100×100×50 mesh cell grid can be created in the X, Y, and Z directions, with at least four mesh nodes set at each TSV location to accurately represent the geometry. For the interlayer interconnects, at least three mesh nodes are set in the thickness direction to capture stress gradient variations. After meshing is completed, isoparametric element interpolation functions are used to establish the spatial mapping relationship between nodes. For hexahedral elements, either linear or quadratic interpolation functions can be used. Linear interpolation uses the eight vertex nodes of the element, while quadratic interpolation uses 20 nodes, including intermediate nodes. Using interpolation functions, physical quantities (such as temperature and displacement) at any point can be represented as linear combinations of the node quantities, thereby establishing a complete finite element analysis mesh.

[0022] The temperature data and stress data are mapped to the grid nodes of the finite element analysis grid. For data obtained from other sources, if the data points do not completely correspond to the grid nodes, interpolation processing is performed. Methods such as nearest neighbor interpolation, linear interpolation, or spline interpolation can be used. For example, suppose that data with a temperature of 85°C is obtained at a certain measurement point, and this point is not on the grid node, the temperature value is distributed to the surrounding grid nodes by distance weighted averaging. For each grid node, its temperature value T and stress tensor σ (including six components σxx, σyy, σzz, τxy, τyz, τxz) are recorded. In actual applications, the temperature may vary from room temperature 25°C to operating temperature 85°C. The corresponding thermal stress can reach 200MPa near the silicon via, while it is only 50MPa in the area far away from the silicon via.

[0023] The temperature-stress response relationship of the grid nodes is calculated to generate a thermomechanical coupling coefficient matrix that characterizes the degree of coupling between the temperature and stress fields. Specifically, for each grid node, the relationship between temperature change and stress change is analyzed. By analyzing multiple sets of temperature-stress data, a linear or nonlinear relationship can be determined. For example, consider multiple intermediate points (such as 35°C, 45°C, 55°C, 65°C, and 75°C) during a temperature change from 25°C to 85°C. The stress values ​​corresponding to each temperature point can then be recorded and fitted to these data points to form a temperature-stress response curve. For a linear relationship, the thermomechanical coupling coefficient matrix C can be expressed as the ratio of the stress change to the temperature change. For example, if the temperature increases from 25°C to 85°C (a 60°C increase), and the σxx of a node increases from 10 MPa to 130 MPa (a 120 MPa increase), then the thermomechanical coupling coefficient of that node in the xx direction is 2 MPa / °C. Similar calculations are performed for all nodes in the integrated circuit to obtain a complete thermomechanical coupling coefficient matrix. The thermomechanical coupling coefficient matrix reflects the degree to which temperature changes at different locations affect the stress field and can be used to guide the thermal management and reliability design of integrated circuits.

[0024] In this embodiment, by obtaining the complete parameters of the heterogeneous integrated circuit, a comprehensive description of the complex heterogeneous integrated structure is achieved. The regular hexahedral grid division method is adopted to avoid the singular solutions that may be caused by irregular grids. By calculating the temperature-stress response relationship of the grid nodes, the degree of coupling between the temperature field and the stress field is quantitatively characterized, providing an accurate data basis for subsequent instability warning.

[0025] In an optional embodiment, A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, and the key causal transmission path of thermal stress evolution is identified through a structured causal discovery method. The spatiotemporal causal mapping relationship of the thermal field is established, including: A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, an adjacency matrix is ​​constructed by setting a coupling strength threshold and performing a temporal causal test, and a directed graph containing thermal stress system nodes is constructed based on the adjacency matrix, wherein the edge set of the directed graph represents the thermal coupling relationship between the nodes; Establishing constraint conditions including the coupling relationship between the temperature field and the stress field based on the directed graph, constructing an optimization objective function including a temperature-stress data fitting error term and a physical constraint violation degree term, and obtaining the thermal stress evolution relationship between the nodes according to the optimization objective function; Based on the thermal stress evolution relationship between the nodes, the conditional probability of the node states at adjacent moments is calculated to obtain a transfer entropy that characterizes the degree of information transmission between the nodes. The product of the transfer entropy and the thermal coupling strength is used as the weight value of the corresponding edge in the directed graph. The maximum flow-minimum cut algorithm is combined with the weight value to identify the key causal transmission path of thermal stress evolution based on the weight value; The key causal transmission path is projected into the space-time domain to establish a space-time causal mapping function. The space-time causal mapping function maps the space-time coordinates of the source point to the space-time coordinates of the target point. The space-time causal mapping function is estimated by a nonlinear regression method to obtain a space-time causal mapping relationship that characterizes the space-time propagation characteristics of thermal stress.

[0026] A 30×30 thermal coupling coefficient matrix was obtained. For a high-temperature equipment monitoring system, data was collected from 20 temperature sensors and 10 stress sensors. This coupling coefficient matrix was embedded in a causal inference network. A coupling strength threshold of 0.35 was set, and weak coupling relationships below this threshold were eliminated. A time-series Granger causality test was performed on node pairs that passed the threshold, with a significance level of 0.05. A binary adjacency matrix was constructed. In this binary adjacency matrix, a value of 1 indicates a causal relationship, and a value of 0 indicates no causal relationship. Based on this adjacency matrix, a directed graph containing 30 nodes was constructed, where nodes represent temperature or stress measurement points and edges represent the thermal coupling relationships between nodes.

[0027] Based on the above directed graph, constraints including the coupling relationship between the temperature field and the stress field are established. Under the guidance of physical knowledge, the constraints established include: heat conduction law constraints, thermal expansion law constraints, and material constitutive relationship constraints. For example, for two adjacent temperature nodes, the heat flow transfer follows Fourier's heat conduction law; for the influence of the temperature node to the stress node, it follows the relationship between thermal expansion and stress generation. An optimization objective function is constructed, which includes a data fitting error term and a physical constraint violation degree term. The data fitting error term calculates the mean square error between the predicted value and the measured value, and the error is 0.082 on the test data set; the physical constraint violation degree term is obtained by calculating the sum of the squares of the residuals of the physical equation, and the weight coefficient is set to 0.5. The alternating direction multiplier method is used to solve the optimization problem, and it converges after 200 iterations, and the thermal stress evolution relationship between the nodes is obtained, which is expressed in the form of a state transfer matrix.

[0028] Based on the obtained inter-node thermal stress evolution relationships, the conditional probabilities of node states at adjacent moments are calculated. For each pair of causally related nodes, the state of the source node at the previous moment is probabilistically analyzed with the state of the target node at the next moment, and the transfer entropy is calculated. The transfer entropy quantifies the amount of information flow from the source node to the target node, with a value ranging from 0 to 1.08. The product of the transfer entropy and the thermal-mechanical coupling coefficient is used as the weight of the corresponding edge in the directed graph to construct a weighted directed graph. The maximum flow-minimum cut algorithm is applied to identify the key causal transfer paths of thermal stress evolution. The network source is set as the temperature node in the initial heated area, and the sink is the stress evolution node of interest. The critical path of thermal stress evolution is determined by finding the minimum cut set. Experimental results show that in the test case, the critical path consists of five temperature nodes and three stress nodes, forming a main thermal stress propagation chain.

[0029] The identified key causal transfer paths were projected into the spatiotemporal domain to establish a spatiotemporal causal mapping function. This mapping function maps the spatiotemporal coordinates (xs, ys, zs, ts) of the source point to the spatiotemporal coordinates (xt, yt, zt, tt) of the target point. Taking into account the nonlinear characteristics of thermal stress transfer, a radial basis function network was used for nonlinear regression to estimate the spatiotemporal causal mapping function. The training data consisted of 200 source-target coordinate pairs obtained from simulations. After 10-fold cross-validation, the average prediction error of the mapping function was 1.2 mm and 0.3 s.

[0030] In this embodiment, the spatiotemporal causal mapping relationship can accurately predict the propagation path and time delay of thermal disturbances in the system. By constructing a causal reasoning network guided by physical knowledge, the accurate identification of the causal structure in the thermal stress system is achieved. The key causal transmission path of thermal stress evolution can be effectively identified, and an accurate spatiotemporal causal mapping relationship can be established, providing strong support for the design optimization, fault diagnosis and predictive maintenance of the thermal system.

[0031] In an optional embodiment, Constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map includes: Constructing an orthogonal variational encoder network structure including a mean encoder and a variance encoder, wherein the mean encoder encodes the thermal stress field state to obtain a mean parameter of the latent space variable, and the variance encoder encodes the thermal stress field state to obtain a variance parameter of the latent space variable, and constructing a probability distribution of the latent space variable based on the mean parameter and the variance parameter; A loss function is constructed based on the orthogonal variational encoder, wherein the loss function includes a reconstruction loss term for reconstructing the thermal stress field state using the latent space variables, a KL divergence loss term between the probability distribution of the latent space variables and the prior distribution, and an orthogonality constraint loss term between latent features, and the orthogonal decoupling feature of the thermal stress field state is obtained by minimizing the loss function; Constructing a thermal field embedding module, wherein the thermal field embedding module performs a first nonlinear transformation on the thermal field tensor using a preset first weight matrix and a first bias vector, performs a second nonlinear transformation on the output of the first nonlinear transformation using a second weight matrix and a second bias vector to obtain a thermal field embedding representation, combines the orthogonal decoupling feature with the thermal field embedding representation, and constructs a stress response function based on the spatiotemporal causal mapping relationship; The stress response function is used to calculate the stress response of each node in the finite element analysis grid, generate multiple sets of thermal stress evolution data including node position, time and stress value, and construct a thermal stress evolution analysis map, wherein the nodes in the thermal stress evolution analysis map represent stress states, the edges represent state transition relationships, and the weights represent transition probabilities.

[0032] The orthogonal variational encoder network preprocesses the input thermal stress field state data, mapping the data to a uniform numerical range through normalization. The mean encoder uses convolutional layers to extract local features of the thermal stress field. This is then input into a fully connected layer after dimensionality reduction through a pooling layer. This layer-by-layer feature compression produces the mean parameter of the latent space. The variance encoder uses the same network structure but with independent parameters to extract features and compress the dimensionality of the input data, outputting the variance parameter of the latent space. The mean and variance parameters are then fed into a reparameterization module, which generates latent variable samples that follow a specific distribution through random sampling.

[0033] During training, the loss function is constructed in multiple steps. The first involves calculating the reconstruction loss, which uses the decoder to reconstruct the original thermal stress field state from the latent variables and calculates the error between the reconstruction and the input data. The second step is to calculate the divergence loss, comparing the distribution of the latent space variables with the standard normal distribution to quantify the difference between the two distributions. Finally, the orthogonality constraint loss is calculated to extract the correlation matrix between the latent features and achieve feature decoupling by minimizing the off-diagonal elements. In each iteration, the gradient of the loss function with respect to the network parameters is calculated using the backpropagation algorithm, and the network parameters are updated using the optimizer to gradually reduce the total loss value.

[0034] The processing flow of the thermal field embedding module consists of two stages: feature extraction and nonlinear transformation. The thermal field tensor is input into the attention mechanism module, which calculates the correlation between thermal features at different locations and generates an attention weight map. Based on the attention weight map, the thermal features are weighted and aggregated to highlight information in important areas. The first nonlinear transformation uses a batch normalization layer and an activation function to perform a nonlinear mapping on the weighted features. The output of the nonlinear mapping is fused with the original features through a residual connection module. The second nonlinear transformation also includes batch normalization and activation operations, but uses a different activation function to enhance the expressiveness of the features.

[0035] When fusing orthogonal decoupled features with the embedded representation of the thermal field, the two sets of features are dimensionally aligned and mapped to the same feature space using a projection layer. A gating mechanism dynamically adjusts the importance weights of the two features to achieve adaptive feature fusion. Based on the fused features and the spatiotemporal causal mapping relationship, a stress response function of a multilayer perceptron structure is constructed, which can predict stress distribution under given conditions.

[0036] When generating thermal stress evolution data, the spatial position coordinates and time information of each node in the finite element analysis grid are collected, and the stress response function is input to obtain the node stress value. By changing the input conditions, the stress response of the node under different working conditions is obtained to form a data set containing position, time and stress values. When constructing the thermal stress evolution analysis map, each stress state is represented as a node in the graph, and the node attributes contain the main components of the stress tensor. The transition relationship between states is determined based on time series correlation analysis, and the states at adjacent moments are connected with directed edges. The weight of the edge is determined by statistically analyzing the frequency and jump probability of state transitions, and the weight value reflects the possibility of state transition. Finally, the map is sparsely processed, and the edges with smaller weights are removed, retaining the main state transition paths.

[0037] For example, thermal stress field data is collected in the through-silicon via (TSV) region of a heterogeneous integrated circuit as input. A mean encoder encodes this input data into a mean vector in a ten-dimensional latent space using a three-layer neural network, while a variance encoder outputs the corresponding variance vector. Based on these parameters, a probability distribution model is constructed using a multidimensional normal distribution as the latent variable.

[0038] When constructing the loss function, the reconstruction loss is expressed as mean squared error, the divergence loss is calculated based on the relative entropy between distributions, and the orthogonality constraint is enforced by the square sum of the off-diagonal elements of the feature correlation matrix. After optimization training, orthogonal decoupled features are obtained that can characterize stress amplitude, gradient, and direction.

[0039] The thermal field embedding module adopts a two-layer neural network structure. The first layer uses a weight matrix and bias vector containing 128 neurons, and the second layer uses a weight matrix and bias vector containing 64 neurons. The final result is a low-dimensional embedding representation of the thermal field. Combined with the orthogonal decoupling feature, a response function is constructed to predict the stress distribution around silicon vias.

[0040] Using a stress response function, the stress response of mesh nodes surrounding TSVs under varying temperature loads is calculated, generating thousands of sets of thermal stress evolution data. In this constructed thermal stress evolution analysis graph, different stress states are represented by nodes, and state transitions between adjacent moments are connected by weighted edges, with the weights determined by the frequency of state transitions.

[0041] In this embodiment, an orthogonal variational encoder network is constructed to achieve efficient feature extraction and dimensionality reduction of the thermal stress field state. Deep feature extraction of thermal field information is achieved through two nonlinear transformations. The spatiotemporal causal mapping relationship is considered when constructing the stress response function, allowing the prediction results to reflect the physical evolution laws of the thermal stress field. A systematic description of the thermal stress evolution process is achieved by constructing a thermal stress evolution analysis map. In existing technologies, the thermal stress field is modeled as a whole, making it difficult to effectively extract key features from the thermal stress field. This results in insufficient prediction accuracy and a lack of consideration of the physical mechanisms of the thermal stress field. This makes it impossible to ensure that the prediction results conform to physical laws, affecting the reliability of the prediction results. Furthermore, the computational efficiency is low when processing high-dimensional thermal stress data, making it difficult to meet the needs of real-time monitoring. This embodiment improves the accuracy and efficiency of thermal stress field feature extraction, reduces the consumption of computing resources, enhances the model's perception of local thermal stress changes, improves prediction accuracy, ensures that the prediction results conform to physical laws, enhances the reliability of the model, realizes the systematic characterization of the thermal stress evolution process, and provides effective support for instability warning.

[0042] Figure 2This is a simulation diagram of the thermal stress field distribution in the through-silicon via region of the heterogeneous integrated circuit thermal mechanical stress collaborative analysis method according to an embodiment of the present invention, comparing the predicted results of this technical solution with the actual measured data. In the figure, concentric ellipses represent the contour distribution of the thermal stress field, and the central gray circular area represents the TSV core position. The blue dotted line series represents the measured stress distribution under 25°C temperature conditions, which are 43.2MPa, 27.5MPa, 12.8MPa and 5.3MPa from the inside to the outside. As the distance from the TSV core increases, the stress value shows a significant attenuation trend, which is consistent with the theoretical expectations of the stress field. The red dotted line represents the measured data under high temperature conditions of 125°C, showing a clear elliptical distribution, indicating that the thermal stress field is directional and is mainly affected by material anisotropy and temperature gradient (8.3°C / mm). The stress values ​​under high temperature conditions are 127.6MPa, 96.3MPa, 62.1MPa and 31.7MPa from the inside to the outside, which is about 3 times the stress value at room temperature, showing a significant effect of temperature on the stress field.

[0043] The green solid line shows the prediction results of this technical solution based on the orthogonal variational encoder network. The four contour lines predicted by this technical solution are highly consistent with the red high-temperature measured data. After calculation, the stress values ​​from the inner layer to the outer layer predicted by this technical solution are 125.3MPa, 94.8MPa, 61.5MPa and 30.9MPa, respectively. The relative errors with the measured values ​​(127.6MPa, 96.3MPa, 62.1MPa and 31.7MPa) are all controlled within 2%, among which the prediction accuracy of the middle and outer layers is as high as 99.0%. This high-precision prediction capability stems from the thermal field embedding module and stress response function constructed in this technical solution, which can effectively capture the nonlinear distribution characteristics of the thermal stress field.

[0044] The purple dots in the figure mark four hotspots located at different locations on the TSV-silicon substrate interface: top, right, bottom, and left. The stress values ​​in these hotspots are significantly higher than those in the surrounding areas, with the top hotspot reaching a high of 156.2 MPa, the bottom hotspot 135.4 MPa, and the right and left hotspots 118.7 MPa and 122.3 MPa, respectively. These hotspots are often the most susceptible to reliability issues during thermal cycling in TSV structures and require special attention.

[0045] The yellow dots represent the node locations in the finite element analysis grid, totaling 25 nodes, and the stress response data of these nodes under different temperature conditions are collected.

[0046] Compared to existing technologies such as the traditional finite element method (average error of 8.7%) and physical information neural networks (average error of 5.3%), this solution achieves an average prediction error of only 1.8%. The improvement in prediction accuracy is particularly pronounced in areas with large stress gradient variations. Combined with the stress response function constructed by embedding the thermal field, this approach provides a powerful tool for thermal stress analysis and life assessment of key integrated circuit structures.

[0047] In an optional embodiment, Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained, and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method, including: The thermal stress field sampling points in the thermal stress evolution data are constructed as a node set, the dynamic associations between the node sets are constructed as an edge set, and the weight matrix between the nodes is calculated by the Gaussian kernel function to obtain the dynamic network of the thermal stress system; Calculating the root mean square deviation of the node phase angle and the average phase angle in the dynamic network according to the inter-node weight matrix to obtain an order parameter, and analyzing the time evolution curve of the order parameter to obtain the evolution characteristics of the thermal stress system; Constructing a phase space reconstruction module for the thermal stress system according to the evolution characteristics, calculating an optimal time delay parameter of the time series by a mutual information method, determining a minimum embedding dimension of a phase space reconstruction vector by combining a phase space false neighbor method, and generating the phase space reconstruction vector based on the optimal time delay parameter and the minimum embedding dimension; The phase space reconstruction vector is input into a recursive quantitative analysis method to calculate the exponential growth rate of the disturbance in the orthogonal direction in the phase space to obtain the local Lyapunov exponent spectrum. The bifurcation point of the thermal stress system is identified by analyzing the sign change of the local Lyapunov exponent spectrum, and the attractor of the thermal stress system is calculated by constructing a correlation integral based on the Heisenberg step function.

[0048] Thermal stress field sampling points are extracted from thermal stress evolution data as a node set, each containing spatial location and stress state information. The dynamic associations between any two nodes are analyzed, and an edge set is established based on the similarity of node states. A Gaussian kernel function is used to nonlinearly map the state differences between nodes, generating a weight matrix that reflects the strength of the associations between nodes. A dynamic network of the thermal stress system is constructed based on the node set, edge set, and weight matrix. This network characterizes the dynamic coupling relationships between sampling points in the thermal stress field.

[0049] Based on the constructed dynamic network, the order parameter is calculated, and the phase angle information of each node is obtained. The phase angles of all nodes are compared with the average phase angle of the network, and the root mean square deviation is calculated to obtain the order parameter. Time series analysis of the order parameter is performed to observe its evolutionary trends and fluctuation characteristics. By analyzing the mutation points, stable ranges, and fluctuation patterns of the order parameter, the evolutionary characteristics of the thermal stress system are extracted, including information such as the degree of synchronization, stability, and critical behavior.

[0050] A phase space reconstruction module is constructed based on the extracted evolutionary features to preprocess the time series data and eliminate noise and trend terms. The autocorrelation characteristics of the time series are analyzed using the mutual information method. By traversing different delay values, the delay corresponding to the first local minimum of the mutual information function is selected as the optimal delay parameter. The phase space false neighbor method is used to determine the dimension of the reconstructed vector. By gradually increasing the embedding dimension, the proportion of false neighbor points under each dimension is calculated. When the false neighbor ratio falls below a threshold, the corresponding dimension is the minimum embedding dimension. Based on the determined optimal delay parameter and minimum embedding dimension, the original time series is reconstructed in phase space to generate a reconstructed vector.

[0051] The phase space reconstruction vector is input into a recursive quantitative analysis method. A reference orbit is selected in phase space. For each point on the reference orbit, a search is performed for adjacent orbits within its neighborhood. The divergence rates of the adjacent orbits in different orthogonal directions are calculated, and the local Lyapunov exponent spectrum is obtained through exponential fitting. The sign change characteristics of the exponential spectrum are analyzed. When the maximum Lyapunov exponent changes from negative to positive, the corresponding state point is the bifurcation point. After identifying the bifurcation point, the convergence behavior of the orbits in phase space is analyzed by constructing a correlation integral based on the Heigenstep function, and the attractor characteristics of the thermal stress system are determined.

[0052] For example, during the thermal stress evolution of a heterogeneous integrated circuit, temperature and stress monitoring points in the area surrounding through-silicon vias (TSVs) were selected as network nodes. A similarity analysis was performed on the stress states between each pair of monitoring points, and a weight matrix was calculated using a Gaussian kernel function. The constructed dynamic network showed strong coupling between nodes in the region surrounding the TSVs.

[0053] By calculating the phase angles of each node in the network, analysis revealed that the order parameter exhibits significant fluctuations when the temperature changes dramatically. The order parameter evolution curve indicates that at certain critical temperature points, the system's coherence undergoes significant changes, potentially signaling early signs of thermal stress instability.

[0054] During the phase space reconstruction process, the optimal time delay parameter was determined using the mutual information method, and the minimum embedding dimension was determined using the phase space pseudo-nearest neighbor method. The reconstructed phase space orbits clearly demonstrate the dynamic behavior of the thermal stress system, including both periodic and sudden changes.

[0055] Recursive quantitative analysis of the local Lyapunov exponent spectrum revealed that within a specific temperature range, the maximum Lyapunov exponent switches between positive and negative, indicating a potential system bifurcation within that range. Combined with correlation integral analysis to determine the system's attractor characteristics, it was found that near the bifurcation point, the system's orbit jumps from one attractor to another. This jump often indicates the onset of thermal stress instability.

[0056] In this embodiment, by constructing a dynamic network based on nodes and edges, an accurate description of the spatial distribution characteristics of the thermal stress field is achieved. The dynamic correlation information extracted from the thermal stress evolution data can reflect the interaction intensity between the various parts of the system, providing a reliable data basis for subsequent analysis. Through the calculation and analysis of order parameters, a quantitative characterization of the overall evolution characteristics of the thermal stress system is achieved. The phase space reconstruction module ensures that the reconstructed phase space can accurately reflect the intrinsic dynamic characteristics of the system through the optimal selection of the optimal time delay parameter and the minimum embedding dimension. The complete technical chain from data acquisition to feature extraction, from local analysis to overall characterization, through multi-level data analysis and feature extraction, realizes a systematic characterization of the dynamic behavior of the thermal stress system, providing a reliable theoretical basis and technical support for thermal stress instability warning.

[0057] In an optional embodiment, Establish thermal stress instability warning criteria, determine key state variables through embedded dimensional analysis, and establish the dynamic prediction equation of the thermal stress system including: Based on pre-acquired bifurcation points and attractors, pre-collected time-varying data of a thermal stress field are acquired, a simplicial complex sequence is constructed based on a distance metric between nodes in the thermal stress field, the rank of the dimensional homology group of the simplicial complex sequence is calculated, a Betti number sequence corresponding to each moment is obtained, and feature matching is performed with the Betti number sequence at the bifurcation point, a state space reconstruction dimension is determined through embedded dimensionality analysis, and the time-varying data is reconstructed to obtain a reconstructed state space, and a weighted summation of the Betti number sequences of different dimensions is performed to obtain key state variables; Recording the generation time and extinction time of the topological feature in the Betti number sequence, calculating the stability index of the persistent homology feature based on the time difference between the generation time and the extinction time in combination with the feature decay function, and constructing a persistent homology graph based on the bifurcation point and the attractor; Extracting a duration distribution of topological features based on the persistence homology graph, and fusing the duration distribution with the key state variables to obtain enhanced state features; Linearly combining the enhanced state feature, the first-order time derivative of the stability index, and the second-order time derivative of the stability index to generate an early warning index and construct a thermal stress instability early warning criterion; The state vector and its derivative in the reconstructed state space are obtained, and a dynamic prediction equation containing a basis function is constructed by combining the bifurcation point and the attractor with the thermal stress instability warning criterion as a constraint condition.

[0058] Preprocess the previously collected time-varying thermal stress field data. Segment the data using a sliding window, with the data within each window treated as a processing unit. Each processing unit is normalized to eliminate dimensionality effects. A distance matrix is ​​calculated between nodes in the thermal stress field, and an adaptive threshold method is used to determine the distance threshold. Based on the distance threshold, a simplicial complex is constructed layer by layer, starting with a zero-dimensional simplex and gradually forming a simplicial complex sequence by adding higher-dimensional surfaces.

[0059] A topological analysis is performed on the constructed simplicial complex sequence. A boundary operator matrix is ​​established, and the rank of the homology group in each dimension, i.e., the Betti number, is calculated by calculating the dimensions of the kernel space and image space. This calculation process is repeated for each time window, resulting in a sequence of Betti number evolution over the time series. The Betti number sequence is matched with a sequence of characteristic features at pre-marked bifurcation points, and the dynamic time warping algorithm is used to calculate sequence similarity and extract the evolutionary characteristics of the topological structure.

[0060] The state-space reconstruction dimensionality analysis is performed, and the time delay parameters are calculated using the mutual information method. The minimum embedding dimension is determined using the false nearest neighbor algorithm. Based on the obtained reconstruction parameters, the original time-varying data is reconstructed using the time-delay coordinate method to generate orbits in the high-dimensional phase space. In the reconstructed space, the reconstruction quality is evaluated using principal component analysis, and the reconstruction parameters are adjusted as necessary to meet the requirements.

[0061] The importance of Betti number sequences in different dimensions is quantitatively analyzed. The information gain criterion is used to calculate the contribution of each dimension's Betti number to the prediction of system evolution. The contributions are normalized and used as weight coefficients. A weighted summation is performed to obtain key state variables, which comprehensively reflect the multidimensional topological characteristics of the system.

[0062] We perform a lifecycle analysis of topological features in the Betti number sequence. We use a persistent homology algorithm to track the evolution of each topological feature, recording the exact moments of its generation and extinction. We construct a feature decay function that accounts for both the duration and intensity of the feature and calculates its stability. We combine this tracked feature evolution with information about bifurcation points and attractors to construct a persistent homology graph, where nodes represent topological features and edges represent the evolutionary relationships between features.

[0063] Feature statistics are extracted from the persistence homology graph. A probability density estimate is performed on the duration of topological features to obtain the duration distribution characteristics. Kernel density estimation is used to avoid the limitations of parameter distribution assumptions. The duration distribution characteristics are nonlinearly fused with key state variables, and a deep autoencoder is used to achieve dimensionality reduction and reconstruction of the feature space to obtain enhanced state characteristics.

[0064] Numerically differentiate the stability indicators. A central difference scheme is used to calculate the first-order derivative, reflecting the rate of change of the system state. A second-order difference scheme is used to calculate the second-order derivative, characterizing the acceleration characteristics of the system state. Polynomial fitting is used to smooth the derivatives to reduce the noise introduced by numerical differentiation.

[0065] Design a combination strategy for early warning indicators. Using a sliding time window approach, calculate the correlation between enhanced state features and derivative features within each window. Determine the weight coefficients of the linear combination based on the correlation analysis results to construct early warning indicators. Analyze historical instability cases to determine thresholds for early warning indicators and establish a multi-level early warning criterion.

[0066] Dynamic information is extracted from the reconstructed state space. The evolution direction and velocity of the state vector are calculated, and a local linear prediction model is constructed. Radial basis functions are selected as basis functions, and their centers are determined using a clustering algorithm. The warning criteria are converted into constraints, which are then incorporated into the optimization objective of the prediction equation using regularization methods. A gradient descent algorithm is used to optimize the parameters of the prediction equation until the prediction error converges below a set threshold.

[0067] The performance of the prediction equation is evaluated through cross-validation. Prediction results at multiple time scales are compared with the actual evolution trajectory to verify the accuracy and stability of the prediction. Based on the validation results, the parameters of the prediction equation are fine-tuned to improve the prediction accuracy.

[0068] For example, in thermal stress monitoring of heterogeneous integrated circuits, stress monitoring data from the area surrounding through-silicon vias (TSVs) were analyzed. By calculating the Euclidean distance between monitoring points, a simplicial complex sequence was constructed, and the Betti numbers for each dimension were calculated. It was observed that when the temperature fluctuated drastically, the zero-dimensional Betti number increased significantly, indicating a change in the system's connectivity.

[0069] The monitoring data were reconstructed into a state space, and analysis confirmed that the reconstruction dimension was four. In the reconstructed space, the weights of the Betti numbers of different dimensions were determined based on their contributions to bifurcation prediction, resulting in key variables reflecting the system state.

[0070] Analysis of the temporal characteristics of the Betti number sequence reveals that topological features persist for a long time near bifurcation points. Stability indices are calculated using characteristic decay functions, and persistent homology diagrams are constructed to clearly demonstrate the evolution of these topological features.

[0071] The persistence homology graph shows that the duration distribution of the topological feature changes significantly before instability. By integrating this feature with key state variables, the enhanced state feature obtained is more sensitive to the signs of system instability.

[0072] By calculating the derivatives of the stability index, it was found that the first-order and second-order derivatives exhibit characteristic changes before instability. A warning indicator was constructed based on the enhanced state characteristics, and a warning threshold was set.

[0073] By leveraging the dynamic information in the reconstructed state space and selecting appropriate basis functions, a prediction equation is constructed. This equation can accurately predict the stress evolution trend around the TSV, providing a basis for instability warning.

[0074] In this embodiment, by constructing a simplicial complex sequence and calculating the Betti number, the topological structure of the thermal stress field is accurately characterized. The distance measurement method with an adaptive threshold is used to improve the accuracy of topological feature extraction. Through state space reconstruction and dimensional homology analysis, the efficient extraction of system dynamic characteristics is achieved. By introducing basis functions and constraints, the generalization ability of the prediction model is improved. Through the deep fusion of topological features and dynamic features, the early warning system's ability to identify precursors of instability is enhanced.

[0075] Figure 3This figure compares the Betti number sequence feature extraction and bifurcation identification performance of the collaborative analysis method for thermomechanical stress in heterogeneous integrated circuits according to an embodiment of the present invention. It demonstrates the evolution of Betti number sequences of different dimensions during thermal stress instability and compares them with traditional bifurcation identification methods. In the stable region (t = 0-120s), the 0-dimensional Betti number β0 remains within the range of 37±3, the 1-dimensional Betti number β1 is within the range of 15±2, and the 2-dimensional Betti number β2 is within the range of 8±1. As the system approaches the first bifurcation point (t = 152s), β0 suddenly increases to 53 (an increase of 43.2%), β1 increases to 27 (an increase of 80%), and β2 increases to 16 (an increase of 100%), indicating a significant change in the topological structure. At the second bifurcation point (t = 205s), the topological characteristics change even more dramatically, with β0 reaching 72 (an increase of 94.6%), β1 reaching 42 (an increase of 180%), and β2 reaching 23 (an increase of 187.5%). In contrast, the traditional Lyapunov exponent method and entropy analysis methods showed little response at the first bifurcation point, with changes of 12.3% (from 26 to 31) and 18.7% (from 33 to 38), respectively. Significant changes only occurred at the second bifurcation point, at 53.6% (from 31 to 48) and 61.2% (from 38 to 52), respectively. This technical solution, through weighted fusion of the Betti number sequence (weights of 0.45, 0.35, and 0.20, obtained a significant change in the key state variables starting 20 seconds before the bifurcation point, with a sensitivity 2.1 times higher than that of traditional methods.

[0076] In 10 experiments under different operating conditions, the proposed technology achieved an average lead time of 18.7 seconds for predicting bifurcation points, with a prediction error of ±2.3 seconds. The traditional Lyapunov exponent method and entropy analysis method achieved lead times of 7.2 seconds and 9.5 seconds, with prediction errors of ±5.8 seconds and ±4.7 seconds, respectively. The Betti number sequence reflects changes in the system's topological structure, exhibiting a more pronounced trend in the early stages of thermal stress instability. Traditional methods, primarily based on the system's dynamic or statistical characteristics, are unable to effectively capture early topological changes. This technology, combining multi-dimensional Betti number features, can comprehensively characterize the system's topological evolution in different dimensions, providing more dimensional feature information for thermal stress instability warning. Experimental results show that under varying operating conditions with a temperature gradient change rate ranging from 0.5°C / s to 3.5°C / s, the proposed technology consistently outperforms traditional methods in bifurcation prediction, demonstrating the advantages of computational topology-based analysis methods in the evolutionary analysis of complex nonlinear systems.

[0077] In an optional embodiment, The target working state is input into the dynamic prediction equation. The system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified by combining the thermal stress evolution analysis map. A dynamic warning report is generated according to the thermal stress instability warning criteria, including: Input the target working state into the dynamic prediction equation to calculate the evolution trajectory of the thermal stress system; Based on matching the topological features of the evolution trajectory with a thermal stress evolution analysis map, the thermal stress evolution analysis map includes the topological feature distribution of historical instability data, and identifying the potential instability region by calculating the similarity between the Betti number sequence of the evolution trajectory and the instability features in the thermal stress evolution analysis map; According to the distribution characteristics of the potential instability area and the thermal stress instability warning criteria, a dynamic warning report including the instability risk level and warning time is generated in combination with the spatiotemporal distribution characteristics of the potential instability area.

[0078] The target operating state data is preprocessed and standardized as input to the dynamic prediction equation. The prediction equation is iteratively calculated, with progressively increasing time steps, to obtain a series of predicted state points. These state points are then connected in chronological order to form the system's evolution trajectory. The trajectory is then interpolated and smoothed to ensure continuity and smoothness.

[0079] Topological features are extracted from the acquired evolution trajectory. Neighborhood relationships are constructed for the trajectory points, and a simplicial complex is formed based on a distance threshold. Betti numbers are calculated for each dimension to obtain a sequence of Betti numbers that characterize the trajectory topology. The topological distribution of historical instability cases is extracted from the thermal stress evolution analysis map. Cluster analysis is performed on the instability data in the map to identify typical instability patterns.

[0080] Perform feature matching analysis. A dynamic programming algorithm is used to calculate the similarity between the predicted trajectory's Betti number sequence and the instability features in the atlas. A similarity threshold is set to filter out trajectory segments that are highly similar to the historical instability features. These trajectory segments are spatially located and time-stamped to determine the extent of the potential instability region.

[0081] Analyze the distribution characteristics of potential unstable areas. Calculate the spatial density of unstable areas to identify hotspots of concentrated instability risk. Conduct time series analysis of unstable areas to determine the temporal evolution of instability risk. Match the distribution characteristics with early warning criteria and conduct a comprehensive assessment based on multiple early warning indicators.

[0082] Construct a risk assessment model across time and space. Combine spatial distribution characteristics with temporal evolution characteristics to establish a quantitative standard for risk levels. Categorize instability risks into multiple levels based on the severity of the risk. Set a corresponding warning time window for each risk level to ensure timely warnings.

[0083] Generate dynamic early warning reports, visualize risk assessment results, and create risk distribution heat maps. Integrate risk level information with warning time information to form hierarchical early warning information. Design different warning prompts based on the urgency of the risk. Organize warning information in chronological order to form a complete early warning report.

[0084] For example, an instability early warning analysis was performed on the TSV region of a heterogeneous integrated circuit. The current temperature distribution and stress state were input into a prediction equation, and the evolution trajectory over the next 24 hours was calculated. The prediction results showed that stress concentration was likely to occur at the top of the TSV.

[0085] A topological analysis of the predicted trajectory was performed, and a series of Betti numbers from zero to two dimensions was calculated. This was then matched with historical instability cases from the thermal stress evolution analysis map. The matching results showed that the predicted trajectory exhibited topological features highly similar to those of historical instability cases in certain time periods.

[0086] Similarity analysis identified three potential instability regions: the top, sidewall, and bottom contact areas of the TSV. The Betti number sequences in these regions were similar to historical instability cases exceeding a preset threshold.

[0087] Analysis revealed that the top area presents the highest risk of instability, with instability predicted to occur within the next eight hours. The sidewall area exhibits a moderate risk, with a warning window of 12 to 16 hours. While the bottom contact area also exhibits signs of instability, the risk level is relatively low.

[0088] Based on the risk assessment results, a dynamic early warning report is generated, including risk levels and warning times. The report uses heat maps to visually display the risk distribution and provides graded early warning information. For high-risk areas, preventive measures are recommended within four hours. The report also provides a timeline of risk evolution to help determine the optimal time for intervention.

[0089] In this embodiment, the system evolution trajectory is calculated using the dynamic prediction equation to accurately predict the future state of the thermal stress system. The analysis method based on topological feature matching is used to achieve efficient comparison with historical instability cases. By calculating the similarity of the Betti number sequence, areas with instability risks are accurately identified, thereby improving the accuracy of the early warning.

[0090] A second aspect of an embodiment of the present invention provides a system for collaboratively analyzing thermal mechanical stress of heterogeneous integrated circuits, comprising: The first unit is used to obtain the three-dimensional geometric parameters, material parameters, temperature data and stress data of the heterogeneous integrated circuit, construct a finite element analysis grid based on the three-dimensional geometric parameters and material parameters, and calculate the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; The second unit is used to construct a causal reasoning network guided by physical knowledge, embed the thermal-mechanical coupling coefficient matrix into the causal reasoning network, identify the key causal transmission paths of thermal stress evolution through a structured causal discovery method, establish a spatiotemporal causal mapping relationship of the thermal field, construct a stress response function based on the spatiotemporal causal mapping relationship, combine the stress response function with the finite element analysis grid to generate multiple sets of thermal stress evolution data, and form a thermal stress evolution analysis map; The third unit is used to construct the thermal stress system into a dynamic network based on the thermal stress evolution data, obtain the evolution characteristics of the thermal stress system and construct a phase space reconstruction module, identify the bifurcation points and attractors of the thermal stress system through a recursive quantitative analysis method and establish a thermal stress instability warning criterion, determine the key state variables through embedded dimensional analysis, and establish a dynamic prediction equation for the thermal stress system; The fourth unit is used to input the target working state into the dynamic prediction equation, determine the system state prediction result by calculating the evolution trajectory of the thermal stress system, identify potential unstable areas in combination with the thermal stress evolution analysis map, and generate a dynamic warning report based on the thermal stress instability warning criteria.

[0091] According to a third aspect of an embodiment of the present invention, an electronic device is provided, including: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to call the instructions stored in the memory to execute the aforementioned method.

[0092] According to a fourth aspect of an embodiment of the present invention, a computer-readable storage medium is provided, on which computer program instructions are stored. When the computer program instructions are executed by a processor, the method described above is implemented.

[0093] The present invention may be a method, an apparatus, a system and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for executing various aspects of the present invention.

[0094] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements 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. A collaborative analysis method for thermal-mechanical stress in heterogeneous integrated circuits, characterized in that: include: Obtaining the three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, constructing a finite element analysis mesh based on the three-dimensional geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; Constructing a causal reasoning network guided by physical knowledge, embedding the thermal-mechanical coupling coefficient matrix into the causal reasoning network, and identifying key causal transmission paths of thermal stress evolution through a structured causal discovery method, establishing a spatiotemporal causal mapping relationship of the thermal field, constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map; Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method and a thermal stress instability warning criterion is established. The key state variables are determined by embedding dimensionality analysis, and a dynamic prediction equation of the thermal stress system is established. The target working state is input into the dynamic prediction equation, and the system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified in combination with the thermal stress evolution analysis map, and a dynamic warning report is generated according to the thermal stress instability warning criterion.

2. The method according to claim 1, characterized in that Constructing a finite element analysis mesh based on 3D geometric parameters and material parameters, and calculating the thermal-mechanical coupling coefficient matrix based on temperature and stress data include: Acquiring three-dimensional geometric parameters, material parameters, temperature data, and stress data of the heterogeneous integrated circuit, wherein the three-dimensional geometric parameters include the dimensions of the interlayer interconnect structure, the distribution position of through-silicon vias, the chip layout, and the package boundary dimensions of the heterogeneous integrated circuit; the material parameters include the elastic modulus, Poisson's ratio, thermal expansion coefficient, and thermal conductivity of each layer of the medium; the temperature data includes the temperature distribution field and the temperature gradient field; and the stress data includes the thermal stress field distribution and the stress gradient field distribution; Dividing the heterogeneous integrated circuit into regular hexahedral grid units according to the three-dimensional geometric parameters and material parameters, setting grid nodes at interlayer interconnect structures and through-silicon via locations, and constructing a spatial mapping relationship between the nodes using an isoparametric unit interpolation function to obtain a finite element analysis grid; The temperature data and stress data are mapped to the grid nodes of the finite element analysis grid, the temperature-stress response relationship of the grid nodes is calculated, and a thermomechanical coupling coefficient matrix representing the degree of coupling between the temperature field and the stress field is generated.

3. The method according to claim 1, characterized in that A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, and the key causal transmission path of thermal stress evolution is identified through a structured causal discovery method. The spatiotemporal causal mapping relationship of the thermal field is established, including: A causal reasoning network guided by physical knowledge is constructed, the thermal coupling coefficient matrix is ​​embedded in the causal reasoning network, an adjacency matrix is ​​constructed by setting a coupling strength threshold and performing a temporal causal test, and a directed graph containing thermal stress system nodes is constructed based on the adjacency matrix, wherein the edge set of the directed graph represents the thermal coupling relationship between the nodes; Establishing constraint conditions including the coupling relationship between the temperature field and the stress field based on the directed graph, constructing an optimization objective function including a temperature-stress data fitting error term and a physical constraint violation degree term, and obtaining the thermal stress evolution relationship between the nodes according to the optimization objective function; Based on the thermal stress evolution relationship between the nodes, the conditional probability of the node states at adjacent moments is calculated to obtain a transfer entropy that characterizes the degree of information transmission between the nodes. The product of the transfer entropy and the thermal coupling strength is used as the weight value of the corresponding edge in the directed graph. The maximum flow-minimum cut algorithm is combined with the weight value to identify the key causal transmission path of thermal stress evolution based on the weight value; The key causal transmission path is projected into the space-time domain to establish a space-time causal mapping function. The space-time causal mapping function maps the space-time coordinates of the source point to the space-time coordinates of the target point. The space-time causal mapping function is estimated by a nonlinear regression method to obtain a space-time causal mapping relationship that characterizes the space-time propagation characteristics of thermal stress.

4. The method according to claim 1, wherein Constructing a stress response function based on the spatiotemporal causal mapping relationship, combining the stress response function with a finite element analysis grid to generate multiple sets of thermal stress evolution data, and forming a thermal stress evolution analysis map includes: Constructing an orthogonal variational encoder network structure including a mean encoder and a variance encoder, wherein the mean encoder encodes the thermal stress field state to obtain a mean parameter of the latent space variable, and the variance encoder encodes the thermal stress field state to obtain a variance parameter of the latent space variable, and constructing a probability distribution of the latent space variable based on the mean parameter and the variance parameter; A loss function is constructed based on the orthogonal variational encoder, wherein the loss function includes a reconstruction loss term for reconstructing the thermal stress field state using the latent space variables, a KL divergence loss term between the probability distribution of the latent space variables and the prior distribution, and an orthogonality constraint loss term between latent features, and the orthogonal decoupling feature of the thermal stress field state is obtained by minimizing the loss function; Constructing a thermal field embedding module, wherein the thermal field embedding module performs a first nonlinear transformation on the thermal field tensor using a preset first weight matrix and a first bias vector, performs a second nonlinear transformation on the output of the first nonlinear transformation using a second weight matrix and a second bias vector to obtain a thermal field embedding representation, combines the orthogonal decoupling feature with the thermal field embedding representation, and constructs a stress response function based on the spatiotemporal causal mapping relationship; The stress response function is used to calculate the stress response of each node in the finite element analysis grid, generate multiple sets of thermal stress evolution data including node position, time and stress value, and construct a thermal stress evolution analysis map, wherein the nodes in the thermal stress evolution analysis map represent stress states, the edges represent state transition relationships, and the weights represent transition probabilities.

5. The method according to claim 1, wherein Based on the thermal stress evolution data, the thermal stress system is constructed as a dynamic network, the evolution characteristics of the thermal stress system are obtained, and a phase space reconstruction module is constructed. The bifurcation points and attractors of the thermal stress system are identified by a recursive quantitative analysis method, including: The thermal stress field sampling points in the thermal stress evolution data are constructed as a node set, the dynamic associations between the node sets are constructed as an edge set, and the weight matrix between the nodes is calculated by the Gaussian kernel function to obtain the dynamic network of the thermal stress system; Calculating the root mean square deviation of the node phase angle and the average phase angle in the dynamic network according to the inter-node weight matrix to obtain an order parameter, and analyzing the time evolution curve of the order parameter to obtain the evolution characteristics of the thermal stress system; Constructing a phase space reconstruction module for the thermal stress system according to the evolution characteristics, calculating an optimal time delay parameter of the time series by a mutual information method, determining a minimum embedding dimension of a phase space reconstruction vector by combining a phase space false neighbor method, and generating the phase space reconstruction vector based on the optimal time delay parameter and the minimum embedding dimension; The phase space reconstruction vector is input into a recursive quantitative analysis method to calculate the exponential growth rate of the disturbance in the orthogonal direction in the phase space to obtain the local Lyapunov exponent spectrum. The bifurcation point of the thermal stress system is identified by analyzing the sign change of the local Lyapunov exponent spectrum, and the attractor of the thermal stress system is calculated by constructing a correlation integral based on the Heisenberg step function.

6. The method according to claim 1, wherein Establish thermal stress instability warning criteria, determine key state variables through embedded dimensional analysis, and establish the dynamic prediction equation of the thermal stress system including: Based on pre-acquired bifurcation points and attractors, pre-collected time-varying data of a thermal stress field are acquired, a simplicial complex sequence is constructed based on a distance metric between nodes in the thermal stress field, the rank of the dimensional homology group of the simplicial complex sequence is calculated, a Betti number sequence corresponding to each moment is obtained, and feature matching is performed with the Betti number sequence at the bifurcation point, a state space reconstruction dimension is determined through embedded dimensionality analysis, and the time-varying data is reconstructed to obtain a reconstructed state space, and a weighted summation of the Betti number sequences of different dimensions is performed to obtain key state variables; Recording the generation time and extinction time of the topological feature in the Betti number sequence, calculating the stability index of the persistent homology feature based on the time difference between the generation time and the extinction time in combination with the feature decay function, and constructing a persistent homology graph based on the bifurcation point and the attractor; Extracting a duration distribution of topological features based on the persistence homology graph, and fusing the duration distribution with the key state variables to obtain enhanced state features; Linearly combining the enhanced state feature, the first-order time derivative of the stability index, and the second-order time derivative of the stability index to generate an early warning index and construct a thermal stress instability early warning criterion; The state vector and its derivative in the reconstructed state space are obtained, and a dynamic prediction equation containing a basis function is constructed by combining the bifurcation point and the attractor with the thermal stress instability warning criterion as a constraint condition.

7. The method according to claim 1, characterized in that The target working state is input into the dynamic prediction equation. The system state prediction result is determined by calculating the evolution trajectory of the thermal stress system. The potential instability area is identified by combining the thermal stress evolution analysis map. A dynamic warning report is generated according to the thermal stress instability warning criteria, including: Input the target working state into the dynamic prediction equation to calculate the evolution trajectory of the thermal stress system; Based on matching the topological features of the evolution trajectory with a thermal stress evolution analysis map, the thermal stress evolution analysis map includes the topological feature distribution of historical instability data, and identifying the potential instability region by calculating the similarity between the Betti number sequence of the evolution trajectory and the instability features in the thermal stress evolution analysis map; According to the distribution characteristics of the potential instability area and the thermal stress instability warning criteria, a dynamic warning report including the instability risk level and warning time is generated in combination with the spatiotemporal distribution characteristics of the potential instability area.

8. A system for collaborative analysis of thermal mechanical stress of heterogeneous integrated circuits, for implementing the method according to any one of claims 1 to 7, characterized in that: include: The first unit is used to obtain the three-dimensional geometric parameters, material parameters, temperature data and stress data of the heterogeneous integrated circuit, construct a finite element analysis grid based on the three-dimensional geometric parameters and material parameters, and calculate the thermal-mechanical coupling coefficient matrix based on the temperature data and stress data; The second unit is used to construct a causal reasoning network guided by physical knowledge, embed the thermal-mechanical coupling coefficient matrix into the causal reasoning network, identify the key causal transmission paths of thermal stress evolution through a structured causal discovery method, establish a spatiotemporal causal mapping relationship of the thermal field, construct a stress response function based on the spatiotemporal causal mapping relationship, combine the stress response function with the finite element analysis grid to generate multiple sets of thermal stress evolution data, and form a thermal stress evolution analysis map; The third unit is used to construct the thermal stress system into a dynamic network based on the thermal stress evolution data, obtain the evolution characteristics of the thermal stress system and construct a phase space reconstruction module, identify the bifurcation points and attractors of the thermal stress system through a recursive quantitative analysis method and establish a thermal stress instability warning criterion, determine the key state variables through embedded dimensional analysis, and establish a dynamic prediction equation for the thermal stress system; The fourth unit is used to input the target working state into the dynamic prediction equation, determine the system state prediction result by calculating the evolution trajectory of the thermal stress system, identify potential unstable areas in combination with the thermal stress evolution analysis map, and generate a dynamic warning report based on the thermal stress instability warning criteria.

9. An electronic device, characterized in that: include: processor; a memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that: When the computer program instructions are executed by a processor, the method according to any one of claims 1 to 7 is implemented.

Images