Lightweight multi-heat-source thermal model construction method and system for integrated microsystem

Abstract

The invention discloses a lightweight multi-heat-source thermal model construction method and system for an integrated microsystem, and the method comprises the following steps: S1, constructing a nonlinear material attribute incidence matrix Mmat through self-adaptive analysis of a grid file, extraction of a topological relation between tetrahedral elements, node coordinates and surface elements; s2, setting a transient heat source and a heat exchange boundary according to a tetrahedron unit space relationship, assembling a unit matrix based on an Mmat matrix, and integrating a system matrix; s3, substituting the system matrix into a preset frequency point to solve a Volterra series response matrix of the system, and performing singular value decomposition on the response matrix to generate an orthogonal temperature rise field base and a nonlinear auxiliary field base lambda s (r); s4, decoupling a nonlinear term mu of the temperature rise field by using an auxiliary field base lambda s (r), and realizing time domain solution of the transient nonlinear temperature rise field through an alternating iteration mechanism of the temperature rise field and an auxiliary field in a time domain; and S5, decoupling the multi-heat-source power density function into independent mapping through temperature rise field base superposition, and realizing lightweight multi-heat-source thermal model construction of the integrated microsystem.

Description

Technical Field

[0001] This invention belongs to the field of thermal management technology for microelectronic systems, particularly the construction of lightweight multi-heat source thermal models for integrated microsystems. It is especially applicable to transient thermal simulation and thermal management design technology for microscale systems such as chip packaging, three-dimensional integrated circuits, and MEMS devices. Specifically, it is a method and system for constructing lightweight multi-heat source thermal models for integrated microsystems. Background Technology

[0002] As microelectronic devices continue to evolve towards advanced process nodes, integrated microsystems face increasingly severe thermal management challenges. The local heat flux density of three-dimensional stacked chips increases significantly, leading to the risk of thermal failure; the dynamic time-varying thermal field of multi-core processors and heterogeneous integrated modules causes significant deviations in traditional single-heat-source models; and the nonlinear thermal conductivity characteristics of phase change materials and nanocomposite materials further increase the complexity of thermal field modeling.

[0003] Existing technologies have three main limitations:

[0004] While the full-order finite element method can guarantee accuracy, it consumes enormous computational resources. Transient simulations of million-degree-of-freedom models are too time-consuming, and the system matrix storage requirements far exceed the load capacity of ordinary workstations, making it difficult to meet the real-time requirements of design iterations.

[0005] While traditional reduced-order models can compress scale, their linear reduced-order basis cannot accurately capture the nonlinear effects of temperature-dependent thermal conductivity, leading to significant errors in strongly nonlinear scenarios. Furthermore, the cross-interference information from multiple heat sources carried by the basis functions can cause local temperature prediction biases.

[0006] While the equivalent thermal path method simplifies calculations, it neglects the spatial thermal gradient characteristics due to oversimplification of the physical model, resulting in transient response distortion. Furthermore, the model parameters rely on empirical calibration, which limits its universality.

[0007] Therefore, current technologies cannot simultaneously achieve efficiency, accuracy, and nonlinear / multi-heat-source characterization capabilities. There is an urgent need in this field to develop a nonlinear reduced-order basis generation method that supports temperature-dependent thermal conductivity, enabling independent mapping of multiple heat sources and decoupling of coupling effects, thus constructing a lightweight thermal model to meet real-time simulation requirements. Summary of the Invention

[0008] To address the aforementioned problems in existing technologies, this invention provides a lightweight multi-heat-source thermal model construction method and system for integrated microsystems. This invention achieves efficient and accurate simulation of complex thermal fields through three core technologies: nonlinear reduced-order basis generation, independent mapping of multiple heat sources, and dual-field iterative solution. First, this invention extracts the topological relationships and node coordinate data of tetrahedral elements using adaptive mesh analysis technology to construct a nonlinear material property correlation matrix. Then, it solves the Volterra series response within a preset frequency domain, generating orthogonalized temperature rise field basis and nonlinear auxiliary field basis. A transient nonlinear solution is completed using an alternating iterative algorithm of the temperature rise field and auxiliary field. Finally, the power density function of multiple heat sources is decoupled into independent mappings by superimposing the temperature rise field basis, thus realizing the construction of a lightweight multi-heat-source thermal model for integrated microsystems. The entire process compresses a million-degree-of-freedom system into a lightweight model of a dozen or so orders, overcoming the technical bottlenecks of traditional methods in terms of computational efficiency and nonlinear characterization.

[0009] The present invention adopts the following technical solution:

[0010] A method for constructing a lightweight multi-heat-source thermal model for integrated microsystems, comprising the following steps:

[0011] S1. By adaptively analyzing the mesh file, extract the tetrahedral elements, nodal coordinates, and surface element topological relationships to construct the nonlinear material property correlation matrix M_mat;

[0012] S2. Set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble element matrices and integrate system matrices based on the M_mat matrix;

[0013] S3, at the preset frequency point (10 -14 Up to 10 8 Substituting Hz into the system matrix, we solve for the Volterra series response of the system. Then, we perform singular value decomposition on the response dataset to generate orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r);

[0014] S4. Utilizing the nonlinear auxiliary field basis λ s (r) Decouple the nonlinear term μ of the temperature rise field and realize the time-domain solution of the transient nonlinear temperature rise field through the alternating iteration mechanism of the temperature rise field and the auxiliary field in the time domain;

