A transient electromagnetic-thermal coupling numerical calculation method based on finite difference format

Abstract

The application discloses a transient electromagnetic-thermal coupling numerical calculation method based on a finite difference format, and belongs to the cross field of integrated circuit technology and computational physics. The calculation method of the application faces the actual demand of a system-level package, utilizes an alternating direction implicit finite difference format, solves state equations of electromagnetic fields and temperature fields in a time domain, and guarantees unconditional stability of an algorithm; a Debye model is added to materials to describe the variation characteristics of medium material parameters with frequency, and adaptability to high-frequency signals is guaranteed; meanwhile, related physical quantities of electromagnetic fields and temperature fields are set in a set of coordinate systems, so that the solving process meets grid consistency; a time factor is added to the time step of the temperature field, and the problem of inconsistency of time scales between two physical fields is effectively coordinated. The formed electromagnetic-thermal coupling solving method can accurately describe the transient distribution of electromagnetic fields and temperature fields of a target, so that the real physical process of the simulation target changing with time during work is restored.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of integrated circuit technology and computational physics, and in particular to a transient electromagnetic-thermal coupling numerical calculation method based on a finite difference scheme. Background Technology

[0002] System-in-package (SiP) is one of the advanced packaging technologies in the integrated circuit field and is considered an important means of continuing Moore's Law. This technology integrates multiple functional bare chips into a single package, offering higher integration density and design flexibility. However, the scale of integration leads to significant heat generation in these products, and the coupling phenomenon between electromagnetic and temperature fields becomes exceptionally prominent. This is especially true for critical electrical interconnect structures such as substrates, redistribution layers, and vertical through-silicon vias (TSVs), where high interconnect density and low thermal conductivity materials easily cause heat accumulation, leading to signal quality degradation, distortion, and even electrical performance failure. Therefore, analyzing these electromagnetic-thermal coupling phenomena is crucial for ensuring the normal operation of SiP products. Multiphysics numerical simulation (MPS) can accurately describe the physical processes of these phenomena, enabling early-stage design evaluation, analysis of potential risks, and timely feedback of feasible optimization suggestions, thus improving product development efficiency. Consequently, this type of MPS simulation analysis has received widespread attention in the industry and has become a research hotspot in recent years.

