Three-dimensional integrated micro-system electro-thermal synergic calculation method based on meshless galerkin method

Abstract

The application discloses a three-dimensional integrated microsystem electric-thermal cooperative calculation method based on a meshless Galerkin method, and comprises the following steps: establishing a background integral grid based on a three-dimensional integrated microsystem model, and establishing integral points by using grid and boundary information; finding integral points in a solution domain, and determining field nodes of a support domain where the integral points are located; calculating a shape function and a derivative thereof at each integral point; obtaining a temperature stiffness matrix and a node thermal load matrix at a boundary, and a basic temperature stiffness matrix and a basic node thermal load matrix based on the shape function and the derivative thereof; adding electric-thermal coupling terms in the basic temperature stiffness matrix and the basic node thermal load matrix to obtain a temperature stiffness matrix and a node thermal load matrix in a discrete domain; and obtaining a temperature field of the three-dimensional integrated microsystem model according to the temperature stiffness matrix and the node thermal load matrix at the boundary, and the temperature stiffness matrix and the node thermal load matrix in the discrete domain. The application eliminates the influence of grid shape on a solution result, and ensures precision.

Description

Technical Field

[0001] This invention belongs to the field of microelectronics and three-dimensional thermal analysis numerical calculation, specifically involving a three-dimensional integrated microsystem electrothermal co-calculation method based on the meshless Galerkin method. Background Technology

[0002] The development of microsystem integration technology has brought many disruptive advancements to the performance of electronic products, offering advantages such as high integration density, short interconnects, high reliability, and low power consumption. However, with the advancement of integration technology, the development of three-dimensional integrated microsystems faces several challenges. For example, multi-chip stacking leads to increased power consumption and reduced cooling efficiency, resulting in higher temperatures within the integrated circuits. As temperature changes, the thermal and electrical conductivity of materials alter simultaneously during heat transfer, affecting hotspot temperatures. Localized hotspots pose numerous challenges to electrical design aspects such as leakage current, electromigration, and signal and power integrity, significantly reducing chip performance, reliability, and lifespan. Furthermore, with continuously increasing integration density, thermal problems become even more severe. In conclusion, conducting electrothermal coupling analysis in the early stages of the design of high-power three-dimensional integrated microsystems plays a crucial role in ensuring system reliability.

[0003] Furthermore, 3D integrated microsystems encompass a wide range of geometric scales, from nanometers to millimeters, from the packaging substrate to the RDL layer. Common finite element analysis methods for meshing these microsystems typically employ two approaches: 1. reducing the mesh size; 2. using irregular shapes to adapt to the multi-scale geometry. The former results in dense mesh generation, increasing the computational load of multiphysics coupling calculations and raising time costs; the latter affects the accuracy of electrothermal coupling calculations, leading to significant errors.

[0004] Therefore, finding a method for electrothermal co-computation that adapts to multiple geometric scales in three-dimensional integrated microsystems has become a key focus. Summary of the Invention

[0005] To address the aforementioned problems in the existing technology, this invention provides a method for electrothermal co-calculation of three-dimensional integrated microsystems based on the meshless Galerkin method.

[0006] The technical problem to be solved by this invention is achieved through the following technical solution:

[0007] A method for electrothermal co-computation of a three-dimensional integrated microsystem based on the meshless Galerkin method, the method comprising:

[0008] Step 1: Obtain the field node information of the three-dimensional integrated microsystem model;

[0009] Step 2: Establish a background integral mesh based on the three-dimensional integrated microsystem model, establish integration points using mesh and boundary information, and calculate the information of the integration points within the background mesh. The information of the integration points includes the coordinates and weights of the integration points.

[0010] Step 3: Locate the integration points within the solution domain and determine the field nodes within the support domain of each integration point;

[0011] Step 4: Calculate the shape function and its derivative at each integration point;

[0012] Step 5: Based on the shape function and its derivative, calculate the temperature stiffness matrix and nodal thermal load matrix at the boundary, and calculate the contribution of each integration point to the temperature stiffness matrix and nodal load matrix to obtain the basic temperature stiffness matrix and basic nodal thermal load matrix.

[0013] Step 6: Based on the influence of temperature on the three-dimensional integrated microsystem model, add electrothermal coupling terms to the basic temperature stiffness matrix and the basic nodal thermal load matrix to obtain the temperature stiffness matrix and nodal thermal load matrix in the discrete domain.