[0015] S5, via the temperature rise field base θ j (r) Superposition decouples the power density function of multiple heat sources into independent mappings, realizing the construction of a lightweight multi-heat source thermal model for integrated microsystems.

[0016] Compared to traditional lightweight thermal model construction methods, this invention can characterize the influence of material nonlinearity on the thermal field, has a simple and high-precision construction process with multiple heat sources, and can realize intelligent compression of systems with tens of thousands of degrees of freedom to the order of tens of orders.

[0017] Preferably, in step S1, the mesh adaptive analysis and material modeling are as follows:

[0018] S1.1 Adaptively analyze the mesh file (.mphtxt) generated by COMSOL to accurately extract the topological relationship of tetrahedral elements (mesh_e), node coordinates (mesh_n), and surface element data (mesh_s);

[0019] S1.2 Construct the nonlinear material property correlation matrix M_mat, and label the domains d n It is correlated with material parameters {c,k0,μ,ρ}, where c is the constant-pressure heat capacity, k0 is the reference thermal conductivity, μ is the temperature coefficient of thermal conductivity, and ρ is the density. Furthermore, k0 supports anisotropy and can independently handle the thermal conductivity k in the x / y / z directions. 0x ,k 0y ,k 0z .

[0020] Preferably, step S2 is as follows:

[0021] S2.1 Independent time-varying power density configuration for each heat source:

[0022] P(t) = [P1(t), ..., P n (t)] (1)

[0023] Where t is the time variable, P(t) is the time-varying power density matrix of the total heat source, n is the number of heat sources, and P n (t) is the independent time-varying power density vector of the nth heat source;

[0024] S2.2. Perform heat source-domain mapping based on the spatial relationships of the tetrahedral elements extracted in S1.1:

[0025]

[0026] Where r is a spatial variable, i is the heat source number, and g i (r) is the spatial distribution matrix of the heat source, V i Let i be the volume (m) of heat source i 3 ), Ω i The scope of heat source i;

[0027] S2.3. Based on the heat transfer boundary surface settings and the surface element data `mesh_s` extracted in S1.1, the boundary surface type (internal / external) is automatically labeled, and the impact of the boundary layer effect on heat dissipation is compensated by equivalent area correction. Equivalent area of ​​the internal surface:

[0028] A eq =A (3)

[0029] Equivalent area of ​​the external surface:

[0030]

[0031] Where A is the geometric area. This is a correction factor for the equivalent heat dissipation area;

[0032] S2.4. Calculate the element stiffness matrix K based on the matrix M_mat constructed in S1.2. e Quality matrix C e and heat source matrix F e The system matrix is ​​constructed using a vectorized accumulation technique for efficient assembly: 4×4 unit matrices are added to their corresponding positions in the system matrix for batch processing; a pre-stored node index array `node_idx` avoids nested loops; `node_idx` ensures memory address continuity and improves cache hit rate; symmetry optimization is employed: anisotropic materials are automatically identified (when k...). 0x ≠k 0y (When symmetric optimization is disabled), only the upper triangular elements of the stiffness matrix are traversed, and symmetric assignment is performed on the off-diagonal elements of the lower triangular matrix.

[0033] Preferably, step S3 is as follows:

[0034] S3.1 Construction of the time-domain equation for nonlinear heat conduction:

[0035] The time-domain form of the nonlinear heat conduction equation is:

[0036]

[0037] Where r is a spatial variable, t is a time variable, ρ is density, c is constant pressure heat capacity, θ(r,t) is the temperature rise caused by power density q(r,t) relative to ambient temperature T0, and k is thermal conductivity.

[0038] The boundary conditions are described as follows:

[0039]

[0040] Where n(r) is the outward-pointing unit normal vector, and h(r) is the heat transfer coefficient, which is usually considered independent of temperature. In nonlinear heat conduction problems, the dependence of thermal conductivity on temperature cannot be ignored; therefore, the expression for thermal conductivity with respect to temperature is:

[0041]

[0042] Where k0 is the reference thermal conductivity and μ is the thermal conductivity temperature coefficient;

[0043] The nonlinearity of thermal conductivity is reformulated by introducing an auxiliary variable λ(r,t):

[0044]

[0045] Substituting equation (8) into equations (5) and (6):

[0046]

[0047] Since λ(r,t) is introduced, the equation now has two unknowns, so a new equation needs to be reconstructed for solving. By differentiating λ(r,t) with respect to time, we obtain:

[0048]

[0049] S3.2, Spatiotemporal Separation:

[0050] The temperature rise θ(r,t) can be approximated by a set of spatial basis functions as follows:

[0051]

[0052] Where, θ j (r) is the temperature rise space basis function. The number of basis functions. Let λ be the temperature rise time coefficient. Similarly, λ(r,t) can also be expressed as:

[0053]

[0054] Where, λ s (r) is the nonlinear auxiliary variable space basis function. The number of basis functions. The time coefficients are nonlinear auxiliary variables. For the input power density q(r,t), since the input power density depends only on the geometry, q(r,t) is a linear variable; therefore, the dimension of the basis functions is... q(r,t) is represented as:

[0055] q(r,t)=g(r)P(t) (14)

[0056] Where g(r) is the spatial distribution matrix of the heat source, and P(t) is the power density.

[0057] S3.3 Solving for the basis functions in the frequency domain using the Volterra kernel:

[0058] Introducing variables

[0059]

[0060] Rewrite equation (9) as:

[0061]