[0003] Based on the mechanism of electromagnetic-thermal coupling, it can be classified into different types: Joule heating caused by conduction current and dielectric loss; induced heating caused by induced current and magnetic loss; and microwave heating caused by resistance, dielectric, and magnetic loss at high frequencies. For system-in-package (SiP), considering its geometry, operating frequency, and material characteristics, Joule heating is the main source of coupling. This is a typical bidirectional coupling problem. Under the influence of an electric field, the ohmic loss of the conduction current and the dielectric loss of the dielectric material induce thermal effects, and the key parameters involved are conductivity σ and complex relative permittivity ε. r (ε r =ε r '+jε r ”, where ε r ' represents the real part, j represents the imaginary unit, and ε r The complex relative permittivity (μ) and thermal conductivity (κ) typically change nonlinearly with temperature, resulting in different Joule heating rates at different temperatures. Meanwhile, the complex relative permittivity exhibits a frequency-dependent function, which makes it difficult to analyze this type of electromagnetic-thermal coupling phenomenon.

[0004] On the other hand, the frequency of signals in electromagnetic fields is typically in the GHz range, with time periods on the nanosecond scale; while temperature changes relatively slowly, potentially taking seconds or minutes. This significant difference makes the computational cost of using transient formulas to solve time-harmonic electromagnetic problems extremely high. Currently, a common approach to address this time-scale inconsistency is to first calculate the periodic average loss of the electromagnetic signal and then use these losses as a constant heat source to solve the steady-state heat transfer problem. However, for critical electrical interconnect structures in system-in-package (SISP), time-varying electrical and thermal characteristics are a more effective means of analyzing design feasibility and reliability, which existing steady-state analysis methods cannot address.

[0005] The above analysis shows that transient electromagnetic-thermal coupling analysis of key electrical interconnect structures in system-in-package (SiP) can reconstruct the coupling process of electromagnetic and temperature fields and their changes over time. This is a crucial basis for optimizing electrical interconnects in the early stages of design. Simultaneously solving Maxwell's equations for electromagnetic fields and the heat conduction equations in the time domain accurately describes the interaction between the two physical fields. Coordinating the time differences between the two physical fields in the time domain improves the algorithm's executability and enhances its ability to analyze transient problems. Therefore, it is necessary to establish an efficient transient numerical analysis method to accurately analyze the electrical and thermal performance of key electrical interconnect structures in SiP, providing necessary theoretical support and efficient simulation tools for experiments. Summary of the Invention

[0006] The purpose of this invention is to provide a transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme to solve the problems in the background art.

[0007] To address the aforementioned technical problems, this invention provides a transient electromagnetic-thermal coupling numerical calculation method based on a finite difference scheme, comprising:

[0008] Step 1: Establish the discrete space of the simulation target, specify the solution domain for the electromagnetic field and temperature field, which completely encloses the calculation target; set the discrete space and time step, both of which simultaneously meet the accuracy requirements of the discrete space in the electromagnetic field and temperature field.

[0009] Step 2: Establish a three-dimensional geometric model of the simulation target and perform mesh generation;

[0010] Step 3: Initialize all physical field components in the computational domain;

[0011] Step 4: Set the electromagnetic signal excitation and background temperature;

[0012] Step 5: Solve the Maxwell equations in the entire computational space using the ADI scheme to obtain the electric field strength, magnetic field strength, and polarization intensity at the current moment.

[0013] Step 6: Solve for the equivalent heat source in the temperature field computation space;

[0014] Step 7: Solve the heat conduction equation in the temperature field computation space using the ADI scheme and update the material parameters;

[0015] Step 8: Repeat steps 5, 6, and 7 until the specified number of iterations is completed, and output the calculation results of each physical field.

[0016] In one embodiment, in step one, a discrete space is established in the Cartesian coordinate system with a right hexahedron as the basic unit, and the calculation region of the temperature field is smaller than that of the electromagnetic field. Within the calculation region, the discrete spatial step size in the x, y and z directions is set to Δx, Δy and Δz respectively, and the time step size is set to Δt. The above spatial and time steps are consistent in the electromagnetic field and the temperature field.

[0017] In one implementation, step two includes:

[0018] Establish a three-dimensional geometric model of the key electrical interconnect structure in the system-level packaging, and use coordinates to describe the position of each node;

[0019] Based on the node coordinates and the discrete space step size, the entire computational target is divided into meshes to achieve the discretization of the geometric model.

[0020] In one implementation, step three includes two aspects:

[0021] First, the relative positions of the electric field strength, magnetic field strength, and temperature in space are set in a formal manner. That is, the electric field strength E is set on the edge of the hexahedral grid, the magnetic field strength H is set on the face center of the hexahedral grid, and the temperature T is set on the centroid of the hexahedral grid. The other physical quantities to be solved are set in the corresponding positions according to their relationship with E, H, and T.

[0022] Secondly, the above physical quantities are initialized, that is, each parameter is assigned a corresponding value, and each physical field quantity is set to zero. Among them, the Debye model is used to describe the relative permittivity ε. r .

[0023] In one implementation, step four includes two aspects:

[0024] First, electromagnetic signal excitation is set up, and the forms of excitation include time-harmonic signals of specific frequencies, Gaussian pulse sources, differential Gaussian pulse sources, and modulated Gaussian pulses;

[0025] Second, set the initial background temperature, which can be a constant or a function that changes over time.

[0026] In one implementation, step five uses the ADI scheme for solving the problem and includes an analysis of the Debye model. Since the ADI scheme requires implicit computation to be performed alternately in different directions, a time step of n+1 / 2 is set as a temporary variable in the n-time step and the n+1-time step.

[0027] In one implementation, step six includes two aspects:

[0028] First, a time factor is set, which is the time step n of the electromagnetic field and the time step n of the temperature field. t Multiples between;

[0029] Secondly, when solving for the equivalent heat source, the equivalent resistivity simultaneously includes both the ohmic loss in the metallic material and the dielectric loss in the non-metallic material.

[0030] In one implementation, step seven includes two aspects:

[0031] First, the heat conduction equation is solved in the temperature field computation space using the ADI scheme;

[0032] Second, based on the temperature results of each node and the relationship between material parameters and temperature, update the material parameters at the current temperature, including dielectric constant, magnetic permeability, resistivity, and thermal conductivity.

[0033] The present invention provides a transient electromagnetic-thermal coupling numerical calculation method based on a finite difference scheme, which has the following advantages:

[0034] (1) The numerical calculation method established in this invention can solve the transient characteristics of the electromagnetic field and temperature field of the target, and the target includes both metallic materials and dispersive media. Based on the traditional method, this invention adds a Debye model to the material to describe the variation characteristics of the medium material parameters with frequency, ensuring adaptability to high-frequency signals; at the same time, in a coordinate system, the relevant physical quantities of the electromagnetic field and temperature field are set to ensure mesh consistency in the solution process; a time factor is added to the time step of the temperature field to effectively coordinate the problem of inconsistent time scales between the two physical fields; the resulting electromagnetic-thermal coupled solution method can accurately describe the transient distribution of the electromagnetic field and temperature field of the target, thereby restoring the real physical process of the simulated target changing with time during operation;

[0035] (2) The established numerical calculation method is unconditionally stable. This invention uses the ADI scheme to solve Maxwell's equations for electromagnetic fields and the heat conduction equations for computational thermodynamics, effectively avoiding the stability constraints of traditional finite difference schemes and ensuring the algorithm remains unconditionally stable. At the same time, it lays the foundation for using different time factors. This difference scheme can also increase the time step, improving the execution efficiency of the algorithm while ensuring the accuracy of mesh generation. Attached Figure Description

[0036] Figure 1 This is a schematic diagram of the specific calculation process of the present invention.

[0037] Figure 2 A schematic diagram of the discrete grid cells used in this invention. Detailed Implementation

[0038] The following detailed description, in conjunction with the accompanying drawings and specific embodiments, provides a further detailed explanation of the transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme proposed in this invention. The advantages and features of this invention will become clearer from the following description. It should be noted that the accompanying drawings are all in a very simplified form and use non-precise scales, and are only used to facilitate and clarify the illustration of the embodiments of this invention.

[0039] To efficiently and accurately describe transient electromagnetic-thermal coupling phenomena in system-in-package (SiP), this invention simultaneously solves the Maxwell equations for the electromagnetic field and the heat conduction equations for the temperature field using a finite difference scheme. It employs the Debye model to describe the dielectric constant of the medium and utilizes the ADI scheme to circumvent stability constraints. A time factor is used to balance the time difference between the two physical fields, significantly improving algorithm efficiency and enabling the analysis of complex models. This is particularly relevant to the electromagnetic-thermal coupling simulation problem of key electrical interconnect structures in SiP within integrated circuits. Figure 1 The following is a detailed description of the specific calculation process of the present invention. The process is intended to explain the present invention and is not intended to limit it.

[0040] Step 1: Establish the discrete space of the simulation target, specify the solution region (i.e., the computational region) of the electromagnetic field and temperature field, which completely encloses the computational target; set the discrete space and time step, both of which simultaneously meet the accuracy requirements of the discrete space in the electromagnetic field and temperature field.

[0041] A discrete space is established in the Cartesian coordinate system, with the right hexahedron as the basic unit. The spatial coordinate range of the electromagnetic field calculation region is set in the x-direction as X. EM _n→X EM _p, set in the y-direction to Y EM _n→Y EM _p, set to Z in the z-direction EM _n→Z EM _p. The spatial coordinate range of the temperature field calculation region is set to X in the x-direction. T _n→X T _p, set in the y-direction to Y T _n→Y T _p, set to Z in the z-direction T _n→ZT _p. The computational domain for the temperature field is smaller than that for the electromagnetic field, i.e.: X EM _n<X T _n、X EM _p>X T _p、Y EM _n<Y T _n、Y EM _p>Y T _p、Z EM _n<Z T _n、Z EM _p>Z T Within the computational domain, the discrete spatial step sizes in the x, y, and z directions are set to Δx, Δy, and Δz, respectively, and the time step size is set to Δt. These spatial and time step sizes remain consistent in both the electromagnetic and temperature fields.

[0042] Step 2: Establish a three-dimensional geometric model of the simulation target and perform mesh generation.

[0043] A three-dimensional geometric model of the key electrical interconnect structure in the system-level packaging is established, and the position of each node is described using coordinates. Based on the node coordinates and the discrete space step size, the entire computational target is meshed to achieve discretization of the geometric model.

[0044] Step 3: Initialize all physical field components in the computational domain.

[0045] Based on all the coordinates of the discrete space obtained in steps one and two, according to Figure 2 The physical quantities required for solving the electromagnetic field and temperature field are set in the form of [formula missing] and initialized. In this invention, (i, j, k) represents the spatial coordinate position; subscripts x, y, z are added to the physical quantities to indicate the spatial direction of the physical field quantity; n represents the unit time (i.e., time step) in the electromagnetic field, n represents the current time, n+1 represents the next time, and n+1 / 2 represents the transition time between the current time and the next time; n... t This represents a unit of time in the temperature field.

[0046] The physical quantities included in an electromagnetic field and their initialization methods are as follows:

[0047] The three directional components of the electric field intensity E are:

[0048] E x | i,j,k =0, E y | i,j,k =0, E z | i,j,k =0;

[0049] The three directional components of the electric flux density D are:

[0050] D x | i,j,k =0, D y | i,j,k =0, D z | i,j,k =0;

[0051] The three directional components of the polarization intensity P are:

[0052] P x | i,j,k =0, P y | i,j,k =0, P z | i,j,k =0;

[0053] The three directional components of the magnetic field strength H are:

[0054] H x | i,j,k =0, H y | i,j,k =0, H z | i,j,k =0;

[0055] The three directional components of the magnetic flux density B are:

[0056] B x | i,j,k =0, B y | i,j,k =0, B z | i,j,k =0;

[0057] The dielectric constant in vacuum is ε₀ = 8.854 × 10⁻⁶. -12 F / m; Permeability in vacuum μ0 = 4π × 10 -7 N / A 2 The physical quantities included in the temperature field and their initialization methods are: power density Po| i,j,k =0, temperature T| i,j,k =0; Different values ​​are assigned to the remaining parameters according to the properties of different materials, including: relative permittivity ε r Relative permeability μ r Electrical conductivity σ, heat capacity C p Thermal conductivity κ, temperature coefficient of magnetic permeability α μ Temperature coefficient of dielectric constant α ε Temperature coefficient of electrical conductivity α σ Temperature coefficient of thermal conductivity α μ Considering the frequency-dependent dielectric properties of critical electrical interconnects in system-in-package (SIPP), this invention employs the Debye model to describe the relative permittivity ε. rThe relationship between the electric field strength, current flux density, and polarization intensity can be expressed by the following formula:

[0058] D(ω)=ε0(ε r (ω)'+jε r (ω)”)E(ω)

[0059] =ε0ε r (ω)E(ω)=ε0[ε ∞ +χ(ω)]E(ω)=ε0ε ∞ E(ω)+P(ω)

[0060]

[0061] Where ω is the angular frequency, ε(ω)' is the real part of the complex permittivity, and ε(ω)” is the imaginary part of the complex permittivity. ∞ ε is the dielectric constant at infinite frequency, Δε is the difference between the static dielectric constant and the dielectric constant at infinite frequency, j is the imaginary unit, τ0 is the relaxation time, and χ(ω) is the polarizability.

[0062] Step 4: Set the electromagnetic signal excitation and background temperature.

[0063] At the electromagnetic port of the simulated target, the excitation electric field is set as follows: Where f is the excitation source wavefunction, which may take forms including, but is not limited to, a time-harmonic signal of a specific frequency, a Gaussian pulse source, a differential Gaussian pulse source, and a modulated Gaussian pulse. Simultaneously, an initial background temperature T0 is set, which may be a constant or a function that varies with time.

[0064] Step 5: Solve the Maxwell equations in the entire computational space using the ADI scheme to obtain the electric field strength, magnetic field strength, and polarization intensity at the current moment.

[0065] This invention employs the ADI scheme to circumvent the inherent stability constraints of FDTD (Finite-Difference Time-Domain) and includes an analysis of the Debye model. Since the ADI scheme requires alternating implicit computations in different directions, a time step of n+1 / 2 is set as a temporary variable in both the n and n+1 time steps. The solution process is as follows.

[0066] (1) Solve for the electric field strength at time n+1 / 2:

[0067] 1)

[0068] at this time, c d4 =ε0ε ∞ -fd4 c d5 =1-f d2 ,

[0069] 2)