[0014] Step 7: Based on the heat conduction control equation, obtain the temperature field of the three-dimensional integrated microsystem model according to the temperature stiffness matrix and nodal thermal load matrix at the boundary, as well as the temperature stiffness matrix and nodal thermal load matrix in the discrete domain.

[0015] Optionally, step 1 includes:

[0016] Step 1.1: Obtain the three-dimensional integrated microsystem model;

[0017] Step 1.2: Discretize the field nodes of the three-dimensional integrated microsystem model to obtain the field node information, which includes: the total number of field nodes, the field node number, the field node coordinates, and the field node spacing.

[0018] Optionally, step 2 includes:

[0019] Step 2.1: Establish the background integral mesh based on the three-dimensional integrated microsystem model, wherein the background integral mesh covers the three-dimensional integrated microsystem model, and the background integral mesh includes several mesh points arranged in a matrix form;

[0020] Step 2.2: Within each background grid, based on the number of integration points in one dimension, the weight of the integration point is determined using a Gaussian function, and the coordinates of the integration point are obtained by multiplying the distance between two adjacent grid points in a background grid with the weight of the integration point.

[0021] Step 2.3: Locate the field node at the boundary of the three-dimensional integrated microsystem model, and using the field node as a reference, determine the coordinates of the integration point at the boundary according to the weight of the integration point at the boundary within the background mesh corresponding to the field node boundary.

[0022] Optionally, step 4 includes:

[0023] Step 4.1: Construct the field function within the support domain. The field function is:

[0024]

[0025] Where u(x) is the field function, a(x) is the m-dimensional coefficient vector, and a i (x) is the j-th dimension coefficient vector, p(x) is the basis function, p i (x) is the basis function of the i-th field node, h is the number of the field node, and T is the transpose operation;

[0026] Step 4.2: Within the solution domain, construct a weighted discrete sum of squares J based on the field function. The weighted discrete sum of squares J is:

[0027]

[0028] Where, ω i (x) is the weight function, n is the number of field nodes in the support domain corresponding to the integration point in the weight function, and u i Let be the field function of the i-th field node;

[0029] Step 4.3: Obtain the linear equation based on the minimum value of the weighted discrete sum of squares J. The linear equation is:

[0030] a(x)=A -1 (x)B(x)u s

[0031] Where A(x) is the first matrix, B(x) is the second matrix, and u s It is a vector consisting of the function values ​​of the field function at n field nodes in the support domain;

[0032] u s =[u1 u2 ... u n ] T

[0033] Step 4.4: Obtain the shape function based on the first matrix and the second matrix;

[0034] Step 4.5: Obtain the derivative of the shape function based on the shape function.

[0035] Optionally, step 4.4 includes:

[0036] The transpose of the form function is obtained from the transpose of the basis functions, the first matrix, and the second matrix. The transpose of the form function is:

[0037] Φ T (x)=p T (x)A -1 (x)B(x)

[0038] Where Φ(x) is a shape function.

[0039] Optionally, step 4.4 includes:

[0040] Step 4.41: Perform LU decomposition on the first matrix to obtain the decomposed result;

[0041] Step 4.42: Obtain the auxiliary matrix based on the first matrix, the decomposition result, and the field function;

[0042] Step 4.43: Obtain the shape function based on the transpose of the auxiliary matrix and the second matrix.

[0043] Optionally, step 5 includes:

[0044] Step 5.1: Based on the shape function and its derivative, obtain the foundation temperature stiffness matrix and the foundation nodal thermal load matrix;

[0045] The basic temperature stiffness matrix is:

[0046]

[0047] The basic node thermal load matrix is:

[0048] f i =∫ Ω Φ i dΩ

[0049] Where, k i Let Ω be the basic temperature stiffness matrix at the i-th field node, Ω be the discrete domain of the three-dimensional microsystem, (x,y,z) be the coordinates of the i-th field node, and f be the temperature stiffness matrix at the i-th field node. i Let Φ be the base temperature stiffness matrix at the i-th field node. i Let be the shape function of the i-th field node;

[0050] Step 5.2: Based on the shape function and its derivative, obtain the temperature stiffness matrix and nodal thermal load matrix at the boundary according to the Dirichlet boundary condition, Neumann boundary condition or Robin boundary condition;

[0051] The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Dirichlet boundary conditions are as follows:

[0052]

[0053]

[0054] Where K1 is the first temperature stiffness matrix at the boundary, F1 is the first nodal thermal load matrix at the boundary, α is the penalty function factor, and T f Let Γ1 be the ambient temperature, Γ1 be the boundary influenced by the Dirichlet boundary conditions, and Φ be the boundary temperature. j Let be the shape function of the j-th field node;

