A Fast Prediction Method for Thermal Field of Electronic Devices with Strong Multiphysics Coupling

Abstract

This invention belongs to the field of thermal management technology for electronic devices. It discloses a rapid prediction method for the thermal field of electronic devices under strong multi-physics coupling. The invention constructs a hierarchical heat source intensity distribution model, develops a grid mapping relationship topology graph technique and a boundary node energy conservation compensation mechanism, designs a multi-timescale coupling transfer function, and introduces a dynamic adjustment strategy based on phase deviation analysis. This invention accurately captures the non-uniform heat source distribution and dynamic characteristics of heat diffusion under high-frequency conditions, ensuring energy consistency in multi-physics coupling and improving the accuracy of thermal field prediction. It has particularly outstanding advantages in predicting the phase and amplitude of transient thermal responses, while significantly reducing computational resource consumption.

Description

Technical Field

[0001] This invention relates to the field of thermal management technology for electronic devices, and more specifically, to a method for rapid prediction of thermal fields in electronic devices with strong multi-physics coupling. Background Technology

[0002] As chip integration increases and operating frequencies reach the GHz range, the limitations of traditional thermal analysis methods become increasingly apparent. At high frequencies, the skin effect causes current to concentrate primarily in an extremely thin layer on the conductor's surface, resulting in strong spatial non-uniformity in heat source distribution. Mainstream thermal models still employ the simplified assumption of uniform heating across the conductor cross-section, ignoring this microscopic distribution difference. While this simplification is acceptable in static analysis, it leads to severe phase lag and amplitude attenuation errors in nanosecond-level transient thermal response prediction, making it impossible to accurately capture the thermal transient characteristics of high-speed electronic systems. Even more challenging is the inherent difference in multiphysics grid systems. Electromagnetic simulations, for accurate calculation of the skin effect, typically use extremely high-density grids in surface regions; while thermal simulations focus more on volumetric heat conduction, resulting in relatively uniform grid distribution. This fundamental difference in grid structure makes the mapping process from electromagnetic loss to heat source exceptionally complex. Traditional interpolation algorithms often exhibit severe violations of energy conservation when dealing with boundary regions with large grid density differences, leading to a significant deviation between the total heat source power and the actual electromagnetic loss. Furthermore, existing methods generally neglect the temporal differences in heat diffusion characteristics of heat sources at different depths and fail to establish dynamic coupling models that consider multi-timescale effects, resulting in systematic deviations in the time series of temperature responses. These problems are particularly prominent in the thermal design of high-power radio frequency equipment, high-speed communication systems, and advanced computing platforms, leading to the adoption of excessive safety margins in engineering practice, increasing the burden on cooling systems, and limiting further improvements in equipment performance.

[0003] In view of this, the present invention proposes a method for rapid prediction of the thermal field of electronic devices with strong multi-physics coupling to solve the above problems. Summary of the Invention

[0004] To overcome the aforementioned deficiencies of the prior art and to achieve the above objectives, the present invention provides the following technical solution: a method for rapid prediction of the thermal field of electronic devices with strong multi-physics coupling, comprising:

[0005] Acquire electromagnetic field simulation data and initial thermal simulation mesh data of electronic devices at their operating frequency. The electromagnetic field simulation data includes the spatial distribution of current density inside the conductor.

[0006] Gradient analysis is performed on the spatial distribution of current density along the radial direction of the conductor to identify the radial depth at which the current density decays to a preset proportion of the surface peak, which is denoted as the skin feature depth.

[0007] Based on the skin feature depth, the conductor cross section is divided into multiple annular sub-regions, the current density integral value of each annular sub-region is calculated, and a layered heat source intensity distribution model is constructed.

[0008] Extract the spatial locations of grid nodes from the initial thermal simulation grid data and electromagnetic field simulation data, calculate the nearest distance and local Jacobian determinant between grid node pairs, and construct a topological graph of grid mapping relationship;

[0009] Based on the mesh mapping relationship topology, the hierarchical heat source intensity distribution model is mapped to the initial thermal simulation mesh data, and the total power deviation before and after mapping is calculated.

[0010] Based on the total power deviation, an energy conservation compensation source term is introduced at the boundary nodes of the grid mapping topology to generate a corrected heat source distribution field;

[0011] A multi-timescale coupled transfer function from electromagnetic loss to temperature rise is constructed. The multi-timescale coupled transfer function includes the time constants of the thermal diffusion characteristics of each level of the hierarchical heat source intensity distribution model.

[0012] A multi-timescale coupled transfer function is applied to the corrected heat source distribution field to generate transient thermal field prediction results.

[0013] The characteristic frequency components of the temperature response in the transient thermal field prediction results are extracted and phase deviation analysis is performed between them and the excitation frequency of the electromagnetic field simulation data.

[0014] Based on the phase deviation analysis results, the heat source delay phase of each annular sub-region in the layered heat source intensity distribution model is dynamically adjusted to update the transient thermal field prediction results, thereby achieving rapid and accurate prediction of the thermal field of electronic devices.

[0015] The technical effects and advantages of this invention for a rapid prediction method of thermal field in electronic devices with strong multi-physics coupling are as follows:

[0016] This invention accurately captures the non-uniform spatial distribution characteristics of heat sources under high-frequency conditions, making thermal field analysis more closely aligned with physical reality. It achieves high-precision mapping between grid systems, effectively eliminating energy leakage or accumulation in boundary regions and ensuring energy conservation from the electromagnetic field to the thermal field, thus laying a reliable physical foundation for thermal field calculations. By considering the differences in thermal diffusion time characteristics of heat sources at different depths, it successfully simulates the multi-timescale effects of the thermal field, improving the temporal accuracy of transient temperature prediction. Furthermore, by optimizing the temporal phase characteristics of the heat source model, this invention effectively eliminates phase and amplitude deviations in the temperature response, providing a precise basis for the thermal reliability analysis of high-frequency electronic products. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the rapid prediction method for the thermal field of electronic devices with strong multi-physics coupling according to the present invention. Detailed Implementation

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

[0019] Please see Figure 1 In this embodiment of the invention, the specific implementation process of the method for rapid prediction of thermal field of electronic devices with strong multi-physics coupling includes:

[0020] Electromagnetic field simulation data and initial thermal simulation mesh data of the electronic device at its operating frequency were acquired. The electromagnetic field simulation data includes the spatial distribution of current density inside the conductor. The electromagnetic field simulation data records key electromagnetic parameters such as the electric field distribution, magnetic field distribution, and current density distribution of the conductor in the electronic device at a specific operating frequency. The initial thermal simulation mesh data contains the geometric mesh generation and material thermal property definitions required for thermal field solving. The spatial distribution of current density reflects the electromagnetic loss distribution characteristics inside the conductor and is the main heat source input for thermal field prediction. These data provide fundamental physical field information for subsequent analysis, ensuring the integrity and accuracy of multiphysics coupling analysis.

[0021] Gradient analysis is performed on the spatial distribution of current density along the radial direction of the conductor to identify the radial depth at which the current density decays to a predetermined proportion of the surface peak value, denoted as the skin depth. The skin depth reflects the distribution characteristics of high-frequency current within the conductor and is crucial for accurately modeling heat source distribution. Under high-frequency conditions, the current density decays exponentially from the conductor surface inwards; gradient analysis can accurately capture this decay characteristic. This step quantitatively analyzes the spatial distribution gradient of current density to identify the depth location where the current density decreases to a specific proportion (typically 1 / e or 36.8%) of the surface peak value, serving as a characteristic scale of the skin effect and providing key parameters for subsequent layered heat source modeling.

[0022] Based on skin depth, the conductor cross-section is divided into multiple annular sub-regions. The integral value of the current density in each annular sub-region is calculated to construct a layered heat source intensity distribution model. This layered heat source intensity distribution model is a key technique for accurately characterizing the non-uniform heat source distribution caused by high-frequency current, capturing the radial decay characteristics of current density through spatial layering. The modeling process divides the conductor cross-section from the surface to the center into multiple concentric annular sub-regions. The division interval is dynamically adjusted according to the current density decay characteristics, ensuring more detailed division in the surface region with large gradient changes. Volume integration is performed on the current density within each annular sub-region to obtain the equivalent current value, and the heat source intensity density of each layer is calculated in conjunction with the material resistivity. This layered model intuitively reflects the spatial non-uniform distribution of heat sources under high-frequency conditions, providing more accurate heat source input for subsequent thermal field calculations.

[0023] The spatial positions of grid nodes in the initial thermal simulation grid data and electromagnetic field simulation data are extracted. The nearest distance and local Jacobian determinant between grid node pairs are calculated to construct a grid mapping topology. This topology serves as a bridge for converting electromagnetic and thermal field data, solving the data mapping problem between different grid systems in multiphysics simulations. The construction process first compares the spatial positions of grid nodes in the electromagnetic and thermal fields to identify correspondences; then, it calculates the spatial distance and local transformation characteristics between node pairs to quantify mapping accuracy; finally, it constructs a complete mapping topology, recording the spatial correspondences and mapping weights between grids. This technique ensures that electromagnetic loss data can be accurately and efficiently converted into the heat source distribution required for thermal field solutions, guaranteeing the accuracy of multiphysics coupling calculations.