[0070] at this time, c d4 =ε0ε ∞ -f d4 c d5 =1-f d2 ,

[0071] 3)

[0072] at this time, c d4 =ε0ε ∞ -f d4 c d5 =1-f d2 ,

[0073] in,

[0074] (2) Solve for the field strength at time n+1 / 2:

[0075] 1)

[0076] 2)

[0077] 3)

[0078] (3) Solve for the field value of polarization intensity at time n+1 / 2:

[0079] 1)

[0080] 2)

[0081] 3)

[0082] in,

[0083] (4) Solve for the electric field strength at time n+1:

[0084] 1)

[0085] at this time, c d4 =ε0ε ∞ -fd4 c d5 =1-f d2 ,

[0086] 2)

[0087] at this time, c d4 =ε0ε ∞ -f d4 c d5 =1-f d2 ,

[0088] 3)

[0089] at this time, c d4 =ε0ε ∞ -f d4 c d5 =1-f d2 ,

[0090] in,

[0091] (5) Solve for the field strength at time n+1:

[0092] 1)

[0093] 2)

[0094]

[0095] 3)

[0096] (6) Solve for the field value of polarization intensity at time n+1:

[0097] 1)

[0098] 2)

[0099] 3)

[0100] in,

[0101] Step 6: Solve for the equivalent heat source in the temperature field computation space.