[0062] Rewrite equation (11) as:

[0063]

[0064] By analyzing the variables λ(r,t), θ(r,t), and Performing a multivariable Laplace transform, we obtain the Volterra series expansion term Λ. m (r,s m ), Θ m (r,s m ) and ψ m (r,s m )as follows:

[0065]

[0066] Among them, s m It is a complex frequency variable, each s k (k=1,…,m) all belong to the complex number domain and are used to characterize the response characteristics of the Volterra kernel function in the frequency domain. m p′ is an m-dimensional vector composed of 0s and 1s. m It is by p m The vector obtained by replacing all 0s with 1s and 1s with 0s is a vector. By selecting s m In and p m The vector obtained from the corresponding elements; δ m1 The value of is 1 only when m = 1, and 0 otherwise. It is used to control the activation condition of g(r) to ensure that the power excitation only acts on the first-order system.

[0067] Around the selected frequency point σ m Expanding the Volterra series, the preset frequency range is σ. m ∈[10 -14 10 8]Hz, covering low-frequency thermal diffusion to high-frequency transient response; frequency points are logarithmically uniformly distributed, dense in the low-frequency region and sparse in the high-frequency region, Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m This can be represented as:

[0068]

[0069] Where, α m For multi-indexed vectors, describe the moments of the Volterra series expansion, where α is the moment when m = 1. m The value of α is 0; when m = 2, α m The value of α is (0,0); when m = r, α m The value of is (0,0,…0), representing an r-dimensional zero vector; For Θ m (r,s m ) moment;

[0070] Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m Substituting the Taylor series of ) into equations (18)-(21), we get:

[0071]

[0072] Where ρ(r)c(r) is the system quality matrix C constructed in S2.4. That is, the system stiffness matrix K constructed in S2.4, g(r) is the heat source matrix F constructed in S2.4; e m,k For the k-th unit vector, The value is only when α m = 1 when = 0, otherwise = 0;

[0073] By solving the linear heat diffusion problem in the frequency domain (Equations (25) and (26)), the temperature rise response matrix Θ(r,σ) is obtained. m ) and nonlinear auxiliary variable matrix Λ(r,σ m Basis function θ j (r) and λ s (r) are respectively and The results of the standard orthogonalization.

[0074] Preferably, step S4 is as follows:

[0075] S4.1 Construction of Time-Domain Iterative Equations:

[0076] Multiply both sides of equation (5) by θ i (r), and integrate over the spatial domain Ω. Then, take equations (6), (9), (10), and (11) for all

[0077]

[0078] in,

[0079]

[0080] Where θ(r) is the basis function Θ(r) obtained in S3.3;

[0081] Similarly, multiplying both sides of equation (11) by λ i (r), and integrate over the spatial domain Ω, substituting equations (9) and (10) into the equations, for all

[0082]

[0083] in,

[0084]

[0085] Where λ(r) is the basis function Λ(r) obtained in S3.3;

[0086] S4.2, Time Discreteness:

[0087] In the time-domain solution, the implicit Euler method is used for time discretization, transforming the differential terms into difference forms; let the time step be Δt, then...

[0088]

[0089] Where n is the time step index;

[0090] The time discretization formula of equation (29) is as follows:

[0091]

[0092] Similarly, the formula after discretizing equation (31) in time is:

[0093]

[0094] S4.3, Two-field iterative solution:

[0095] The temperature rise field and auxiliary field of the nonlinear model are solved iteratively in the time domain. Equation (34) is then rearranged to obtain:

[0096]

[0097] Where m is the iteration step index;

[0098] Simplifying equation (35) yields:

[0099]

[0100] Using the auxiliary field freezing technique, in each time step iteration, the nonlinear auxiliary variable is... Treating the variables as known quantities (initially 0), nonlinear terms are decoupled. The old and new iterative residuals of the temperature rise field and auxiliary field are monitored synchronously until convergence. Then, the iterative solution is performed for the next time step, finally completing the time-domain solution of the transient nonlinear temperature rise field and obtaining the temperature rise time coefficient. and auxiliary variable time coefficient

[0101] Preferably, step S5 is as follows:

[0102] Construction of a multi-heat-source thermal model:

[0103] A basis function can only carry information about one heat source. The θ obtained in S3.3... j (r) rewritten as θ ij (r) achieves the synthesis of multi-heat-source temperature rise fields through the linear superposition of the results of multiplying multiple temperature rise spatial basis functions and temperature rise time coefficients:

[0104]

[0105] in, For the complete temperature rise variable, θ ij (r) is the temperature rise basis function corresponding to the i-th heat source, and n is the number of heat sources;

[0106] The complete temperature field T(r,t) of the lightweight multi-heat-source thermal model of the final integrated microsystem is calculated by the following formula:

[0107]

[0108] This invention also discloses a lightweight multi-heat-source thermal model construction system for integrated microsystems, used to perform the above method, comprising the following modules:

[0109] The nonlinear material property correlation matrix construction module is used to extract tetrahedral elements, nodal coordinates, and surface element topological relationships from adaptively parsed mesh files, and construct the nonlinear material property correlation matrix M_mat.

[0110] System matrix integration module: used to set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble element matrices and integrate system matrices based on the M_mat matrix;

[0111] Field basis generation module: Used to solve the Volterra series response matrix of the system by substituting the system matrix at a preset frequency point, perform singular value decomposition on the response matrix, and generate orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r);