[0024] Based on the topology of the mesh mapping relationship, the layered heat source intensity distribution model is mapped to the initial thermal simulation mesh data, and the total power deviation before and after mapping is calculated. Heat source distribution mapping is a core step in multiphysics coupling, ensuring that electromagnetic losses are correctly converted into heat source input for the thermal field. The mapping process employs spatial interpolation techniques based on the topology graph to accurately convert the data of the layered heat source model onto the thermal simulation mesh while maintaining total energy conservation. By calculating the comparison of total power before and after mapping, the accuracy of energy conversion is quantified, providing a basis for subsequent compensation and correction. This step solves the data conversion problem between different mesh systems of the electromagnetic field and the thermal field, ensuring the spatial accuracy and energy consistency of the heat source distribution.

[0025] Based on the total power deviation, energy conservation compensation source terms are introduced at the boundary nodes of the mesh mapping topology to generate a corrected heat source distribution field. Energy conservation compensation is a key technique to ensure the accuracy of multiphysics coupling, ensuring energy conservation constraints through local heat source intensity adjustment. The compensation process first identifies the thermal mesh cells in the mesh boundary region as the main carriers of compensation; then, based on the total power deviation and cell volume ratio, the compensation amount for each boundary cell is calculated; finally, the compensated heat source is added to the original heat source distribution to generate a corrected heat source field that satisfies energy conservation. This boundary compensation method not only ensures total energy conservation but also maintains the spatial characteristics of the heat source distribution, improving the physical accuracy of multiphysics coupling calculations.

[0026] A multi-timescale coupled transfer function (MTF) is constructed to describe the dynamic relationship between electromagnetic loss and temperature rise. This MTF incorporates the characteristic time constants of heat diffusion at each level of the hierarchical heat source intensity distribution model. The MTF is a mathematical model that describes the dynamic relationship between electromagnetic loss and temperature response, considering the multi-scale characteristics of the heat diffusion process. The construction process first calculates the characteristic time constants of heat diffusion in each annular sub-region to characterize the dynamic characteristics of heat transfer. Then, an equivalent circuit of a multi-order RC thermal network is established to simulate the diffusion process of heat sources at different depths. Finally, the time-domain form of the coupled transfer function is obtained through Laplace transform and inverse transform. This function accurately describes the dynamic process of temperature rise caused by electromagnetic loss, considering the multi-timescale effects of heat diffusion, and provides a dynamic basis for transient thermal field prediction.

[0027] A multi-timescale coupled transfer function is applied to the corrected heat source distribution field to generate transient thermal field prediction results. Transient thermal field prediction is the core step in realizing dynamic thermal field analysis, calculating the dynamic response of heat sources to temperature through transfer functions. The prediction process first groups the heat source distribution into annular sub-regions to form a time series input; then, the corresponding transfer function is applied to the heat sources in each region to calculate the local temperature rise response; finally, the temperature rise responses of each region are superimposed and combined with environmental boundary conditions to generate a complete transient temperature field. This transfer function-based method greatly improves the efficiency of thermal field calculation while maintaining physical accuracy, enabling rapid prediction of temperature distribution and hotspot evolution processes of electronic devices under complex operating conditions.

[0028] Characteristic frequency components of the temperature response are extracted from transient thermal field prediction results, and phase deviation analysis is performed between these components and the excitation frequency of the electromagnetic field simulation data. Phase deviation analysis is an important tool for evaluating the accuracy of thermal field predictions, revealing the phase relationship between the temperature response and electromagnetic excitation through frequency domain analysis. The analysis process first extracts the spectral characteristics of the temperature time series to identify the dominant frequency components; then, it compares and matches these components with the electromagnetic excitation frequency to calculate the phase difference; finally, it analyzes the spatial distribution of the phase deviation to identify anomalous regions. This frequency domain analysis method can accurately assess the dynamic characteristics of the thermal response, discover potential inaccuracies in the transfer function, and provide an important basis for model optimization.

[0029] Based on phase deviation analysis results, the heat source delay phase of each annular sub-region in the layered heat source intensity distribution model is dynamically adjusted to update the transient thermal field prediction results, achieving rapid and accurate prediction of the thermal field of electronic devices. Dynamic phase adjustment is a closed-loop optimization technique to improve the accuracy of thermal field prediction by adjusting the heat source model parameters. The adjustment process first calculates the time delay of each annular sub-region based on the phase deviation; then, it applies the corresponding time-shift operation to the heat source time series; finally, it recalculates the temperature rise response to generate updated thermal field prediction results. This feedback-based dynamic adjustment mechanism greatly improves the accuracy of thermal field prediction, especially under complex multi-frequency excitation and nonlinear material properties, accurately capturing the dynamic evolution characteristics of the thermal field and providing a reliable basis for the thermal design and management of electronic devices.

[0030] In this embodiment of the invention, gradient analysis is performed on the spatial distribution of current density along the radial direction of the conductor to identify the radial depth at which the current density decays to a preset proportion of the surface peak, denoted as the skin feature depth, including:

[0031] Starting from nodes on the conductor surface, a current density sampling sequence is extracted along the radial direction of the conductor. Radial sampling is a fundamental step in capturing the skin effect of high-frequency current, revealing the current distribution pattern through systematic sampling along a specific direction. The sampling process first identifies representative nodes on the conductor surface as sampling starting points; then, a sampling path is set along the radial direction, typically choosing a radial path perpendicular to the electric field direction to maximize the capture of skin effect characteristics; finally, dense sampling points are set along the path to extract the current density values ​​at each point, forming a complete spatial distribution sequence. The sampling interval is dynamically adjusted according to the conductor size and the expected skin depth to ensure denser sampling in areas with large gradient changes; typically, the sampling density in the surface region is 3-5 times that in the inner region to capture rapidly changing current distribution characteristics.

[0032] An exponential fitting is performed on the current density sampling sequence to obtain the decay constant of the current density. Exponential fitting is a key technique for quantifying the skin effect characteristics, accurately describing the spatial decay law of current density through a mathematical model. The fitting process is based on the physical principle of the skin effect and employs an exponential function model. The sampled data were fitted using least squares, where Distance from surface Current density at that point Surface current density, The attenuation constant to be fitted is denoted as . The fitting algorithm employs iterative weighted least squares, assigning higher weights to high-gradient regions near the surface to improve fitting accuracy. The attenuation constant obtained through fitting directly reflects the penetration depth characteristics of current within the conductor and is a fundamental parameter for skin depth calculation.

[0033] Based on the attenuation constant, the radial position corresponding to the current density decaying to 1 / e times the surface peak value is calculated and denoted as the initial skin depth. Calculating the initial skin depth is a direct method for quantitatively determining the characteristic scale of the skin effect, based on the standard definition in classical skin effect theory. The calculation process directly applies an exponential model to determine the radial position where the current density drops to 36.8% (i.e., 1 / e times) of the surface value; this position corresponds to the theoretical skin depth. This theoretical depth provides a preliminary assessment of the skin effect, but it needs to be verified and corrected through subsequent steps combined with actual current distribution to ensure the practicality of the depth definition.

[0034] Centered on the initial skin depth, a predetermined radial distance is extended inwards and outwards to form a skin effect evaluation window. Setting the evaluation window is a preparatory step for a comprehensive skin effect evaluation, ensuring its comprehensiveness by defining an appropriate spatial range. The window is set with the initial skin depth as the center, extending a certain proportion in both the surface and inward directions, typically 50% outwards and 100% inwards, forming a complete evaluation range. This asymmetric extension considers the more drastic changes in current density in the surface region, ensuring the window can fully capture the current distribution in key areas. The size of the evaluation window is dynamically adjusted, adaptively determined based on the operating frequency and material properties, ensuring the window range can fully capture the skin effect characteristics without introducing interference from irrelevant areas.

[0035] The skin concentration is calculated as the ratio of the integral of the current density within the skin effect assessment window to the total current in the conductor. The skin concentration is a quantitative indicator for assessing the significance of the skin effect, revealing the degree of current concentration through energy distribution characteristics. The calculation process first involves volume integration of the current density within the assessment window to obtain the total current within the window; then, the total current across the entire conductor cross-section is calculated; finally, the ratio of the two is calculated to obtain the skin concentration. The integration ranges are the assessment window and the entire conductor volume, respectively. A higher concentration indicates that the current is more concentrated in the surface region, and the skin effect is more significant; conversely, a lower concentration indicates a relatively uniform current distribution and a less pronounced skin effect.