[0102] The time factor is set to m, which means that the electromagnetic field time step cycle is repeated m times for one temperature field time step cycle. Therefore, n t=m·n. When the time conditions are met for solving the temperature field, the dissipated power is calculated by averaging the effective values ​​of the electric field intensity at the center of the grid cells within the previous temperature field time step, as follows:

[0103] The effective value of the electric field intensity within one electromagnetic field time step is

[0104]

[0105] Average power dissipation in the previous temperature field time step:

[0106]

[0107] Therefore, the power dissipation is:

[0108] Where, σ eff The equivalent resistivity, which includes ohmic losses in metallic materials and dielectric losses in non-metallic materials, is denoted as σ. eff =σ+ωε0ε″

[0109] The power consumption at each grid node is then calculated as follows:

[0110] Step 7: Solve the heat conduction equation in the temperature field computation space using the ADI scheme and update the material parameters.

[0111] The heat conduction equation is solved using the Analytic Difference (ADI) scheme to obtain the spatial temperature distribution at the current moment, and the material parameters at each node are updated based on the temperature. Since the ADI scheme requires implicit calculations to be performed alternately in different directions, time steps n+1 / 3 and n+2 / 3 are set as temporary variables in time steps n and n+1. The solution process is as follows.

[0112]

[0113] in,