[0055] The nodal thermal load matrix at the boundary based on the Neumann boundary condition is as follows:

[0056]

[0057] Where F2 is the second node thermal load matrix at the boundary, q is the heat flux density, and Γ2 is the boundary affected by the Neumann boundary condition;

[0058] The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Robin boundary conditions are as follows:

[0059]

[0060]

[0061] Where K3 is the third temperature stiffness matrix at the boundary, F3 is the third nodal thermal load matrix at the boundary, h is the convective heat transfer coefficient, and Γ3 is the boundary affected by the Robin boundary condition.

[0062] Optionally, the temperature stiffness matrix and the nodal thermal load matrix in the discrete domain are respectively:

[0063] K i =k i ·K(T)

[0064] F i =f i ·P(T)

[0065]

[0066]

[0067] Among them, K i Let F be the temperature stiffness matrix of the i-th field node. i Let c be the nodal thermal load matrix of the i-th field node. k V is the temperature interpolation coefficient for the material, where T is the temperature of the material, T0 to T1 is the interpolation interval, and V is the temperature interpolation coefficient for the material.dd R is the voltage of the microsystem. TSV (T) is the TSV equivalent resistance, ρ(T) is the resistivity, and l TSV For TSV height, r TSV For the TSV radius, f is the frequency of the voltage, μ is the permeability of the metal, and σ(T) is the conductivity.

[0068] Optionally, step 7 includes:

[0069] Step 7.1: Obtain the temperature field T based on the ambient temperature and input power consumption. old Then, the temperature stiffness matrix and nodal thermal load matrix within the corresponding discrete domain are added to the temperature stiffness matrix and nodal thermal load matrix at the boundary to obtain the temperature stiffness matrix K at the current temperature. new and thermal load matrix F new ;

[0070] Step 7.2: Based on the temperature stiffness matrix K at the current temperature. new and thermal load matrix F new Obtain the temperature field T new ;

[0071] Step 7.3: Compare the temperature field T old and temperature field T new If the difference is less than the set precision tol, then output the temperature field T. new If not, then use the temperature field T. new The thermal conductivity and electrical conductivity are recalculated, and a new round of temperature field iteration calculations is performed until the accuracy is less than the set precision tol, so as to obtain the temperature field of the three-dimensional integrated microsystem model.

[0072] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0073] This invention provides a method for electrothermal coupling analysis of three-dimensional integrated microsystems based on the Galerkin method. By discretizing the solution domain into field nodes, a background integration mesh and integration points are established. The weights and Jacobian matrices of the integration points within the mesh are calculated, and the shape function values ​​of the field nodes, the global temperature stiffness matrix, the global nodal load matrix, and the stiffness and load matrices on the boundaries are gradually obtained. Finally, the heat transfer equation is established. Based on the synergistic effect of temperature on the electric field and material properties, temperature-related coupling terms are added to the heat transfer equation to achieve electrothermal coupling. Compared with the finite element method, this invention eliminates the influence of mesh shape on the solution results while ensuring computational accuracy, making it better suited for three-dimensional integrated microsystems with multi-scale characteristics.

[0074] The present invention will now be described in further detail with reference to the accompanying drawings. Attached Figure Description

[0075] Figure 1 This is a flowchart illustrating a three-dimensional integrated microsystem electrothermal co-computation method based on the meshless Galerkin method provided in an embodiment of the present invention.

[0076] Figure 2 This is a flowchart of another electrothermal co-computation method for three-dimensional integrated microsystems based on the meshless Galerkin method provided in this embodiment of the invention;

[0077] Figure 3 This is a schematic diagram of the background integral grid provided in an embodiment of the present invention;

[0078] Figure 4 This is the equivalent circuit diagram of TSV under high-frequency conditions provided in the embodiments of the present invention;

[0079] Figure 5 This is a flowchart of obtaining the temperature field of a three-dimensional integrated microsystem provided in an embodiment of the present invention. Detailed Implementation

[0080] The present invention will be further described in detail below with reference to specific embodiments, but the implementation of the present invention is not limited thereto.

[0081] Example 1

[0082] Please see Figure 1 and Figure 2 , Figure 1 This is a flowchart illustrating a three-dimensional integrated microsystem electrothermal co-computation method based on the meshless Galerkin method provided in an embodiment of the present invention. Figure 2 This is a flowchart of another electrothermal co-calculation method for three-dimensional integrated microsystems based on the meshless Galerkin method provided in this embodiment of the invention. The invention provides a three-dimensional integrated microsystem electrothermal co-calculation method based on the meshless Galerkin method, which includes:

[0083] Step 1: Obtain the field node information of the three-dimensional integrated microsystem model. The field node information includes: the total number of field nodes, the field node number, the field node coordinates, and the field node spacing. The field node spacing is the distance between two adjacent nodes.

[0084] Step 1.1: Obtain the three-dimensional integrated microsystem model.

[0085] It is understandable that the construction process of the three-dimensional integrated microsystem model is the same as that of existing technologies, and will not be elaborated here.

[0086] Step 1.2: Discretize the field nodes of the three-dimensional integrated microsystem model to obtain the information of each field node.

[0087] Step 2: Establish a background integral grid based on the 3D integrated microsystem model, establish integral points using grid and boundary information, and calculate the information of integral points within the background grid. The information of integral points includes the coordinates and weights of the integral points, as well as the total number of integral points, the number of integral points in each background grid, and the grid number of the integral point.

[0088] In this embodiment, the background integral grid and boundary information include: total number of grid points, grid point number, grid point coordinates, total number of grids, connection method of grid points forming the grid, grid number, coordinates of boundary grid points, and grid number of boundary points.

[0089] The integration weights determine the different influences of each integration point in each grid, such as... Figure 3 As shown, the blue circles represent background grid points, and the red * shapes represent integration points. Nodes farther from the grid points have a smaller impact, while nodes closer to the grid points have a larger impact. The weight of the integration points in the background grid is a fixed value.

[0090] In one specific embodiment, step 2 may include:

[0091] Step 2.1: Establish a background integral mesh based on the three-dimensional integrated microsystem model. The background integral mesh covers the three-dimensional integrated microsystem model and includes several mesh points arranged in a matrix.

[0092] Specifically, the background integration grid is a six-sided cube (i.e., a cubic) in the form of a background integration grid. The background integration grid consists of multiple background grids, each of which is composed of 6 grid points connected together. The spacing between two adjacent grid points is 2 to 3 times the spacing between field nodes.

[0093] Step 2.2: Within each background grid, based on the number of integration points in one dimension, the weight of the integration point is determined using a Gaussian function, and the coordinates of the integration point are obtained by multiplying the distance between two adjacent grid points in a background grid with the weight of the integration point.

[0094] Specifically, first, the number of integration points in one dimension is set. Then, the weight value corresponding to each integration point is found in the Gaussian curve for the corresponding number of integration points. This determines the weight of each integration point. After that, for a background grid, the product result is obtained by multiplying the distance between two adjacent grid points and the weight of the integration point, and the product result is used as the coordinate of the integration point.

[0095] Optionally, each background grid can be divided into 3 to 9 integration points in one dimension.

[0096] Step 2.3: Locate the field node at the boundary of the 3D integrated microsystem model, and use the field node as a reference. Within the background mesh corresponding to the field node boundary, determine the coordinates of the integration point at the boundary according to the weight of the integration point at the boundary. That is, for the background mesh, use the product of the distance between two adjacent mesh points and the weight of the integration point to obtain the coordinates of the integration point at the boundary.

[0097] Step 3: Locate the integration points within the solution domain and determine the field nodes within the support domain of each integration point.

[0098] In this embodiment, the solution domain is the region obtained by drawing a circle with the field node as the center and the average spacing between the field nodes in the background grid as the radius. The support domain is the region obtained by drawing a circle with the integration point in the background grid as the center and the average spacing between the integration points as the radius. The average spacing is the average of the spacing between two adjacent points.

[0099] The information supported for field nodes within the domain includes: the field node's ID and the field node's coordinates.

[0100] In one specific embodiment, step 3 may include:

[0101] Step 3.1: Determine the size and shape of the solution domain for each field node.

[0102] Step 3.2: Traverse the integration points. If the integration point is within the solution domain, find field nodes whose distance from the integration point is less than the radius of the support domain.

[0103] Step 4: Calculate the shape function and its derivative at each integration point.

[0104] Step 4.1: Construct the field function within the support domain.

[0105] The shape function is constructed using the moving least squares method. Within the solution domain Ω, it is assumed that the field values ​​u of the field function u(x) at the n nodes are known, i.e.:

[0106] u i =u(x i )

[0107] Where u(x) i Let be the field function of the i-th field node.

[0108] Within the support domain of the calculated field nodes, the approximate form of the field function is constructed as follows:

[0109]