[0112] Transient nonlinear temperature rise field time-domain solution module: utilizing nonlinear auxiliary field basis λ s (r) Decouple the nonlinear term μ of the temperature rise field and realize the time-domain solution of the transient nonlinear temperature rise field through the alternating iteration mechanism of the temperature rise field and the auxiliary field in the time domain;

[0113] Lightweight multi-heat-source thermal model building module: used to construct thermal models based on temperature rise field θ j (r) Superposition decouples the power density function of multiple heat sources into independent mappings, realizing the construction of a lightweight multi-heat source thermal model for integrated microsystems.

[0114] Compared to traditional lightweight thermal model construction methods, this invention can characterize the influence of material nonlinearity on the thermal field, has a simple and high-precision construction process with multiple heat sources, and can realize intelligent compression of systems with tens of thousands of degrees of freedom to the order of tens of orders. Attached Figure Description

[0115] Figure 1 This is a diagram of the actual measured system of the thermal testing chip microsystem;

[0116] Figure 2 This is a 3D model diagram of the thermal testing chip microsystem;

[0117] Figure 3 This is a flowchart illustrating the preferred embodiment of the present invention for constructing a lightweight multi-heat-source thermal model for an integrated microsystem.

[0118] Figure 4 It is a node temperature map of the test, simulation (FEM), and thermal model (ROM) of the thermal test chip microsystem;

[0119] Figure 5 This is a block diagram of a lightweight multi-heat-source thermal model construction system for integrated microsystems, according to a preferred embodiment of the present invention. Detailed Implementation

[0120] The following specific implementation examples illustrate the embodiments of the present invention. Those skilled in the art can intuitively understand other advantages of the present invention from the content described in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features in the embodiments can be combined with each other.

[0121] The RQS-TTC-SC01 thermal test chip is used as an example model below. This chip accurately simulates the heating characteristics of an actual working chip through a heating element, and uses a temperature measuring element to monitor the surface temperature of the chip under transient temperature changes of the heating element. Three thermal test chips are attached to the PCB board with thermally conductive adhesive to simulate the situation where multiple heat sources generate heat simultaneously. Detailed information on the test chip and PCB board can be found in [reference needed]. Figure 1 .

[0122] Reference Figure 2 The 3D models of the thermal test PCB and the bonded chip were constructed using the commercial software COMSOL.

[0123] Reference Figure 3 This implementation case illustrates the specific steps of a lightweight multi-heat-source thermal model construction method for integrated microsystems:

[0124] S1. Mesh Analysis and Material Matrix Construction: Through adaptive analysis of the COMSOL mesh file (with 23454 degrees of freedom), the topological relationships of tetrahedral elements (mesh_e), node coordinates (mesh_n), and surface element data (mesh_s) are accurately extracted; a nonlinear material property correlation matrix M_mat is constructed, with domain identifiers d... n It is correlated with material parameters {c,k0,μ,ρ}, where c is the constant-pressure heat capacity, k0 is the reference thermal conductivity, μ is the temperature coefficient of thermal conductivity, and ρ is the density. Furthermore, k0 supports anisotropy, and the thermal conductivity k in the x / y / z directions is handled independently. 0x ,k 0y ,k 0z The material parameters for this embodiment are shown in Table 1, which includes silicon, a material with nonlinear thermal conductivity, and FR4, a material with anisotropic thermal conductivity.

[0125] Table 1 Material parameters of the thermal testing chip microsystem

[0126]

[0127] S2. Set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble the element matrix and integrate the system matrix based on the M_mat matrix:

[0128] S2.1 Independent time-varying power density configuration for each heat source:

[0129]

[0130] S2.2. Perform heat source-domain mapping based on the spatial relationships of the tetrahedral elements extracted in S1.1:

[0131]

[0132] Among them, g i (r) is the spatial distribution matrix of the heat source, V i Let i be the volume (m) of heat source i 3 ), Ω i The domain of heat source i; Ω1, Ω2, and Ω3 respectively correspond to Figure 1 The regions of chip 1, chip 2, and chip 3;

[0133] S2.3 The heat transfer boundary surface is the outer surface of the entire model, and the heat transfer coefficient is 16 [W / m]. 2 ·K]. Based on the heat transfer boundary surface settings and the surface element data mesh_s extracted in S1.1, the boundary surface type (internal / external) is automatically labeled, and the impact of the boundary layer effect on heat dissipation is compensated by equivalent area correction. Equivalent area of ​​the internal surface:

[0134] A eq =A(3)

[0135] Equivalent area of ​​the external surface:

[0136]

[0137] Where A is the geometric area. This is a correction factor for the equivalent heat dissipation area;

[0138] S2.4 Construct the element stiffness matrix K based on the matrix M_mat constructed in step S1.2. e Quality matrix C e and heat source matrix F e Vectorized accumulation technology is used to efficiently assemble the global matrices of element stiffness matrix, mass matrix, and heat source matrix to obtain system matrices K, C, and F. The 4×4 element matrices are added to the corresponding positions in the global matrix for batch processing. A pre-stored node index array `node_idx` avoids nested loops, and `node_idx` ensures memory address continuity, improving cache hit rate. Symmetric optimization is employed: anisotropic materials are automatically identified (when k...). 0x ≠k 0y (When symmetric optimization is disabled), only the upper triangular elements of the stiffness matrix are traversed, and symmetric assignment is performed on the off-diagonal elements of the lower triangular matrix.