[0114] 4) Update the resistivity and dielectric constant of the material as follows:

[0115]

[0116]

[0117]

[0118]

[0119] Step 8: Repeat steps 5, 6, and 7 until the specified number of iterations is completed, and output the calculation results of each physical field.

[0120] The above description is merely a description of preferred embodiments of the present invention and is not intended to limit the scope of the present invention in any way. Any changes or modifications made by those skilled in the art based on the above disclosure shall fall within the protection scope of the claims.

Claims

1. A numerical calculation method for transient electromagnetic-thermal coupling based on a finite difference scheme, characterized in that, include: Step 1: Establish the discrete space of the simulation target, specify the solution domain for the electromagnetic field and temperature field, which completely encloses the calculation target; set the discrete space and time step, both of which simultaneously meet the accuracy requirements of the discrete space in the electromagnetic field and temperature field. Step 2: Establish a three-dimensional geometric model of the simulation target and perform mesh generation; Step 3: Initialize all physical field components in the computational domain; Step 4: Set the electromagnetic signal excitation and background temperature; Step 5: Solve the Maxwell equations in the entire computational space using the ADI scheme to obtain the electric field strength, magnetic field strength, and polarization intensity at the current moment. Step 6: Solve for the equivalent heat source in the temperature field computation space; Step six includes two aspects: First, a time factor is set, which is the time step of the electromagnetic field. n With temperature field time step n t Multiples between; Secondly, when solving for the equivalent heat source, the equivalent resistivity simultaneously includes both the ohmic loss in the metallic material and the dielectric loss in the non-metallic material. Step 7: Solve the heat conduction equation in the temperature field computation space using the ADI scheme and update the material parameters; The heat conduction equation is solved using the ADI scheme to obtain the spatial temperature distribution at the current moment, and the material parameters of each node are updated according to the temperature. Since the ADI scheme requires implicit calculations to be performed alternately in different directions, time steps n+1 / 3 and n+2 / 3 are set as temporary variables in time steps n and n+1. The solution process is as follows: in, ; (4) Update the resistivity and dielectric constant of the material as follows: Step 8: Repeat steps 5, 6, and 7 until the specified number of iterations is completed, and output the calculation results of each physical field.