[0036] When the skin concentration exceeds a preset concentration threshold, the outer boundary depth of the skin effect evaluation window is marked as the skin feature depth. Determining the feature depth is the final step in skin effect analysis, defining a practical skin feature scale using a concentration criterion. The determination process is based on a comparison between the skin concentration and a preset threshold. When the concentration exceeds the threshold (typically set to 0.85-0.95, indicating that the window contains 85%-95% of the total current), the outer boundary of the evaluation window (i.e., the boundary away from the conductor surface) is defined as the skin feature depth. This energy distribution-based definition method better meets the needs of practical engineering applications and can more accurately guide subsequent layered modeling compared to theoretical depth. If the concentration does not reach the threshold, the evaluation window size needs to be adjusted, and the calculation process repeated until the concentration requirement is met. The finally determined skin feature depth will serve as a key input parameter for subsequent layered heat source modeling.

[0037] In this embodiment of the invention, based on the skin feature depth, the conductor cross-section is divided into multiple annular sub-regions, the current density integral value of each annular sub-region is calculated, and a layered heat source intensity distribution model is constructed, including:

[0038] Based on the skin depth, the conductor surface to the center is divided into N ring-shaped sub-regions according to a geometric progression. The common ratio of the geometric progression is determined by the attenuation constant of the current density. Geometric progression is an efficient layering method that adapts to the gradient distribution of current density, improving model accuracy through non-uniform partitioning. The partitioning process uses the skin depth as a reference, setting multiple concentric ring-shaped regions from the surface to the center, with the layer thickness increasing geometrically, meaning each layer's thickness is a fixed multiple of the previous layer. (Common ratio) The determination is based on the current density decay constant. The calculation formula is:

[0039] ;

[0040] in, Where is the conductor radius, The adjustment factor is typically set to 2-3, dynamically adjusted according to accuracy requirements and computational efficiency. The number of layers N is usually set to 5-10 to ensure that the spatial characteristics of the current distribution can be fully described within a reasonable computational load. This geometrical hierarchical partitioning strategy ensures finer partitioning in the surface region with a large current density gradient, while the partitioning in the interior region with a relatively gentle change is relatively coarse, thus optimizing the allocation of computational resources.

[0041] The equivalent current value of each annular sub-region is obtained by volume integration of the current density. Equivalent current calculation is a crucial step in obtaining the comprehensive electromagnetic properties of each layer, accurately characterizing the local current distribution through spatial integration. The calculation process performs three-dimensional integration of the current density field within each annular sub-region, considering the spatial variation of the current density within the region to obtain a more accurate equivalent current. The integration method employs Gaussian integration techniques, dynamically adjusting the number and distribution of integration points based on the region shape and current distribution characteristics to ensure integration accuracy. The equivalent current value directly reflects the electromagnetic activity level of each annular sub-region and serves as the fundamental input for subsequent heat source intensity calculations; its accuracy directly affects the precision of thermal field prediction.

[0042] Based on the equivalent current value and resistivity of the conductor material in each annular sub-region, the heat source intensity per unit volume of each annular sub-region is calculated. Calculating the heat source intensity is the core step in converting electromagnetic losses into heat source input, quantifying the heat energy generated by the current based on the Joule heating principle. The calculation process first calculates the total power loss of each annular sub-region based on the equivalent current value; then, considering the region volume, it is converted into heat source intensity per unit volume. The formula for calculating the heat source intensity is:

[0043] ;

[0044] in, For the first Heat source intensity per unit volume of each annular sub-region The resistivity of the conductor material (considering temperature and frequency dependence). This is the equivalent current value for this region. Let be the volume of this region. This calculation takes into account the non-uniform distribution of current and the influence of material properties, and can accurately characterize the complex heat source distribution under high-frequency conditions.

[0045] A hierarchical heat source intensity distribution model is constructed, which includes the spatial location identifiers of each annular sub-region, the heat source intensity per unit volume, and the geometric volume parameters of the annular sub-regions. The construction of the hierarchical heat source model is the final step in integrating multi-layered heat source information, forming a complete description of the spatial distribution of heat sources. The construction process organizes the key parameters of each annular sub-region into a structured model, including spatial location information (inner and outer radii, axial range), thermophysical properties (heat source intensity per unit volume), and geometric information (volume, shape factor). The model adopts a hierarchical data structure, supporting efficient spatial query and interpolation operations, facilitating subsequent mesh mapping processing. This hierarchical heat source model comprehensively describes the non-uniform distribution characteristics of heat sources inside conductors under high-frequency conditions, providing refined heat source input for accurate thermal field prediction and overcoming the limitations of traditional uniform heat source models in high-frequency applications.

[0046] In this embodiment of the invention, the spatial positions of grid nodes in the initial thermal simulation grid data and electromagnetic field simulation data are extracted, the nearest distance and local Jacobian determinant between grid node pairs are calculated, and a grid mapping relationship topology graph is constructed, including:

[0047] Electromagnetic grid node coordinate sets and thermal grid node coordinate sets were extracted from electromagnetic field simulation data and initial thermal simulation grid data, respectively. Coordinate set extraction is the first step in constructing the grid mapping, obtaining precise spatial location information through data parsing. The extraction process first accesses the node information blocks of the electromagnetic field simulation data and thermal simulation grid data, parsing the correspondence between node indices and coordinate values; then, the coordinate data is organized into a structured coordinate set, including node IDs, spatial coordinates, and topological connection information; finally, data preprocessing is performed, including outlier detection and missing value handling, to ensure the integrity and accuracy of the coordinate set. This step provides the foundational data for subsequent spatial mapping analysis, and the coordinate accuracy directly affects the spatial accuracy of the final thermal field prediction.

[0048] Coordinate system consistency verification is performed on the electromagnetic grid node coordinate set and the thermal grid node coordinate set. When coordinate system deviations exist, a rigid body transformation matrix is ​​used to unify them into the same coordinate system. Coordinate system consistency verification is a crucial step in ensuring the comparability of spatial data and resolving the issue of different reference frames being used in simulations of different physics fields. The verification process first analyzes the spatial distribution characteristics of the two coordinate sets, including boundary extent, axial distribution, and geometric center; then, through feature point matching or spatial correlation analysis, it detects whether there are translational, rotational, or scaling differences in the coordinate systems; if a system deviation is detected, a rigid body transformation matrix (containing rotational and translational components) is constructed, and the two coordinate systems are aligned through matrix transformation. The transformation formula is:

[0049] ;

[0050] in, The coordinates are after transformation. Original coordinates For rotation matrix, This coordinate system serves as the translation vector. This unified coordinate system ensures the precise spatial correspondence of different physical field data, a fundamental prerequisite for multiphysics coupling analysis, and eliminates systematic errors caused by coordinate system differences.

[0051] A KD-tree spatial index structure is constructed based on the coordinate set of electromagnetic grid nodes. KD-tree construction is a key technology for accelerating spatial nearest neighbor search, achieving efficient querying through hierarchical spatial partitioning. The construction process adopts a top-down recursive method. First, the coordinate axis with the largest variance is selected as the partition dimension; then, the median coordinate on this dimension is selected as the partition point, dividing the dataset into left and right parts; then, the above process is repeated for each subset, recursively constructing subtrees; finally, a balanced KD-tree structure is formed, where each node contains information about an electromagnetic grid node and the definition of the partition plane. The time complexity of KD-tree construction is O(n log n), where n is the number of electromagnetic grid nodes. This spatial index structure reduces the time complexity of nearest neighbor search from O(n) of the naive method to O(log n), significantly improving the mapping efficiency of large-scale grid systems and is an essential optimization technique for processing complex electronic device models.

[0052] For each hot grid node in the hot grid node coordinate set, the KD-tree spatial index structure is used to query the K nearest electromagnetic grid nodes. The K-nearest neighbor query is the core operation for establishing grid mapping relationships, determining the basic nodes for data transformation through spatial proximity. The query process first inputs the hot grid node coordinates as query points into the KD-tree; then, a priority queue-based nearest neighbor search algorithm is executed, recursively traversing the KD-tree structure and maintaining the currently found K nearest points; finally, the sorted K nearest electromagnetic grid nodes and their distances are returned. The choice of K depends on the grid fineness and mapping accuracy requirements, and is usually set to 4-8, ensuring sufficient spatial interpolation accuracy while avoiding excessive computational overhead. The nearest neighbor query results directly determine the locally dependent nodes of the heat source mapping and are key inputs for constructing accurate heat source distribution.

[0053] Calculate the Euclidean distances from the thermal grid node to each of the K electromagnetic grid nodes, and select the minimum value as the nearest mapping distance. The nearest distance calculation is a direct indicator of the accuracy of spatial mapping, reflecting the tightness of the grid matching. The calculation uses the standard Euclidean distance formula; the distance value intuitively reflects the spatial fit between the two grid systems at that point, and is an important indicator for evaluating mapping quality. For complex geometries, the distribution characteristics of the nearest mapping distance can reveal the spatial variation patterns of grid matching, guiding the optimization and adjustment of subsequent mapping strategies.