[0139] S3, at the preset frequency point (10 -14 Up to 10 8 Substituting Hz into the system matrix, we solve for the Volterra series response of the system. Then, we perform singular value decomposition on the response dataset to generate orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r). Details are as follows:

[0140] S3.1 Construction of the time-domain equation for nonlinear heat conduction:

[0141] The time-domain form of the nonlinear heat conduction equation is:

[0142]

[0143] Where r is a spatial variable, t is a time variable, ρ is density, c is constant pressure heat capacity, θ(r,t) is the temperature rise caused by power density q(r,t) relative to ambient temperature T0, and k is thermal conductivity.

[0144] The boundary conditions are described as follows:

[0145]

[0146] Where n(r) is the outward-pointing unit normal vector, and h(r) is the heat transfer coefficient; the expression for thermal conductivity with respect to temperature is:

[0147]

[0148] Where k0 is the reference thermal conductivity and μ is the thermal conductivity temperature coefficient;

[0149] The nonlinearity of thermal conductivity is reformulated by introducing an auxiliary variable λ(r,t):

[0150]

[0151] Substituting equation (8) into equations (5) and (6):

[0152]

[0153] By taking the derivative of λ(r,t) with respect to time, we obtain:

[0154]

[0155] S3.2, Spatiotemporal Separation:

[0156] The temperature rise θ(r,t) can be approximated by a set of spatial basis functions as follows:

[0157]

[0158] Where, θ j (r) is the temperature rise space basis function. The number of basis functions. The temperature rise time coefficient; λ(r,t) can also be expressed as:

[0159]

[0160] Where, λ s (r) is the nonlinear auxiliary variable space basis function. The number of basis functions. The time coefficients are nonlinear auxiliary variables; for the input power density q(r,t), since the input power density depends only on the geometry, q(r,t) is a linear variable, therefore the dimension of the basis functions is... q(r,t) is represented as:

[0161] q(r,t)=g(r)P(t) (14)

[0162] Where g(r) is the spatial distribution matrix of the heat source, and P(t) is the power density;

[0163] S3.3 Solving for the basis functions in the frequency domain using the Volterra kernel:

[0164] Introducing variables

[0165]

[0166] Rewrite equation (9) as:

[0167]

[0168] Rewrite equation (11) as:

[0169]

[0170] By analyzing the variables λ(r,t), θ(r,t), and Performing a multivariable Laplace transform, we obtain the Volterra series expansion term Λ. m (r,s m ), Θ m (r,s m ) and ψ m (r,s m )as follows:

[0171]

[0172]

[0173] Among them, s mIt is a complex frequency variable, each s k All belong to the complex number domain, k = 1, ..., m, and are used to characterize the frequency domain response of the Volterra kernel function; p m p′ is an m-dimensional vector composed of 0s and 1s. m It is by p m The vector obtained by replacing all 0s with 1s and 1s with 0s is a vector. By selecting s m In and p m The vector obtained from the corresponding elements; δ m1 The value of is 1 only when m = 1, and 0 otherwise, and is used to control the activation condition of g(r);

[0174] Around the selected frequency point σ m Expanding the Volterra series, the preset frequency range is σ. m ∈[10 -14 10 8 ]Hz; Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m ) is represented as:

[0175]

[0176] Where, α m For multi-indexed vectors, describe the moments of the Volterra series expansion, where α is the moment when m = 1. m The value of α is 0; when m = 2, α m The value of α is (0,0); when m = r, α m The value of is (0,0,…0), representing an r-dimensional zero vector; For Θ m (r,s m ) moment;

[0177] Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m Substituting the Taylor series of ) into equations (18)-(21), we get:

[0178]

[0179] Where ρ(r)c(r) is the system quality matrix C constructed in step S2.4. Let K be the system stiffness matrix constructed in step S2.4, and g(r) be the heat source matrix F constructed in step 2.4; m,k For the k-th unit vector, The value is only when α m = 1 when = 0, otherwise = 0;

[0180] By solving the linear heat diffusion problem in the frequency domain, i.e., equations (25) and (26), the temperature rise response matrix Θ(r,σ) is obtained. m ) and nonlinear auxiliary variable matrix Λ(r,σ m ); basis function θ j (r) and λ s (r) are respectively and The results of the standard orthogonalization.

[0181] S4. Utilizing the nonlinear auxiliary field basis λ s (r) The nonlinear term μ of the temperature rise field is decoupled, and the transient nonlinear temperature rise field is solved in the time domain through an alternating iterative mechanism of the temperature rise field and the auxiliary field in the time domain. Specifically:

[0182] S4.1 Construction of Time-Domain Iterative Equations:

[0183] Multiply both sides of equation (5) by θ i (r), and integrate over the spatial domain Ω; then, take equations (6), (9), (10), and (11) for all

[0184]

[0185] in,

[0186]

[0187] Where θ(r) is the basis function Θ(r) obtained in step S3.3;

[0188] Multiply both sides of equation (11) by λ i (r), and integrate over the spatial domain Ω, substituting equations (9) and (10) into the equations, for all

[0189]

[0190] in,

[0191]

[0192] Where λ(r) is the basis function Λ(r) obtained in step S3.3;

[0193] S4.2, Time Discreteness:

[0194] In the time-domain solution, the implicit Euler method is used for time discretization, transforming the differential terms into difference forms; let the time step be Δt, then:

[0195]

[0196] Where n is the time step index;

[0197] The time discretization formula of equation (29) is as follows:

[0198]