[0110] Where u(x) is the field function, a(x) is the m-dimensional coefficient vector, and a j(x) is the j-th dimension coefficient vector, h is the field node number, T is the transpose operation, p(x) is the basis function, p i (x) is the basis function of the i-th field node. To ensure its minimum completeness, the basis function is usually a monomial determined by the Pascal triangle.

[0111] Step 4.2: Within the solution domain, construct the weighted discrete sum of squares J based on the field function. The weighted discrete sum of squares J is:

[0112]

[0113] Where, ω i (x) is the weight function, ω i (x)=ω(xx i ), with field node x i As the distance between integration points increases or decreases, the value of the weight function gradually decreases, where n is the number of field nodes in the support domain corresponding to the integration point in the weight function, and u i Let be the field function of the i-th field node.

[0114] Step 4.3: Obtain the linear equation based on the minimum value of the weighted discrete sum of squares J. The linear equation is:

[0115] a(x)=A -1 (x)B(x)u s

[0116] Where A(x) is the first matrix, B(x) is the second matrix, and u s It is a vector consisting of the function values ​​of the field function at n field nodes in the support domain.

[0117] The first matrix A(x) is:

[0118]

[0119] The second matrix B(x) is:

[0120]

[0121] The first matrix A(x) is an m×m square matrix, and the second matrix B(x) is an m×n matrix. Substituting a(x) into the approximate expression of the field function, we obtain:

[0122]

[0123] Where, Φ i (x) is the shape function at the i-th field node, and is a function of the coordinates (x,y,z).

[0124] Step 4.4: Obtain the shape function based on the first and second matrices.

[0125] This invention provides two methods for shape functions, which will be described below.

[0126] The first method is:

[0127] The transpose of the form function is obtained from the transpose of the basis functions, the first matrix, and the second matrix. The transpose of the form function is:

[0128] Φ T (x)=p T (x)A -1 (x)B(x)

[0129] Where Φ(x) is a shape function.

[0130] The second method is:

[0131] Step 4.41: Perform LU decomposition on the first matrix to obtain the decomposed result, i.e., A = LU.

[0132] Step 4.42: Obtain the auxiliary matrix based on the first matrix, the decomposition result, and the field function.

[0133] A(x)c(x)=p(x)

[0134] c(x)=U -1 L -1 p(x)

[0135] Where c(x) is an auxiliary matrix.

[0136] Step 4.43: Obtain the shape function based on the transpose of the auxiliary matrix and the second matrix. The shape function is:

[0137] Φ(x)=c T (x)B(x).

[0138] Step 5: Based on the shape function and its derivative, calculate the temperature stiffness matrix and nodal thermal load matrix at the boundary, and calculate the contribution of each integration point to the temperature stiffness matrix and nodal load matrix to obtain the basic temperature stiffness matrix and basic nodal thermal load matrix.

[0139] The fundamental governing equation for steady-state heat transfer is:

[0140]

[0141] Step 5.1: Based on the shape function and its derivative, obtain the basic temperature stiffness matrix and the basic nodal thermal load matrix.

[0142] Specifically, by selecting temperature variation as the test function and using the weighted residual method, the basic temperature stiffness matrix and basic nodal thermal load matrix in the Galerkin discrete control equations for steady-state heat transfer can be obtained.

[0143] The basic temperature stiffness matrix is:

[0144]

[0145] The basic node thermal load matrix is:

[0146] f i =∫ Ω Φ i dΩ

[0147] Where, k i Let Ω be the basic temperature stiffness matrix at the i-th field node, Ω be the discrete domain of the three-dimensional microsystem, (x,y,z) be the coordinates of the i-th field node, and f be the temperature stiffness matrix at the i-th field node. i Let Φ be the base temperature stiffness matrix at the i-th field node. i Let be the shape function of the i-th field node;

[0148] Step 5.2: Based on the shape function and its derivative, obtain the temperature stiffness matrix and nodal thermal load matrix at the boundary according to the Dirichlet boundary condition, Neumann boundary condition or Robin boundary condition.

[0149] The Dirichlet heat transfer boundary is:

[0150] The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Dirichlet boundary conditions are as follows:

[0151]

[0152]

[0153] Where K1 is the first temperature stiffness matrix at the boundary, F1 is the first nodal thermal load matrix at the boundary, α is the penalty function factor, and T f Let Γ1 be the ambient temperature, Γ1 be the boundary influenced by the Dirichlet boundary conditions, and Φ be the boundary temperature. j Let be the shape function of the j-th field node.

[0154] The Neumann heat transfer boundary is:

[0155] The nodal thermal load matrix at the boundary based on the Neumann boundary condition is then:

[0156]

[0157] Where F2 is the second node thermal load matrix at the boundary, q is the heat flux density, and Γ2 is the boundary affected by the Neumann boundary condition;

[0158] The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Robin boundary conditions are as follows:

[0159]

[0160]

[0161] Where K3 is the third temperature stiffness matrix at the boundary, F3 is the third nodal thermal load matrix at the boundary, h is the convective heat transfer coefficient, and Γ3 is the boundary affected by the Robin boundary condition.

[0162] Step 6: Based on the influence of temperature on the three-dimensional integrated microsystem model, add electrothermal coupling terms to the basic temperature stiffness matrix and the basic nodal thermal load matrix to obtain the temperature stiffness matrix and nodal thermal load matrix in the discrete domain, where the discrete domain is the discrete three-dimensional integrated microsystem model.

[0163] Specifically, the influence of temperature on the electrothermal coupling process of a microsystem is mainly reflected in two physical quantities: the thermal conductivity and electrical conductivity of the material. Both physical quantities are calculated based on temperature interpolation coefficients, namely:

[0164]

[0165] Where K(T) is the thermal conductivity, and c k is the temperature difference coefficient, T is the temperature of the material, and T0 to T1 is the interpolation interval.

[0166] according to Figure 4 The equivalent circuit of the three-dimensional integrated interconnect structure shown has a TSV containing only resistors, capacitors, and inductors, making it a purely resistive circuit. The power dissipation through the TSV in high-frequency circuits is:

[0167]

[0168] Where P(T) is the power consumption, V dd R is the voltage of TSV. TSV (T) is the equivalent resistance of TSV, ρ(T) is the resistivity, and l TSV r is the height of TSV TSV Let TSV be the radius. f is the power supply frequency, μ is the permeability of the metal, and σ(T) is the conductivity.

[0169] Therefore, the temperature stiffness matrix and the nodal thermal load matrix in the discrete domain are:

[0170] K i =k i ·K(T)

[0171] Fi =f i ·P(T)

[0172] Among them, K i Let F be the temperature stiffness matrix of the i-th field node. i Let be the nodal thermal load matrix of the i-th field node.

[0173] Step 7: Based on the heat conduction control equation, obtain the temperature field of the three-dimensional integrated microsystem model according to the temperature stiffness matrix and nodal thermal load matrix at the boundary, as well as the temperature stiffness matrix and nodal thermal load matrix in the discrete domain.

[0174] In one specific embodiment, please refer to Figure 5 Step 7 may include:

[0175] Step 7.1: Obtain the temperature field T based on the ambient temperature and input power consumption. old Then, the temperature stiffness matrix and nodal thermal load matrix in the corresponding discrete domain are added to the temperature stiffness matrix and nodal thermal load matrix at the boundary to obtain the temperature stiffness matrix K at the current temperature. new and thermal load matrix F new .

[0176] Specifically, the ambient temperature and power consumption are first input to obtain the current temperature field T. old According to the temperature field T old Calculate the current thermal conductivity K(T) old ) and conductivity σ(T) old Then, the basic temperature stiffness matrix k i and node load matrix f i respectively with thermal conductivity K(T) old ), conductivity σ(T) old Multiplying the power consumption P(T) with the corresponding power consumption P(T) yields the temperature stiffness matrix K in the discrete domain. i and nodal thermal load matrix F i Then, the temperature stiffness matrix and nodal thermal load matrix in the discrete domain are added to the temperature stiffness matrix and nodal thermal load matrix at the boundary to obtain the temperature stiffness matrix K at the current temperature. new and thermal load matrix F new .

[0177] Step 7.2: Based on the temperature stiffness matrix K at the current temperature. new and thermal load matrix F new Obtain the temperature field T new K new T new =F new .

[0178] Step 7.3: Compare the temperature field Told and temperature field T new If the difference is less than the set precision tol, then output the temperature field T. new If not, then use the temperature field T. new The thermal conductivity and electrical conductivity are recalculated, and a new round of temperature field iteration calculations is performed until the accuracy is less than the set precision tol, in order to obtain the temperature field of the three-dimensional integrated microsystem model.