2. The transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme as described in claim 1, characterized in that, In step one, a discrete space with a right hexahedron as the basic unit is established in the Cartesian coordinate system, and the calculation area of ​​the temperature field is smaller than that of the electromagnetic field. Within the computational domain, the discrete spatial step sizes in the x, y, and z directions are set to Δx, Δy, and Δz, respectively, and the time step size is set to Δt. These spatial and time step sizes remain consistent in both the electromagnetic and temperature fields.

3. The transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme as described in claim 2, characterized in that, Step two includes: Establish a three-dimensional geometric model of the key electrical interconnect structure in the system-level packaging, and use coordinates to describe the position of each node; Based on the node coordinates and the discrete space step size, the entire computational target is divided into meshes to achieve the discretization of the geometric model.

4. The transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme as described in claim 3, characterized in that, Step three includes two aspects: First, the relative positions of electric field strength, magnetic field strength, and temperature in space are defined by the form of the electric field strength. E The magnetic field strength is set on the edges of the rectangular hexahedral grid. H Set at the face center and temperature of the hexahedral mesh T The centroid of the hexahedral mesh is set, and the remaining physical quantities to be solved are based on... E , H and T The relationships are set in the corresponding positions; Secondly, the above physical quantities are initialized, that is, each parameter is assigned a corresponding value, and each physical field quantity is set to zero. Among them, the Debye model is used to describe the relative permittivity ε. r .

5. The transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme as described in claim 4, characterized in that, Step four includes two aspects: First, electromagnetic signal excitation is set up, and the forms of excitation include time-harmonic signals of specific frequencies, Gaussian pulse sources, differential Gaussian pulse sources, and modulated Gaussian pulses; Second, set the initial background temperature, which can be a constant or a function that changes over time.

6. The transient electromagnetic-thermal coupling numerical calculation method based on the finite difference scheme as described in claim 5, characterized in that, Step seven includes two aspects: First, the heat conduction equation is solved in the temperature field computation space using the ADI scheme; Second, based on the temperature results of each node and the relationship between material parameters and temperature, update the material parameters at the current temperature, including dielectric constant, magnetic permeability, resistivity, and thermal conductivity.

Images