[0199] Equation (31) after time discretization is as follows:

[0200]

[0201] S4.3, Two-field iterative solution:

[0202] The temperature rise field and auxiliary field of the nonlinear model are solved iteratively in the time domain. Equation (34) is then rearranged to obtain:

[0203]

[0204] Where m is the iteration step index;

[0205] Simplifying equation (35) yields:

[0206]

[0207] Using the auxiliary field freezing technique, in each time step iteration, the nonlinear auxiliary variable is... Treating the variables as known quantities (initially 0), nonlinear terms are decoupled. The old and new iterative residuals of the temperature rise field and auxiliary field are monitored synchronously until convergence. Then, the iterative solution is performed for the next time step, finally completing the time-domain solution of the transient nonlinear temperature rise field and obtaining the temperature rise time coefficient. and auxiliary variable time coefficient

[0208] In summary, this step iteratively solves the temperature rise field and auxiliary field of the nonlinear model in the time domain, employing an auxiliary field freezing technique. In each iteration, the nonlinear auxiliary field is treated as a known quantity (using the value of the current iteration step), thus decoupling the nonlinear terms. The new and old iteration residuals of the temperature rise field and auxiliary field are synchronously monitored until convergence (residual less than 10). -12 Finally, the time-domain solution of the transient nonlinear temperature rise field is completed, and the temperature rise time coefficient is obtained. and auxiliary variable time coefficient

[0209] S5, via the temperature rise field base θ j(r) Superposition decouples the power density functions of multiple heat sources into independent mappings, enabling the construction of a lightweight multi-heat source thermal model for integrated microsystems. Specifically:

[0210] Multi-heat-source thermal model construction achieved by superimposing basis functions:

[0211] A basis function can only carry information about one heat source. The θ obtained from S3... j (r) rewritten as θ ij (r) achieves the synthesis of multi-heat-source temperature rise fields through the linear superposition of the results of multiplying three temperature rise spatial basis functions and the temperature rise time coefficient:

[0212]

[0213] The complete temperature field T(r,t) of the lightweight multi-heat-source thermal model of the final integrated microsystem is calculated by the following formula:

[0214]

[0215] The initial temperature T0 is 300K;

[0216] The dual-field iterative solution took only 0.083025 seconds, and the lightweight multi-heat-source thermal model construction took 0.0882011 seconds.

[0217] Reference Figure 4 The figure shows the node temperatures for the test, simulation (FEM), and thermal model (ROM). Figure 4 In the figure, (a) shows the curves of power applied to the three chips; (b)-(d) are the temperature curves of Chip1, Chip2, and Chip3, respectively. The results in the figure show that the thermal model construction method of the present invention can well characterize the actual thermal situation and finite element calculation results of the integrated microsystem of nonlinear multi-heat source of the material.

[0218] like Figure 5 As shown, this embodiment discloses a lightweight multi-heat-source thermal model construction system for integrated microsystems, used to execute the above method, which includes the following modules:

[0219] The nonlinear material property correlation matrix construction module is used to extract tetrahedral elements, nodal coordinates, and surface element topological relationships from adaptively parsed mesh files, and construct the nonlinear material property correlation matrix M_mat.

[0220] System matrix integration module: used to set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble element matrices and integrate system matrices based on the M_mat matrix;

[0221] Field basis generation module: Used to solve the Volterra series response matrix of the system by substituting the system matrix at a preset frequency point, perform singular value decomposition on the response matrix, and generate orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r);

[0222] Transient nonlinear temperature rise field time-domain solution module: utilizing nonlinear auxiliary field basis λ s (r) Decouple the nonlinear term μ of the temperature rise field and realize the time-domain solution of the transient nonlinear temperature rise field through the alternating iteration mechanism of the temperature rise field and the auxiliary field in the time domain;

[0223] Lightweight multi-heat-source thermal model building module: used to construct thermal models based on temperature rise field θ j (r) Superposition decouples the power density function of multiple heat sources into independent mappings, realizing the construction of a lightweight multi-heat source thermal model for integrated microsystems.

[0224] Other aspects of this embodiment can be found in the above method embodiments.

[0225] In summary, this invention is applied to the field of thermal management technology for microelectronic systems, providing a lightweight multi-heat-source thermal model construction method and system for integrated microsystems. This invention overcomes the technical bottlenecks of traditional thermal simulation through three core technologies: First, it employs a nonlinear reduced-order basis intelligent generation technology, scanning the thermal diffusion response in a wide frequency domain based on Volterra series, and generating temperature rise field basis and auxiliary field basis through orthogonalization processing, achieving a significant compression of system degrees of freedom and characterizing nonlinear thermal conductivity characteristics; second, it introduces a multi-heat-source independent mapping mechanism, eliminating cross-interference between heat sources through precise location of the heat source domain and basis function superposition, supporting independent power configuration; finally, it develops a dual-field iterative solution engine, reconstructing the nonlinear thermal conductivity into an equivalent linear form through auxiliary field freezing technology, achieving alternating iterative solution of the temperature rise field and auxiliary field.

[0226] Compared to the traditional full-order finite element method, the technical solution of this invention has higher computational efficiency and lower memory usage. Compared to the traditional reduced-order model, this invention has stronger characterization capabilities and higher accuracy in nonlinear thermal conductivity characterization and multi-heat source coupling modeling. In advanced process node integrated microsystem models, this solution exhibits excellent transient response accuracy. This invention overcomes the long-standing contradiction between efficiency and accuracy in the field of microsystem thermal simulation. Through lightweight modeling and nonlinear decoupling techniques, it provides a real-time, high-precision thermal model building tool for high-power-density integrated microsystems, promoting the evolution of microelectronic thermal management technology towards intelligence and efficiency.