[0179] To address the challenges of large structural size variations and time-consuming finite element analysis in 3D microsystems, the Galerkin method discretizes the solution domain into field nodes, establishing a background integration mesh and integration points, thus avoiding the mesh generation process and improving computational efficiency. Furthermore, it calculates the weights and Jacobian matrices of the integration points within the mesh, progressively deriving the shape function values ​​of the field nodes, the global temperature stiffness matrix, the global nodal load matrix, and the stiffness and load matrices at the boundaries, finally establishing the heat transfer equation. Based on the synergistic effect of temperature on the electric field and material properties, temperature-related coupling terms are added to the heat transfer equation to achieve electrothermal coupling. Compared to the finite element method, this invention eliminates the influence of mesh shape on the solution results while ensuring computational accuracy, making it better suited for 3D integrated microsystems with multi-scale characteristics.

[0180] Compared to the finite element method, this invention enables rapid electrothermal synergistic calculation of three-dimensional integrated microsystems, satisfying the electrothermal synergistic analysis of structures with arbitrary geometric scales in three-dimensional integrated microsystems. Simultaneously, this invention achieves similar accuracy to the finite element method, meeting the accuracy requirements of electrothermal coupling analysis. Furthermore, compared to the finite element method, the greater the difference in the geometric dimensions of the microsystem structure, the shorter the computation time required by this invention, significantly saving time and utilizing computational resources, effectively realizing the electrothermal synergistic distribution calculation of multi-scale structures in three-dimensional integrated microsystems.

[0181] It should be noted that the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, features defined as "first" or "second" may explicitly or implicitly include one or more features. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0182] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features or characteristics described may be combined in any suitable manner in one or more embodiments or examples. In addition, those skilled in the art can combine and integrate the different embodiments or examples described in this specification.

[0183] Although the invention has been described herein in conjunction with various embodiments, those skilled in the art will understand and implement other variations of the disclosed embodiments by reviewing the accompanying drawings and the disclosure in carrying out the claimed invention. In this specification, the word "comprising" does not exclude other components or steps, and "a" or "an" does not exclude a plurality. While certain measures are described in different embodiments, this does not mean that these measures cannot be combined to produce good results.

[0184] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for electrothermal co-calculation of a three-dimensional integrated microsystem based on the meshless Galerkin method, characterized in that, The three-dimensional integrated microsystem electrothermal co-computation method includes: Step 1: Obtain the field node information of the three-dimensional integrated microsystem model; Step 2: Establish a background integral mesh based on the three-dimensional integrated microsystem model, establish integration points using mesh and boundary information, and calculate the information of the integration points within the background mesh. The information of the integration points includes the coordinates and weights of the integration points. Step 3: Locate the integration points within the solution domain and determine the field nodes within the support domain of each integration point; Step 4: Calculate the shape function and its derivative at each integration point; Step 5: Based on the shape function and its derivative, calculate the temperature stiffness matrix and nodal thermal load matrix at the boundary, and calculate the contribution of each integration point to the temperature stiffness matrix and nodal load matrix to obtain the basic temperature stiffness matrix and basic nodal thermal load matrix. Step 6: Based on the influence of temperature on the three-dimensional integrated microsystem model, add electrothermal coupling terms to the basic temperature stiffness matrix and the basic nodal thermal load matrix to obtain the temperature stiffness matrix and nodal thermal load matrix in the discrete domain; the temperature stiffness matrix and nodal thermal load matrix in the discrete domain are as follows: in, For the first Temperature stiffness matrix of each field node Based on the basic temperature stiffness matrix, For the first The nodal thermal load matrix of each field node Based on the thermal load matrix of the basic nodes, This is the temperature difference coefficient. For the temperature of the material, arrive For the interpolation interval, The voltage of TSV. This is the equivalent resistance of TSV. Resistivity For the height of TSV, Let TSV be the radius. , For power supply frequency, The permeability of the metal, Electrical conductivity; Step 7: Based on the heat conduction control equation, obtain the temperature field of the three-dimensional integrated microsystem model according to the temperature stiffness matrix and nodal thermal load matrix at the boundary, as well as the temperature stiffness matrix and nodal thermal load matrix in the discrete domain.

2. The method for electrothermal co-computation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 1, characterized in that, Step 1 includes: Step 1.1: Obtain the three-dimensional integrated microsystem model; Step 1.2: Discretize the field nodes of the three-dimensional integrated microsystem model to obtain the field node information, which includes: the total number of field nodes, the field node number, the field node coordinates, and the field node spacing.