[0054] Extract the geometric vertex coordinates of the electromagnetic mesh elements to which each of the K electromagnetic mesh nodes belongs, and construct a geometric mapping function from the standard reference element to the actual physical element. Constructing the geometric mapping function is a key step in achieving high-order accuracy interpolation, capturing the geometric characteristics of the mesh elements through parameterized representation. The construction process first determines the mesh element to which each electromagnetic mesh node belongs and extracts the complete geometric information of the element; then, based on the element type (e.g., tetrahedral, hexahedral, etc.), an appropriate basis function system is selected; finally, a mapping function from the standard reference element (usually a [-1,1] cube or unit tetrahedron) to the actual physical element is constructed. For common tetrahedral elements, the mapping function adopts a linear form:

[0055] ;

[0056] in, For the parameter coordinates in the reference element, For shape functions, The coordinates of the physical unit vertex are given. For higher-order units, corresponding higher-order polynomial basis functions are used. This parameterized mapping provides the mathematical foundation for subsequent Jacobian calculations and numerical integration, and is a key technology for achieving accurate geometric representation and high-precision interpolation.

[0057] The Jacobian matrix of the geometric mapping function at the projected location of the hot mesh node is calculated, and its determinant is denoted as the local Jacobian determinant. Jacobian matrix calculation is a core step in quantifying mesh deformation characteristics, evaluating the local properties of mesh mapping through differential geometric analysis. The calculation process first determines the projected location of the hot mesh node in the reference cell, achieved through inverse mapping or nearest-point projection; then, the partial derivative matrix of the geometric mapping function at that point, i.e., the Jacobian matrix, is calculated. Finally, calculate the determinant of the matrix. .

[0058] Determinant It intuitively reflects the proportion of local volume change during the mapping process and is an important indicator for evaluating the degree of mesh distortion and mapping quality. Larger... The value indicates that the grid is relatively sparse in this area, and a smaller value indicates that the grid is relatively small. The value indicates that the grid is relatively dense, and this characteristic directly affects the weight allocation and accuracy control of numerical interpolation.

[0059] A topological graph of mesh mapping relationships is constructed, with hot mesh nodes as vertices and directed edges representing the mapping relationships between hot mesh nodes and their corresponding K electromagnetic mesh nodes. The weight of each directed edge is determined by the product of the reciprocal of the nearest mapping distance and the local Jacobian determinant. Topological graph construction is the final step in integrating spatial relationships, expressing the mapping dependencies between meshes through a graph structure. The construction process first creates a graph structure with hot mesh nodes as vertices; then, directed edges pointing to their K nearest neighbor electromagnetic mesh nodes are added to each hot mesh node; finally, the weight of each edge is calculated and assigned using the following formula:

[0060] ;

[0061] in, To obtain from hot grid nodes To electromagnetic grid node edge weights, This is the corresponding Euclidean distance. for The Jacobian determinant at each point. This weighting design comprehensively considers spatial distance and mesh deformation characteristics; edges with closer distances and larger Jacobian determinants receive higher weights, more accurately reflecting the physical dependencies of the mapping. The completed topology graph comprehensively describes the spatial correspondence and mapping weight distribution of the two mesh systems, providing systematic data structure support for subsequent heat source distribution mapping and ensuring the accuracy and efficiency of multiphysics data conversion.

[0062] In this embodiment of the invention, based on the mesh mapping relationship topology diagram, the hierarchical heat source intensity distribution model is mapped to the initial thermal simulation mesh data, and the total power deviation before and after mapping is calculated, including:

[0063] For each thermal grid cell in the initial thermal simulation mesh data, the nearest electromagnetic grid nodes to that thermal grid cell's centroid are found in the mesh mapping topology graph. Nearest neighbor search is a fundamental operation in mesh mapping, establishing correspondences between different mesh systems through spatial relationships. The search process first calculates the centroid coordinates of the thermal grid cell as spatial location points; then, a nearest neighbor search is performed in the mapping topology graph to identify the nearest electromagnetic grid nodes. The search algorithm employs spatial indexing techniques such as KD-trees or octrees, significantly improving the search efficiency for large-scale mesh systems. The number of nearest neighbor nodes is typically set to 4-8, dynamically adjusted according to the mesh density ratio and accuracy requirements. This spatial nearest neighbor-based mapping strategy maintains physical spatial continuity while adapting to structural differences between different mesh systems, providing a reliable foundation for subsequent data interpolation.

[0064] Based on the hierarchical heat source intensity distribution model data of multiple electromagnetic grid nodes, the preliminary heat source intensity of the thermal grid cell is calculated using the inverse distance weighted interpolation method. Inverse distance weighted interpolation is the core technology for data transformation between grids, achieving spatial continuity mapping of data through distance weight allocation. The interpolation process first calculates the Euclidean distance from the centroid of the thermal grid cell to each nearest neighbor electromagnetic node; then, it calculates the normalized weights based on the reciprocal of the distances; finally, it weights the heat source intensity of the annular sub-region to which each node belongs, and then calculates the preliminary heat source intensity of the thermal grid cell by weighted averaging. The calculation formula is:

[0065] ;

[0066] in, The initial heat source intensity for the thermal grid cell. For the first The weights of the nearest neighbor electromagnetic nodes, This represents the heat source intensity per unit volume of the annular sub-region to which the node belongs. When the centroid is very close to a node (less than a preset threshold), the heat source intensity value of that node is used directly to avoid interpolation singularity issues. This spatial distance-based interpolation method ensures both the smoothness of the data mapping and the preservation of the physical characteristics of the heat source distribution, making it suitable for complex geometries and irregular grid systems.

[0067] The initial heating power of a thermal grid cell is obtained by multiplying the initial heat source intensity by the volume of that cell. This heating power calculation is a necessary step in converting intensity density into absolute power, taking into account the volume differences within the grid cells. The calculation process directly multiplies the heat source intensity per unit volume by the volume of the grid cell to obtain the total heating power of that cell. This calculation considers the non-uniformity of the grid division, ensuring that regions with the same heat source intensity generate heat proportional to volume in grid cells of different sizes, thus maintaining the conservation of physical quantities. Heating power is fundamental data for energy balance analysis, and its accurate calculation is crucial for subsequent power deviation analysis and energy compensation.

[0068] The initial heating power of all thermal grid cells is summed and recorded as the total power after mapping. This total power summation is a crucial step in evaluating the global conservation of energy mapping, verifying the accuracy of energy conversion through aggregate calculation. The summation process traverses all thermal grid cells, adding the initial heating power of each cell to obtain the total input power of the entire thermal field system. This global summation operation not only provides a benchmark value for energy conservation verification but also reflects potential systematic errors in the grid mapping process, providing a quantitative basis for subsequent energy compensation.

[0069] The total heating power of all annular sub-regions in the layered heat source intensity distribution model is calculated and denoted as the total power before mapping. This pre-mapping power calculation is a crucial step in establishing an energy reference baseline, providing the total power loss in electromagnetic analysis. The calculation process traverses all annular sub-regions in the layered heat source model, multiplying the heat source intensity per unit volume of each region by the corresponding volume, and then summing the results to obtain the total power. This reference value represents the theoretical total loss determined by electromagnetic analysis and serves as the benchmark for energy conservation assessment; its accuracy directly affects the precision of subsequent energy compensation.

[0070] The difference between the total power before and after mapping is calculated and denoted as the total power deviation. Power deviation calculation is a direct method for quantifying the accuracy of energy mapping, and the conservation error is evaluated through difference analysis. The calculation uses simple algebraic differences:

[0071] ;

[0072] in, For total power deviation, The total power before mapping. This represents the total power after mapping. The deviation can be positive (energy loss during mapping) or negative (energy gain during mapping), and its absolute value and relative proportion directly reflect the energy conservation performance of the mesh mapping process. When the relative deviation exceeds a preset threshold (usually 1%-5%), energy compensation is required to ensure the physical accuracy of the thermal field calculation. Power deviation analysis is a crucial step in ensuring energy consistency in multiphysics coupling, providing a quantitative basis for subsequent compensation and correction.

[0073] In this embodiment of the invention, an energy conservation compensation source term is introduced at the boundary nodes of the grid mapping topology based on the total power deviation to generate a corrected heat source distribution field, including:

[0074] In the mesh mapping topology diagram, thermal mesh cells located on the conductor surface are identified and marked as boundary thermal mesh cells. Boundary cell identification is a preliminary step for energy compensation, locating suitable compensation areas based on spatial location characteristics. The identification process first analyzes the spatial location characteristics of thermal mesh cells in the mesh topology diagram; then, combined with the conductor geometry model, it determines whether the cell is located near the conductor surface; finally, cells that meet the conditions are marked as boundary cells. The judgment criterion is usually based on the shortest distance from the cell to the surface. When the distance is less than a preset threshold (usually a certain proportion of the average mesh size), the cell is considered to be located in the surface region. As the main carrier of energy compensation, boundary cells can achieve global energy conservation through surface adjustment while maintaining the internal heat source distribution characteristics. This strategy conforms to actual physical phenomena because electromagnetic losses themselves have a concentrated characteristic in the surface region.