[0227] The above description is merely a preferred embodiment of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments, and substitutions can be made without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and may include many other equivalent embodiments without departing from the concept of the present invention, the scope of which is determined by the scope of the appended claims.

Claims

1. A method for constructing a lightweight multi-heat-source thermal model for integrated microsystems, characterized in that, Includes the following steps: S1. By adaptively analyzing the mesh file, extract the tetrahedral elements, nodal coordinates, and surface element topological relationships to construct the nonlinear material property correlation matrix M_mat; S2. Set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble element matrices and integrate system matrices based on the M_mat matrix; S3. Substitute the system matrix into the preset frequency point to solve for the Volterra series response matrix of the system. Perform singular value decomposition on the response matrix to generate the orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r); S4. Utilizing the nonlinear auxiliary field basis λ s (r) Decouple the nonlinear term μ of the temperature rise field and realize the time-domain solution of the transient nonlinear temperature rise field through the alternating iteration mechanism of the temperature rise field and the auxiliary field in the time domain; S5, via the temperature rise field base θ j (r) Superposition decouples the power density function of multiple heat sources into independent mappings, realizing the construction of a lightweight multi-heat source thermal model for integrated microsystems.

2. The method for constructing a lightweight multi-heat-source thermal model for integrated microsystems according to claim 1, characterized in that, Step S1 specifically includes: S1.1 Adaptively analyze the COMSOL-generated mesh file .mphtxt to extract the tetrahedral element topology mesh_e, node coordinates mesh_n, and surface element data mesh_s; S1.2 Construct the nonlinear material property correlation matrix M_mat, and label the domains d n It is correlated with the material parameters {c,k0,μ,ρ}, where c is the constant pressure heat capacity, k0 is the reference thermal conductivity, μ is the thermal conductivity temperature coefficient, and ρ is the density.

3. The method for constructing a lightweight multi-heat-source thermal model for integrated microsystems according to claim 2, characterized in that, Step S2 specifically includes: S2.1 Independent time-varying power density configuration for each heat source: P(t)=[P1(t),…,P n (t)] (1) Where t is the time variable, P(t) is the time-varying power density matrix of the total heat source, n is the number of heat sources, and P n (t) is the independent time-varying power density vector of the nth heat source; S2.

2. Perform heat source-domain mapping based on the spatial relationships of the tetrahedral elements extracted in S1.1: Where r is a spatial variable, i is the heat source number, and g i (r) is the spatial distribution matrix of the heat source, V i Let i be the volume (m) of heat source i 3 ), Ω i The scope of heat source i; S2.

3. Based on the heat transfer boundary surface settings and the surface element data mesh_s extracted in step S1.1, the boundary surface type is automatically marked, and the impact of the boundary layer effect on heat dissipation is compensated by equivalent area correction; the equivalent area of ​​the internal surface: A eq =A (3) Equivalent area of ​​the external surface: Where A is the geometric area. This is a correction factor for the equivalent heat dissipation area; S2.

4. Based on the matrix M_mat constructed in step S1.2, calculate the element stiffness matrix K. e Quality matrix C e and heat source matrix F e The unit stiffness matrix, mass matrix, and heat source matrix are constructed by using vectorized accumulation technology to assemble them into system matrices K, C, and F, respectively.

4. The method for constructing a lightweight multi-heat-source thermal model for integrated microsystems according to claim 3, characterized in that, Step S3 specifically includes: S3.1 Construction of the time-domain equation for nonlinear heat conduction: The time-domain form of the nonlinear heat conduction equation is: Where r is a spatial variable, t is a time variable, ρ is density, c is constant pressure heat capacity, θ(r,t) is the temperature rise caused by power density q(r,t) relative to ambient temperature T0, and k is thermal conductivity. The boundary conditions are described as follows: Where n(r) is the outward-pointing unit normal vector, and h(r) is the heat transfer coefficient; the expression for thermal conductivity with respect to temperature is: Where k0 is the reference thermal conductivity and μ is the thermal conductivity temperature coefficient; The nonlinearity of thermal conductivity is reformulated by introducing an auxiliary variable λ(r,t): Substituting equation (8) into equations (5) and (6): By taking the derivative of λ(r,t) with respect to time, we obtain: S3.2, Spatiotemporal Separation: The temperature rise θ(r,t) can be approximated by a set of spatial basis functions as follows: Where, θ j (r) is the space basis function for temperature rise. The number of basis functions. The temperature rise time coefficient; λ(r,t) can also be expressed as: Where, λ s (r) is the nonlinear auxiliary variable space basis function. The number of basis functions. The time coefficients are nonlinear auxiliary variables; for the input power density q(r,t), since the input power density depends only on the geometry, q(r,t) is a linear variable, therefore the dimension of the basis functions is... q(r,t) is represented as: q(r,t)=g(r)P(t) (14) Where g(r) is the spatial distribution matrix of the heat source, and P(t) is the power density; S3.3 Solving for the basis functions in the frequency domain using the Volterra kernel: Introducing variables Rewrite equation (9) as: Rewrite equation (11) as: By analyzing the variables λ(r,t), θ(r,t), and Performing a multivariable Laplace transform, we obtain the Volterra series expansion term Λ. m (r,s m ), Θ m (r,s m ) and ψ m (r,s m )as follows: Among them, s m It is a complex frequency variable, each s k All belong to the complex number domain, k = 1, ..., m, and are used to characterize the frequency domain response of the Volterra kernel function; p m p is an m-dimensional vector composed of 0s and 1s. ′ m It is by p m The vector obtained by replacing all 0s with 1s and 1s with 0s is a vector. By selecting s m In and p m The vector obtained from the corresponding elements; δ m1 The value of is 1 only when m = 1, and 0 otherwise, and is used to control the activation condition of g(r); Around the selected frequency point σ m Expanding the Volterra series, the preset frequency range is σ. m ∈[10 -14 10 8 ]Hz; Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m ) is represented as: Where, α m For multi-indexed vectors, describe the moments of the Volterra series expansion, where α is the moment when m = 1. m The value of α is 0; when m = 2, α m The value of α is (0,0); when m = r, α m The value of is (0,0,…0), representing an r-dimensional zero vector; For Θ m (r,s m ) moment; Θ m (r,s m ), ψ m (r,s m ) and Λ m (r,s m Substituting the Taylor series of ) into equations (18)-(21), we get: Where ρ(r)c(r) is the system quality matrix C constructed in step S2.