3. The method for electrothermal co-computation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 1, characterized in that, Step 2 includes: Step 2.1: Establish the background integral mesh based on the three-dimensional integrated microsystem model, wherein the background integral mesh covers the three-dimensional integrated microsystem model, and the background integral mesh includes several mesh points arranged in a matrix form; Step 2.2: Within each background grid, based on the number of integration points in one dimension, the weight of the integration point is determined using a Gaussian function, and the coordinates of the integration point are obtained by multiplying the distance between two adjacent grid points in a background grid with the weight of the integration point. Step 2.3: Locate the field node at the boundary of the three-dimensional integrated microsystem model, and using the field node as a reference, determine the coordinates of the integration point at the boundary according to the weight of the integration point at the boundary within the background mesh corresponding to the field node boundary.

4. The method for electrothermal co-calculation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 1, characterized in that, Step 4 includes: Step 4.1: Construct the field function within the support domain. The field function is: in, For the field function, for dimensional coefficient vector, For the first dimensional coefficient vector, As basis functions, For the first basis functions of each field node This refers to the field node number. This is a transpose operation; Step 4.2: Within the solution domain, construct a weighted discrete sum of squares based on the field function. The weighted discrete sum of squares for: in, For the weight function, This represents the number of field nodes in the support domain corresponding to the integration points in the weight function. For the first The field function of each field node; Step 4.3: Based on the weighted discrete sum of squares... The minimum value of is used to obtain a linear equation, which is: in, For the first matrix, For the second matrix, It is a vector consisting of the function values ​​of the field function at n field nodes in the support domain; Step 4.4: Obtain the shape function based on the first matrix and the second matrix; Step 4.5: Obtain the derivative of the shape function based on the shape function.

5. The method for electrothermal co-calculation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 4, characterized in that, Step 4.4 includes: The transpose of the form function is obtained from the transpose of the basis functions, the first matrix, and the second matrix. The transpose of the form function is: in, It is a shape function.

6. The method for electrothermal co-calculation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 4, characterized in that, Step 4.4 includes: Step 4.41: Perform LU decomposition on the first matrix to obtain the decomposed result; Step 4.42: Obtain the auxiliary matrix based on the first matrix, the decomposition result, and the field function; Step 4.43: Obtain the shape function based on the transpose of the auxiliary matrix and the second matrix.

7. The method for electrothermal co-calculation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 1, characterized in that, Step 5 includes: Step 5.1: Based on the shape function and its derivative, obtain the foundation temperature stiffness matrix and the foundation nodal thermal load matrix; The basic temperature stiffness matrix is: The basic node thermal load matrix is: in, For the first The basic temperature stiffness matrix at each field node. For the discrete domain of a three-dimensional microsystem, For the first The coordinates of each field node, For the first The basic temperature stiffness matrix at each field node. For the first Shape functions of each field node; Step 5.2: Based on the shape function and its derivative, obtain the temperature stiffness matrix and nodal thermal load matrix at the boundary according to the Dirichlet boundary condition, Neumann boundary condition or Robin boundary condition; The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Dirichlet boundary conditions are as follows: in, Here is the first temperature stiffness matrix at the boundary. The thermal load matrix of the first node at the boundary. As the penalty function factor, For ambient temperature, The boundary affected by the Dirichlet boundary conditions. For the first Shape functions of each field node; The nodal thermal load matrix at the boundary based on the Neumann boundary condition is as follows: in, The thermal load matrix of the second node at the boundary. For heat flux density, The boundary affected by Neumann boundary conditions; The temperature stiffness matrix and nodal thermal load matrix at the boundary based on the Robin boundary conditions are as follows: in, This is the third temperature stiffness matrix at the boundary. The thermal load matrix of the third node at the boundary. The convective heat transfer coefficient, The boundary affected by Robin's boundary conditions.

8. The method for electrothermal co-computation of a three-dimensional integrated microsystem based on the meshless Galerkin method according to claim 1, characterized in that, Step 7 includes: Step 7.1: Obtain the temperature field based on the ambient temperature and input power consumption. Then, the temperature stiffness matrix and nodal thermal load matrix within the corresponding discrete domain are added to the temperature stiffness matrix and nodal thermal load matrix at the boundary to obtain the temperature stiffness matrix at the current temperature. and thermal load matrix ; Step 7.2: Based on the temperature stiffness matrix at the current temperature. and thermal load matrix Obtain the temperature field ; Step 7.3: Compare the temperature fields and temperature field Is the difference less than the set precision? If so, then output the temperature field. If not, then use the temperature field. Recalculate the thermal conductivity and electrical conductivity, and perform a new round of iterative temperature field calculations until the accuracy is less than the set target. In order to obtain the temperature field of the three-dimensional integrated microsystem model.

Images