[0075] The total number and volume of boundary thermal mesh cells are statistically analyzed. Boundary cell statistics are a necessary step in determining the compensation basis, providing a basis for subsequent energy allocation through quantity and volume calculations. The statistical process traverses all thermal mesh cells marked as boundary cells, accumulating the number and summing the volumes. These statistical statistics directly reflect the mesh resources available for compensation, providing a quantitative basis for the rational allocation of compensation. The number and distribution characteristics of boundary cells also reflect the complexity of the geometric model and the level of detail in the mesh generation, indirectly affecting the selection of compensation strategies.

[0076] The total power deviation is distributed according to the volume ratio of the boundary thermal grid cells to obtain the power compensation amount for each boundary thermal grid cell. Power compensation allocation is the core operation for achieving energy conservation, and proportional allocation ensures the rationality of the compensation. The allocation process is based on the volume ratio of the boundary cells, dividing the total power deviation according to the proportion of each cell to the total boundary volume, ensuring that each cell receives a compensation amount proportional to its volume. This volume-based allocation method considers both the non-uniformity of the mesh division and the relative uniformity of the compensated heat source intensity, avoiding numerical problems that may be caused by local overcompensation. Power compensation allocation is a key step in energy conservation correction, and its rationality directly affects the accuracy and stability of the thermal field calculation.

[0077] The compensated heat source intensity of each boundary thermal grid cell is obtained by dividing the power compensation by the volume of that cell. Calculating the compensation intensity is a necessary step in converting power compensation into heat source intensity, ensuring consistency with the original heat source format. The calculation process directly divides the allocated power compensation by the cell volume, converting it into heat source intensity per unit volume. This conversion ensures that the compensation terms have the same physical dimensions and numerical range as the original heat source distribution, facilitating direct superposition and unified processing. The calculation of the compensated heat source intensity completes the conversion from global power deviation to local heat source adjustment, providing direct input for the final thermal field correction.

[0078] The initial heat source intensity of each thermal grid cell is added to its corresponding compensated heat source intensity to generate a corrected heat source distribution field, which satisfies the global energy conservation constraint. Heat source field correction is the final step in energy compensation, achieving energy calibration of the thermal field model through source term superposition. The correction process adds the initial heat source intensity of each boundary thermal grid cell to the calculated compensated heat source intensity, forming the corrected heat source intensity; for non-boundary cells, the initial heat source intensity remains unchanged. This boundary compensation strategy achieves global energy conservation while maintaining the original characteristics of the internal heat source distribution, and also conforms to the characteristic of concentrated surface losses in actual physical phenomena. The corrected heat source distribution field satisfies the energy conservation constraint, and the total power remains consistent with the electromagnetic analysis results, providing accurate heat source input for subsequent thermal field calculations and ensuring the physical consistency and computational accuracy of multiphysics coupling analysis.

[0079] In this embodiment of the invention, a multi-timescale coupled transfer function from electromagnetic loss to temperature rise is constructed. This multi-timescale coupled transfer function includes the time constants of the thermal diffusion characteristics at each level of the hierarchical heat source intensity distribution model, including:

[0080] For each annular sub-region in the layered heat source intensity distribution model, the squared ratio of the characteristic size of the annular sub-region to the material's thermal diffusivity is calculated and denoted as the characteristic time constant of the annular sub-region. Calculating the characteristic time constant is a crucial step in quantifying the dynamic characteristics of thermal diffusion, determining the thermal response rate through material properties and geometric dimensions. The calculation process first determines the characteristic size of the annular sub-region, typically using radial thickness or equivalent diameter; then, combining the material's thermal diffusivity, the characteristic time is derived and calculated using the classical thermal diffusion equation. The calculation formula is:

[0081] ;