4. The system stiffness matrix K is constructed in step S2.4, and g(r) is the heat source matrix F constructed in step S2.4; m,k For the k-th unit vector, The value is only when α m = 1 when = 0, otherwise = 0; By solving the linear heat diffusion problem in the frequency domain, i.e., equations (25) and (26), the temperature rise response matrix Θ(r,σ) is obtained. m ) and nonlinear auxiliary variable matrix Λ(r,σ m ); basis function θ j (r) and λ s (r) are respectively and The results of the standard orthogonalization.

5. The method for constructing a lightweight multi-heat-source thermal model for integrated microsystems according to claim 4, characterized in that, Step S4 is as follows: S4.1, Construction of time-domain iterative equations: Multiply both sides of equation (5) by θ i (r), and integrate over the spatial domain Ω; then, take equations (6), (9), (10), and (11) for all in, Where θ(r) is the basis function Θ(r) obtained in step S3.3; Multiply both sides of equation (11) by λ i (r), and integrate over the spatial domain Ω, substituting equations (9) and (10) into the equations, for all in, Where λ(r) is the basis function Λ(r) obtained in step S3.3; S4.2, Time Discreteness: In the time-domain solution, the implicit Euler method is used for time discretization, transforming the differential terms into difference forms; let the time step be Δt, then: Where n is the time step index; The time discretization formula of equation (29) is as follows: Equation (31) after time discretization is as follows: S4.3, Two-field iterative solution: The temperature rise field and auxiliary field of the nonlinear model are solved iteratively in the time domain. Equation (34) is then rearranged to obtain: Where m is the iteration step index; Simplifying equation (35) yields: Using the auxiliary field freezing technique, in each time step iteration, the nonlinear auxiliary variable is... Treating these as known quantities, nonlinear terms are decoupled. The old and new iterative residuals of the temperature rise field and auxiliary field are monitored synchronously until convergence. Then, the iterative solution is performed for the next time step, finally completing the time-domain solution of the transient nonlinear temperature rise field and obtaining the temperature rise time coefficient. and auxiliary variable time coefficient 6. The method for constructing a lightweight multi-heat-source thermal model for integrated microsystems according to claim 5, characterized in that, Step S5 is as follows: Construction of a multi-heat-source thermal model: The θ obtained in step S3.3... j (r) rewritten as θ ij (r) achieves the synthesis of multi-heat-source temperature rise fields through the linear superposition of the results of multiplying multiple temperature rise spatial basis functions and temperature rise time coefficients: in, For the complete temperature rise variable, θ ij (r) is the temperature rise basis function corresponding to the i-th heat source, and n is the number of heat sources; The complete temperature field T(r,t) of the lightweight multi-heat-source thermal model of the final integrated microsystem is calculated by the following formula:

7. A lightweight multi-heat-source thermal modeling system for integrated microsystems, used to perform the method as described in any one of claims 1-6, characterized in that, Includes the following modules: The nonlinear material property correlation matrix construction module is used to extract tetrahedral elements, nodal coordinates, and surface element topological relationships from adaptively parsed mesh files, and construct the nonlinear material property correlation matrix M_mat. System matrix integration module: used to set transient heat sources and heat transfer boundaries according to the spatial relationship of tetrahedral elements, assemble element matrices and integrate system matrices based on the M_mat matrix; Field basis generation module: Used to solve the Volterra series response matrix of the system by substituting the system matrix at a preset frequency point, perform singular value decomposition on the response matrix, and generate orthogonal temperature rise field basis θ. j (r) and the nonlinear auxiliary field basis λ s (r); Transient nonlinear temperature rise field time-domain solution module: utilizing nonlinear auxiliary field basis λ s (r) Decouple the nonlinear term μ of the temperature rise field and realize the time-domain solution of the transient nonlinear temperature field through the alternating iteration mechanism of the temperature field and the auxiliary field in the time domain; Lightweight multi-heat-source thermal model building module: used to construct thermal models based on temperature rise field θ j (r) Superposition decouples the power density function of multiple heat sources into independent mappings, realizing the construction of a lightweight multi-heat source thermal model for integrated microsystems.

Images