[0082] in, For the first The thermal diffusion characteristic time constant of each annular sub-region The characteristic size of this region, The thermal diffusivity of the material ( , Thermal conductivity, For density, (Specific heat capacity). The characteristic time constant intuitively reflects the characteristic time required for heat to propagate in the region and is a fundamental parameter for constructing dynamic thermal response models. Due to the skin effect, the heat source distribution varies at different depths, and the thermal diffusion characteristics of each annular sub-region also differ significantly. By calculating the characteristic time constant layer by layer, this multi-scale thermal response characteristic can be accurately captured.

[0083] Based on the thermal diffusion characteristic time constant of each annular sub-region, an equivalent circuit of a multi-order RC thermal network is constructed, where each RC node corresponds to an annular sub-region. RC thermal network construction is an effective method for simulating the dynamic process of thermal diffusion, characterizing the dynamic response of the thermal system through a circuit model. The construction process maps the dynamic behavior of the thermal system to a circuit model, where heat sources correspond to current sources, temperature to voltage, thermal resistance to resistance, and thermal capacitance to capacitance. Each annular sub-region is modeled as an RC node, and the values ​​of R and C are calculated based on the principles of thermal resistance and thermal capacitance.

[0084] ;

[0085] ;

[0086] in, For the first Thermal resistance of each annular sub-region The length of the heat conduction path. Thermal conductivity, The area through which the heat flows; For the first The heat capacity of each annular sub-region For density, For specific heat capacity, The nodes are connected by thermal coupling to form a complete network structure. This multi-order RC network model can effectively capture the thermal response characteristics of heat sources at different depths, reflecting the multi-timescale effects of heat diffusion.

[0087] The Laplace transform is applied to the equivalent circuit of a multi-order RC thermal network to obtain its frequency domain transfer function expression. The Laplace transform is a standard method for converting a time-domain model into a frequency-domain expression, facilitating system analysis and mathematical processing. The transformation process first establishes the state equations of the RC network, describing the dynamic changes in temperature (voltage) at each node; then, the Laplace transform is applied to the equations to convert the time-domain differential equations into algebraic equations; finally, the frequency domain transfer function from input (heat source) to output (temperature) is solved. For multi-order RC networks, the transfer function is typically expressed in rational fraction form:

[0088] ;

[0089] in, For transfer functions, For Laplace variables, For the first Gain coefficients of each annular sub-region This is the characteristic time constant. This frequency domain representation intuitively reflects the dynamic characteristics and frequency response of the system, providing a mathematical basis for subsequent phase correction and time domain transformation.

[0090] A phase correction term is introduced into the frequency domain transfer function expression. This term is determined based on the ratio of the skin depth to the thermal wave penetration depth. Phase correction is a crucial step in improving the accuracy of the dynamic response, considering the coupled phase relationship between the electromagnetic and thermal physical processes. The correction process is based on a comparison of the skin depth and the thermal wave penetration depth, analyzing the spatial scale relationship between the electromagnetic and thermal fields, and introducing appropriate phase adjustments. The correction term is typically expressed in complex exponential form. The phase angle Calculated based on the depth ratio:

[0091] ;

[0092] in, For skin feature depth, The thermal wave penetration depth is defined by the frequency characteristics of the heat diffusion equation. This phase correction takes into account the difference in propagation velocities between the electromagnetic and thermal fields, compensates for the time delay effect between the two physical fields, and improves the accuracy of transient response prediction, especially under high-frequency alternating heat source conditions.

[0093] The inverse Laplace transform (ILT) is applied to the frequency domain transfer function expression to obtain the time-domain form of the multi-timescale coupled transfer function. The ILT is a necessary step in converting the frequency domain model back to its time-domain expression, yielding a mathematical model directly usable for transient calculations. The transformation process applies the ILT to the frequency domain transfer function, which includes phase corrections, converting the s-domain function into a time function. For typical multi-order RC networks, the time-domain transfer function is usually expressed as a linear combination of exponential terms:

[0094] ;

[0095] in, For time-domain transfer function, For time, This is the gain coefficient. The characteristic time constant, For unit step function, The time delay is determined by the phase correction term. This multi-timescale transfer function form can accurately describe the dynamic characteristics of the thermal response, reflect the thermal diffusion process of heat sources at different depths, provide an efficient and accurate mathematical tool for transient thermal field prediction, and effectively solve the problem of dynamic response deviation caused by neglecting multi-timescale effects in traditional thermal analysis.

[0096] In this embodiment of the invention, a multi-timescale coupled transfer function is applied to the corrected heat source distribution field to generate transient thermal field prediction results, including:

[0097] The corrected heat source distribution field is grouped according to the annular sub-regions to which the thermal grid cells belong, and the corresponding heat source time series for each annular sub-region is obtained. Heat source grouping is a fundamental step in multi-timescale analysis, achieving hierarchical processing through spatial region division. The grouping process first determines the annular sub-region to which each thermal grid cell belongs through spatial location comparison and region lookup; then, the heat source intensities of thermal grid cells within the same annular sub-region are grouped together to form a comprehensive heat source for that region; finally, the heat source time series for each annular sub-region is constructed based on its time-varying characteristics. This spatially hierarchical heat source organization method directly corresponds to the structure of the hierarchical heat source model and the multi-timescale transfer function, ensuring the consistency between model input and calculation method, and laying the foundation for accurately capturing the dynamic response characteristics of heat sources at different depths.

[0098] The heat source time series of each annular sub-region is applied with a sub-transfer function from the multi-timescale coupled transfer function corresponding to that sub-region to obtain the temperature rise response time series of each annular sub-region. Transfer function application is the core step in calculating dynamic temperature rise, achieving the response calculation from heat source to temperature through mathematical convolution. The application process processes each annular sub-region separately, performing convolution operations between the heat source time series of that region and the corresponding sub-transfer function to calculate the temperature rise response time series of that region. The convolution calculation is accelerated using Fast Fourier Transform (FFT) technology, converting time-domain convolution into frequency-domain multiplication, significantly improving computational efficiency. The parameters of the sub-transfer function (such as characteristic time constants and time delays) are customized according to the physical characteristics and location of the region to ensure the physical accuracy of the temperature rise calculation. This regional processing method fully considers the differences in thermal response of heat sources at different depths and can accurately capture the multi-timescale characteristics of heat diffusion.

[0099] The temperature rise response time series of each annular sub-region is mapped back to the thermal grid cells of the initial thermal simulation grid data according to their spatial location. Temperature rise mapping is a key step in restoring the spatial temperature distribution, and the grouping results are restored to grid data through back projection. The mapping process first determines the annular sub-region to which each thermal grid cell belongs; then, the temperature rise response time series of that region is extracted; finally, the time series data is assigned to the corresponding grid cells, completing the conversion from region-level response to cell-level temperature. For grid cells located at the region boundary, a distance-based weighted interpolation method is used to ensure the spatial continuity of the temperature distribution. This back mapping process maintains the spatial integrity of the physical field while preserving the dynamic characteristics of multi-timescale analysis, providing spatially structured temperature rise data for subsequent temperature overlay.

[0100] The temperature rise responses received by each thermal grid cell from different annular subregions are time-domain superimposed to obtain the total temperature rise of that thermal grid cell. Temperature rise superposition is a key step in synthesizing multi-source thermal effects, and the overall temperature response is calculated using the principle of linear superposition. The superposition process treats each thermal grid cell individually, adding the time series of temperature rise responses received by that cell from each annular subregion point-to-point to obtain the comprehensive temperature rise effect. This superposition is based on the linear assumption of the heat conduction system and is applicable to most engineering applications. For materials with strong nonlinear temperature dependence, piecewise linearization is used to ensure the accuracy of the superposition results. Temperature rise superposition is the final calculation step in the synthesis of multi-source thermal responses, and its results directly reflect the dynamic temperature change characteristics under complex heat source distributions.

[0101] By combining the ambient temperature boundary conditions of the electronic device, the total temperature rise of each thermal grid cell is converted into absolute temperature to generate transient thermal field prediction results. Boundary condition integration is the final step in completing the thermal field calculation, achieving the conversion from relative temperature rise to absolute temperature by adding a reference temperature. The integration process first determines the ambient temperature conditions of the electronic device, including ambient temperature, cooling conditions, and thermal boundary constraints; then, these conditions are combined with the calculated temperature rise data to convert them into absolute temperature values; finally, complete transient thermal field data is generated, including temperature distribution information in both time and space dimensions. For complex boundary conditions, a thermal boundary layer model is used for refinement to improve the calculation accuracy of the boundary region. The transient thermal field prediction results are presented in time series form, intuitively demonstrating the dynamic evolution of temperature distribution, providing accurate technical basis for thermal management and reliability analysis of electronic devices, and is particularly suitable for thermal analysis of electronic devices under high-frequency operating conditions.

[0102] In this embodiment of the invention, the characteristic frequency components of the temperature response in the transient thermal field prediction results are extracted, and phase deviation analysis is performed between them and the excitation frequency of the electromagnetic field simulation data, including:

[0103] Temperature monitoring points on the conductor surface are selected from the transient thermal field prediction results, and their temperature time series are extracted. Monitoring point selection is a preliminary step in frequency analysis, obtaining representative data through sampling at key locations. The selection process first identifies key locations with significant temperature gradients on the conductor surface; then, temperature monitoring points are set up, and the complete temperature time response is recorded; finally, standardized time series data is generated for subsequent spectral analysis. Monitoring points are typically selected in hotspot regions or on surfaces with concentrated current density, as these areas are most sensitive to electromagnetic excitation and provide temperature signals with high signal-to-noise ratios. Selecting multiple monitoring points distributed in different locations provides more comprehensive temperature response characteristics, enhancing the representativeness and reliability of the analysis results.

[0104] A short-time Fourier transform (SFT) is performed on the temperature time series to obtain the time-frequency spectral distribution of the temperature response. The SFT is an effective tool for analyzing the spectral characteristics of non-stationary signals, providing a two-dimensional representation of the signal in both time and frequency dimensions. The transform process first selects an appropriate window function (usually a Hanning or Hamming window); then, the temperature time series is segmented; finally, a Fourier transform is performed on the data within each time window to obtain the time-varying spectrum. The time-frequency spectral distribution visually demonstrates the time evolution characteristics of the frequency components in the temperature response, enabling the identification of non-stationary characteristics and frequency modulation phenomena, and providing a spectral perspective for in-depth analysis of the relationship between electromagnetic excitation and thermal response.

[0105] In the time-frequency spectrum distribution, the frequency component with the largest amplitude is identified and denoted as the dominant characteristic frequency. Characteristic frequency identification is the core step in spectrum analysis, determining the main periodic components of the signal through peak detection. The identification process first performs peak detection on the time-frequency spectrum to find local maxima of the spectral amplitude; then, it sorts the frequencies according to their amplitude to determine the dominant frequency components; finally, it extracts the amplitude and phase information of these frequency components to form a set of characteristic frequencies. The dominant characteristic frequency directly reflects the main periodic characteristics of the temperature response and usually has a direct correspondence with the electromagnetic excitation frequency, making it a key parameter for analyzing the dynamic characteristics of the thermal response. For complex multi-frequency responses, there may be multiple significant characteristic frequencies, requiring comprehensive consideration of their relative amplitudes and physical meanings.

[0106] The excitation frequency and its harmonics are extracted from electromagnetic field simulation data and frequency matched with the dominant characteristic frequency. Frequency matching is a crucial step in establishing the correlation between electromagnetic excitation and thermal response, determining the correspondence between response and excitation through frequency comparison. The matching process first extracts the fundamental excitation frequency from the electromagnetic field simulation data; then calculates possible harmonics and combination frequencies; and finally compares them with the dominant characteristic frequency of the temperature response to determine the matching relationship. The matching criteria are based on the frequency ratio and allowable error range, typically requiring a relative error of less than 1-5% to be considered a valid match. This frequency matching analysis can reveal the frequency transfer characteristics from electromagnetic excitation to thermal response, identify linear response components and nonlinear modulation effects, and provide a frequency domain perspective for a deeper understanding of thermoelectric coupling mechanisms.

[0107] For successfully matched frequency components, the phase difference between the temperature time series and the electromagnetic loss power time series at that frequency is calculated and denoted as the phase deviation. Phase deviation calculation is a crucial step in evaluating the dynamic response delay, quantifying the time relationship through phase comparison. The calculation process first extracts complex spectral components from the temperature and electromagnetic loss time series at the matched frequency; then, it calculates the phase angle between the two complex phasors; finally, it obtains the phase difference value, thus yielding the phase deviation. The phase deviation directly reflects the time delay of the thermal response lagging behind the electromagnetic excitation and is an important indicator for evaluating the dynamic characteristics of heat transfer. Due to the inertial characteristics of the heat diffusion process, the temperature response typically lags behind the power excitation. The magnitude of the phase deviation reflects the response speed of the thermal system and is key reference data for optimizing the dynamic characteristics of the heat source model.

[0108] The spatial distribution pattern of phase deviation in different annular sub-regions is statistically analyzed to identify anomalous regions where the phase deviation exceeds a preset phase threshold. Spatial distribution analysis is a crucial step in comprehensively evaluating model accuracy, identifying model limitations through spatial patterns of phase deviation. The analysis process first maps phase deviation data to spatial annular sub-regions; then, it statistically analyzes the average and standard deviation of the deviation in each region; finally, it identifies anomalous regions by comparing them with preset thresholds. The anomaly judgment criteria are usually based on statistical significance; when the deviation of a region exceeds two standard deviations of the average or a preset absolute threshold (usually 30°-45°), it is marked as an anomalous region. These anomalous regions usually indicate that the current heat source model has inaccuracies in the temporal dynamic characteristics of that region, requiring correction through phase adjustment. Spatial distribution analysis not only provides a global perspective for model evaluation but also provides spatial location basis for subsequent targeted optimization.

[0109] In this embodiment of the invention, based on the phase deviation analysis results, the heat source delay phase of each annular sub-region in the layered heat source intensity distribution model is dynamically adjusted to update the transient thermal field prediction results, including:

[0110] For each annular sub-region, the corresponding phase deviation is extracted. Regional phase extraction is a preliminary step in phase adjustment, obtaining representative deviation values ​​through regional aggregation. The extraction process first determines the spatial extent of the annular sub-region; then, it collects phase deviation data from all monitoring points within that region; finally, it calculates a representative value using statistical methods (such as median or weighted average) as the characteristic phase deviation of the region. When multiple monitoring points exist within a region, different weights are assigned based on noise level and signal quality assessment to ensure the representativeness and robustness of the extraction results. Regional phase extraction is a crucial step in transforming point-based monitoring data into regional characteristics, providing fundamental parameters for spatial layering in subsequent phase adjustment.

[0111] Phase deviation is converted into a time delay, which is equal to the phase deviation divided by the angular frequency of the corresponding frequency. Phase-time conversion is a necessary step for time-domain adjustment, obtaining the absolute time value through phase angle conversion. The conversion uses a standard formula:

[0112] ;

[0113] in, This is the amount of time delay. Phase deviation (radians). This represents the angular frequency corresponding to the given frequency. This transformation concretizes the abstract phase relationship into a time delay, facilitating intuitive understanding and operational implementation. The amount of time delay directly reflects the amount of time advancement or retardation required to adjust the heat source model, providing a quantitative basis for phase adjustment of the heat source time series. For multi-frequency systems, the dominant frequency is usually selected for transformation, or a weighted average is used to obtain the comprehensive time delay.

[0114] In the hierarchical heat source intensity distribution model, a time-shift operation with a time delay is applied to the heat source time series of the annular sub-region. The time-shift operation is the core step in achieving phase adjustment, dynamically correcting the phase through time series rearrangement. The operation shifts the heat source time series in the time domain, adjusting the time points forward or backward based on the calculated time delay to advance or delay the phase. For discrete time series, the time shift is handled by linear interpolation or spline interpolation to handle non-integer time points, ensuring a smooth and continuous operation. The time-shifted heat source time series, while maintaining unchanged waveform and amplitude, adjusts the phase characteristics, corrects the temporal dynamics of the heat source model, and provides more accurate input for improving temperature response prediction.

[0115] The time series of the heat source after the time shift operation is applied is resubmitted into the multi-timescale coupled transfer function to recalculate the temperature rise response of the annular subregion. Temperature rise recalculation is a direct step to verify the adjustment effect, obtaining an improved temperature output through updated heat source input. The calculation process uses the adjusted heat source time series as new input, applies the multi-timescale coupled transfer function for convolution calculation, and generates an updated temperature rise response. Compared to the initial calculation, this step only updates the temporal characteristics of the heat source input; the transfer function itself remains unchanged, ensuring the consistency of the physical model. Temperature rise recalculation directly verifies the effect of phase adjustment, provides the optimized dynamic temperature response, and provides regional-level foundational data for subsequent global temperature updates.

[0116] After completing the above adjustments for all annular sub-regions, the temperature rise time-domain superposition process is re-executed to generate updated transient thermal field prediction results. Temperature rise re-superposition is the final step in completing the global update, generating the final optimized result through comprehensive integration. The superposition process comprehensively considers the updated temperature rise responses of all annular sub-regions, recalculating the total temperature rise of each thermal grid cell according to the same spatial mapping and superposition principles, and converting it to absolute temperature based on boundary conditions to form a complete updated transient thermal field. This updated result fully reflects the comprehensive effect of phase adjustment, improving the dynamic accuracy of temperature response by optimizing the temporal characteristics of heat sources in each region, especially significantly improving the temporal accuracy of thermal field prediction under high-frequency alternating conditions where phase relationship is critical. The final transient thermal field prediction results achieve precise coupling between the electromagnetic field and the thermal field in both time and space dimensions, providing a reliable technical basis for the thermal design and management of electronic devices.

[0117] In this embodiment of the invention, based on the hierarchical heat source intensity distribution model data of multiple electromagnetic grid nodes, the preliminary heat source intensity of the thermal grid cell is calculated using the inverse distance weighted interpolation method, including:

[0118] The Euclidean distance from the centroid of the thermal grid cell to each electromagnetic grid node is calculated and denoted as the node distance. Node distance calculation is a fundamental step in spatial interpolation, determining the strength of spatial association through geometric metrics. The calculation uses the standard Euclidean distance formula. Distance calculation directly reflects the geometric proximity between spatial points, providing a basic metric for subsequent weight calculation. For large-scale grid systems, spatial indexing techniques (such as KD-trees or octrees) are used to accelerate nearest neighbor search and distance calculation, significantly improving algorithm efficiency. Node distance is a direct input to weight calculation, and its accuracy directly affects the spatial continuity and physical plausibility of the interpolation results.

[0119] The distances to each node are reciprocalized and normalized to obtain the distance weights for each electromagnetic grid node. Weight calculation is the core step of inverse distance interpolation, determining the node contribution ratio through distance transformation. The calculation process first takes the reciprocal of the node distances to form initial weights, reflecting the interpolation principle that closer distances have greater influence; then, all weights are normalized to ensure the total weight sum is 1, satisfying the completeness requirement of interpolation. The calculation formula is:

[0120] ;

[0121] in, For the first Normalized weights of each node, The distance from the node to the centroid is calculated, and the summation range includes all considered nearest neighbor nodes. The power of the reciprocal of the distance can be adjusted according to the interpolation smoothness requirements, usually set to 1 or 2. The higher the power, the more prominent the influence of nearby points, and the more localized the interpolation result. The normalized weight directly reflects the relative importance of each node in the interpolation calculation and is the core control parameter for data interpolation.

[0122] The heat source intensity per unit volume of each electromagnetic grid node within its annular sub-region is extracted. Heat source intensity extraction is a crucial step in obtaining interpolation source data, and the physical quantity value is obtained through model lookup. The extraction process first determines the annular sub-region to which each electromagnetic grid node belongs, achieved through spatial location comparison and region lookup. Then, the heat source intensity per unit volume of that region is read from the hierarchical heat source intensity distribution model as the source data for interpolation calculations. For accurate region attribution determination, the spatial location of the nodes and the geometric boundaries of the annular sub-regions are considered to ensure the accuracy of the extracted data. Heat source intensity is the core physical quantity for interpolation calculations, and its accurate extraction is a fundamental condition for ensuring the physical rationality of the thermal field.

[0123] The heat source intensity per unit volume at each electromagnetic grid node is multiplied by its corresponding distance weight, and then summed using a weighted average to obtain the preliminary heat source intensity. This weighted summation is the final step in spatial interpolation, generating an estimate of the physical quantities at the target point through comprehensive calculation. The calculation uses the standard weighted average formula:

[0124] ;

[0125] in, The initial heat source intensity for the thermal grid cell. For the first Distance weights of each electromagnetic grid node This represents the heat source intensity per unit volume of the annular sub-region to which the node belongs. This spatial distance-based weighted averaging method considers the contributions of all neighboring nodes and assigns different influence weights according to spatial distance, generating a smooth and continuous heat source distribution field while maintaining the spatial continuity of the physical field. The weighted summation is the final calculation step of the inverse distance interpolation, and its result is directly used for subsequent thermal field calculations.

[0126] When the distance between the centroid of a thermal grid cell and a node of an electromagnetic grid is less than a preset distance threshold, the heat source intensity per unit volume of the annular sub-region to which the electromagnetic grid node belongs is directly assigned as the preliminary heat source intensity, avoiding interpolation singularities. Singularity handling is a special step to ensure interpolation stability, using a threshold to handle extreme cases. The process first checks the distance from the centroid to each node to determine if there are extremely close points less than the threshold. If so, the regular interpolation calculation is skipped, and the heat source intensity of the region to which the close point belongs is directly used as the result. The distance threshold is usually set to a small fraction of the average grid size (e.g., 10%) to ensure that the direct assignment mechanism is triggered only in extremely close cases. This singularity handling avoids the numerical instability of weight calculations when node distances are close to zero, and is also physically more reasonable because the physical characteristics of extremely close points are usually the most representative. Singularity handling is an important supplement to improving the robustness of the interpolation algorithm, ensuring stable and reasonable interpolation results under various grid configurations.

[0127] This invention achieves rapid and accurate prediction of the thermal field of electronic devices through cross-physics coupling analysis of electromagnetic and thermal fields, refined modeling of the skin effect, construction of multi-layer heat source distribution, mesh mapping and energy compensation, construction of multi-timescale transfer functions, and dynamic phase adjustment. The multi-physics coupling method of this invention can accurately simulate electromagnetic-thermal coupling effects under high-frequency conditions and effectively predict the thermal field distribution of complex electronic devices.

[0128] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0129] It should be noted that all formulas in this manual are calculated by removing dimensions and taking their numerical values. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0130] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A method for fast prediction of thermal field of electronic equipment facing strong coupling of multiple physical fields, characterized in that, The method comprises the following steps: obtain electromagnetic field simulation data and initial thermal simulation grid data of the electronic device at the working frequency, the electromagnetic field simulation data including the spatial distribution of current density inside the conductor; perform gradient analysis on the spatial distribution of current density along the radial direction of the conductor to identify the radial depth at which the current density decays to a preset proportion of the surface peak value, which is recorded as the skin characteristic depth; based on the skin characteristic depth, divide the conductor cross-section into multiple layers of annular sub-regions, calculate the current density integral value of each annular sub-region, and construct a layered heat source intensity distribution model; extract the grid node spatial position of the initial thermal simulation grid data and the electromagnetic field simulation data, calculate the nearest distance and local Jacobian determinant between the grid node pairs, and construct a grid mapping relationship topology graph; based on the grid mapping relationship topology graph, map the layered heat source intensity distribution model to the initial thermal simulation grid data, and calculate the total power deviation before and after mapping; introduce an energy conservation compensation source term at the boundary nodes of the grid mapping relationship topology graph according to the total power deviation to generate a corrected heat source distribution field; construct a multi-time scale coupling transfer function of electromagnetic loss to temperature rise, apply the multi-time scale coupling transfer function to the corrected heat source distribution field to generate a transient thermal field prediction result; extract the characteristic frequency component of the temperature response in the transient thermal field prediction result, and perform phase deviation analysis with the excitation frequency of the electromagnetic field simulation data; based on the phase deviation analysis result, dynamically adjust the heat source delay phase of each annular sub-region in the layered heat source intensity distribution model, and update the transient thermal field prediction result.

2. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The method comprises the following steps: starting from the surface node of the conductor, extract a current density sampling sequence along the radial direction of the conductor; perform exponential fitting on the current density sampling sequence to obtain the decay constant of the current density; calculate the radial position corresponding to 1 / e times the surface peak value of the current density decay constant, which is recorded as the initial skin depth; centering on the initial skin depth, extend a radial distance of a preset proportion inward and outward to form a skin effect evaluation window; calculate the proportion of the integral of the current density in the skin effect evaluation window to the total current of the conductor, which is recorded as the skin concentration degree; when the skin concentration degree is greater than a preset concentration threshold, mark the outer boundary depth of the skin effect evaluation window as the skin characteristic depth.

3. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The method comprises the following steps: based on the skin characteristic depth, divide the conductor cross-section into multiple layers of annular sub-regions, calculate the current density integral value of each annular sub-region, and construct a layered heat source intensity distribution model, which comprises the following steps: divide the conductor surface layer to the center into N layers of annular sub-regions in geometric progression; perform volume integration on the current density in each annular sub-region to obtain the equivalent current value of the annular sub-region; calculate the unit volume heat source intensity of each annular sub-region according to the equivalent current value of each annular sub-region and the resistivity of the conductor material; The layered heat source intensity distribution model comprises spatial position identification of each annular sub-region, unit volume heat source intensity, and geometric volume parameters of the annular sub-region.

4. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The mapping of the layered heat source intensity distribution model to the initial thermal simulation grid data based on the grid mapping relationship topology and the calculation of the total power deviation before and after the mapping comprise: Each heat grid cell in the initial thermal simulation grid data is traversed, and a plurality of electromagnetic grid nodes closest to the center of the heat grid cell are searched in the grid mapping relationship topology; Based on the layered heat source intensity distribution model data of the plurality of electromagnetic grid nodes, a preliminary heat source intensity of the heat grid cell is calculated by using an inverse distance weighted interpolation method; The preliminary heat source intensity is multiplied by the volume of the heat grid cell to obtain a preliminary heat generation power of the heat grid cell; The preliminary heat generation powers of all heat grid cells are accumulated and recorded as a total power after the mapping; The total heat generation power of all annular sub-regions in the layered heat source intensity distribution model is calculated and recorded as a total power before the mapping; The difference between the total power before the mapping and the total power after the mapping is calculated and recorded as the total power deviation.

5. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 4, characterized in that, The introduction of an energy conservation compensation source term at the boundary nodes of the grid mapping relationship topology according to the total power deviation to generate a corrected heat source distribution field comprises: Boundary heat grid cells on the conductor surface are identified in the grid mapping relationship topology and marked as boundary heat grid cells; The total number and total volume of the boundary heat grid cells are counted; The total power deviation is distributed according to the volume proportion of the boundary heat grid cells to obtain a power compensation amount of each boundary heat grid cell; The power compensation amount of each boundary heat grid cell is divided by the volume of the boundary heat grid cell to obtain a compensation heat source intensity of the boundary heat grid cell; The preliminary heat source intensity of each heat grid cell is added to the corresponding compensation heat source intensity to generate the corrected heat source distribution field.

6. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The construction of the multi-time scale coupling transfer function from electromagnetic loss to temperature rise comprises: For each annular sub-region in the layered heat source intensity distribution model, the square of the ratio of the characteristic size of the annular sub-region to the material thermal diffusivity coefficient is calculated and recorded as a thermal diffusion characteristic time constant of the annular sub-region; According to the thermal diffusion characteristic time constants of the annular sub-regions, a multi-order RC thermal network equivalent circuit is constructed, each RC node in the multi-order RC thermal network equivalent circuit corresponding to an annular sub-region; The multi-order RC thermal network equivalent circuit is subjected to Laplace transformation to obtain a frequency domain transfer function expression; A phase correction term is introduced into the frequency domain transfer function expression; The frequency domain transfer function expression is subjected to inverse Laplace transformation to obtain the multi-time scale coupling transfer function in the time domain.

7. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The application of the multi-time scale coupling transfer function to the corrected heat source distribution field to generate a transient heat field prediction result comprises: The corrected heat source distribution field is grouped according to the annular sub-regions to which the heat grid cells belong to obtain a heat source time sequence corresponding to each annular sub-region; The heat source time sequence corresponding to each annular sub-region is subjected to the multi-time scale coupling transfer function to obtain a transient heat field prediction result of the annular sub-region. applying a sub-transfer function corresponding to the annular sub-region in the multi-time scale coupling transfer function to a heat source time sequence of the annular sub-region respectively to obtain a temperature rise response time sequence of each annular sub-region; mapping the temperature rise response time sequence of each annular sub-region back to a thermal grid cell of the initial thermal simulation grid data according to the spatial position; performing time domain superposition on the temperature rise responses from different annular sub-regions received by each thermal grid cell to obtain the total temperature rise of the thermal grid cell; combining the environmental temperature boundary condition of the electronic equipment, converting the total temperature rise of each thermal grid cell into an absolute temperature to generate the transient thermal field prediction result.

8. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The phase deviation analysis of the characteristic frequency component of the temperature response in the transient thermal field prediction result and the excitation frequency of the electromagnetic field simulation data includes: selecting a temperature monitoring point on the surface of the conductor in the transient thermal field prediction result to extract the temperature time sequence of the monitoring point; performing short-time Fourier transform on the temperature time sequence to obtain the time-frequency spectrum distribution of the temperature response; identifying the frequency component with the maximum amplitude in the time-frequency spectrum distribution as the dominant characteristic frequency; extracting the excitation frequency and its multiple frequency components of the electromagnetic field simulation data, and performing frequency matching with the dominant characteristic frequency; for the frequency components that match successfully, calculating the phase difference between the temperature time sequence and the electromagnetic loss power time sequence at this frequency, denoted as the phase deviation; statistically analyzing the spatial distribution of the phase deviation in different annular sub-regions to identify abnormal regions with phase deviation exceeding a preset phase threshold.

9. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 1, characterized in that, The dynamic adjustment of the heat source delay phase of each annular sub-region in the hierarchical heat source intensity distribution model according to the phase deviation analysis result, and the updating of the transient thermal field prediction result include: for each annular sub-region, extracting the phase deviation corresponding to the annular sub-region; converting the phase deviation into a time delay amount, which is equal to the phase deviation divided by the angular frequency of the corresponding frequency; in the hierarchical heat source intensity distribution model, applying a time shift operation of the time delay amount to the heat source time sequence of the annular sub-region; re-substituting the heat source time sequence after the time shift operation into the multi-time scale coupling transfer function to recalculate the temperature rise response of the annular sub-region; after adjusting all annular sub-regions as described above, re-executing the temperature rise time domain superposition process to generate the updated transient thermal field prediction result.

10. The multi-physical field strong coupling oriented electronic device thermal field fast prediction method according to claim 4, characterized in that, The inverse distance weighted interpolation method is used to calculate the preliminary heat source intensity of the thermal grid cell based on the hierarchical heat source intensity distribution model data of the plurality of electromagnetic grid nodes, including: calculating the Euclidean distance from the centroid of the thermal grid cell to each electromagnetic grid node, denoted as the node distance; taking the inverse of each node distance and performing normalization processing to obtain the distance weight of each electromagnetic grid node; extracting the unit volume heat source intensity of the annular sub-region to which each electromagnetic grid node belongs; multiplying the unit volume heat source intensity of each electromagnetic grid node by the corresponding distance weight and performing weighted summation to obtain the preliminary heat source intensity; When the centroid of the thermal grid unit is less than a preset distance threshold from the node distance of a certain electromagnetic grid node, the unit volume heat source intensity of the annular sub-region to which the electromagnetic grid node belongs is directly assigned as the preliminary heat source intensity.

Images