A frequency domain electromagnetic calculation method

Abstract

The present invention relates to the field of electromagnetic computing, and specifically to a frequency domain electromagnetic computing method, which can realize efficient electromagnetic computing of complex integrated circuit structures. The technical scheme includes: dividing the calculated area into two parts, an outer area and an inner area, according to the material distribution characteristics of the integrated circuit structure; using a two-dimensional finite element method to analyze and solve the outer area, and when solving, the electromagnetic influence of the inner area on it is equivalent to an equivalent current term containing a surface admittance operator, which is added to the electric field partial differential equation obtained by the non-homogeneous Helmholtz equation in the outer area to obtain the outer area electric field equation, and the outer area field distribution is obtained according to the outer area electric field equation; when solving the inner area, according to the outer area field distribution, the field distribution of the inner area boundary is obtained by the condition of longitudinal electric field continuity, and finally the inner area field distribution is obtained by combining the electric field integral equation. The present invention is suitable for frequency domain electromagnetic efficient computing of complex integrated circuit structures.

Description

Technical Field

[0001] The present invention relates to the field of frequency domain electromagnetic calculation, and in particular to a frequency domain electromagnetic calculation method. Background Art

[0002] With the development of packaging-level system integration technology, the device density on the chip is getting higher and higher, and the circuit structure is becoming more and more complex. The electromagnetic interference caused is inevitable. How to achieve good electromagnetic modeling to optimize the packaging design is becoming more and more important. The core of electromagnetic modeling is the electromagnetic numerical calculation method. At present, the traditional electromagnetic numerical calculation method has low accuracy and low calculation efficiency due to its universality and generality requirements, and is not very suitable for the calculation of complex integrated circuit structures.

[0003] Therefore, there is a need for an efficient frequency-domain electromagnetic computation method suitable for complex integrated circuit structures. Summary of the invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a frequency domain electromagnetic calculation method that can realize efficient electromagnetic calculation of complex integrated circuit structures and improve calculation efficiency.

[0005] The present invention adopts the following technical scheme to achieve the above-mentioned purpose. The present invention provides a frequency domain electromagnetic calculation method, which includes:

[0006] S1. Dividing the calculated area into an outer area and an inner area according to the material distribution characteristics of the integrated circuit structure;

[0007] S2. Use the two-dimensional finite element method to analyze and solve the outer region. When solving, the electromagnetic influence of the inner region on it is equivalent to the equivalent current term containing the surface admittance operator, which is added to the electric field partial differential equation obtained by the non-homogeneous Helmholtz equation in the outer region to obtain the electric field equation of the outer region. According to the electric field equation of the outer region, the field distribution of the outer region is obtained;

[0008] S3. When solving the inner region, according to the field distribution in the outer region, the field distribution at the boundary of the inner region is obtained by the condition of longitudinal electric field continuity, and finally the field distribution in the inner region is obtained by combining the electric field integral equation.

[0009] Furthermore, S1 specifically includes dividing the non-uniform medium, background medium and anisotropic medium in the integrated circuit structure into an outer region, and dividing the high conductivity medium and other structures into an inner region.

[0010] Furthermore, S2 specifically includes solving the field distribution in the outer region. First, the equation for the outer region and the boundary satisfying the standard form of two-dimensional finite element analysis is obtained by the non-homogeneous Helmholtz equation. The equation is as follows:

[0011]

[0012] in φ represents the electric field in the outer region is the external current source that may exist in the outer region and the boundary, k0 represents the free space wave number, k2 represents the outer region wave number, Z2 represents the outer region wave impedance, ε r2 , μ r2 represents the relative permittivity and relative permeability of the medium in the outer region;

[0013] The equivalent current term J containing the surface admittance operator is added to the equivalent electromagnetic influence of the inner region on the outer region eq After correction, the final equation of two-dimensional finite element analysis satisfied by the outer region and the boundary is obtained. The final equation is as follows:

[0014]

[0015] in is the longitudinal electric field on the boundary of the inner region, and κ is the differential surface admittance operator.

[0016] Furthermore, S3 specifically includes: solving the final equation of the two-dimensional finite element analysis to obtain the longitudinal electric field of the outer region including its boundary According to the boundary condition of longitudinal electric field continuity, the longitudinal electric field at the inner region boundary is obtained According to the coupling relationship between electric field and magnetic field, The tangential magnetic field on the inner region boundary is calculated Finally, the electric field distribution in the inner area is calculated according to the electric field integral equation The calculation is as follows:

[0017] Where T = 1, represents the Green's function, represents the imaginary unit, represents the second-kind zero-order Hankel function, n′ represents the external normal vector at the source point, Represents an arbitrary vector at the defined source point.

[0018] The beneficial effects of the present invention are:

[0019] The present invention combines the advantages of the finite element method and the method of moments in processing different circuit structures, realizes efficient electromagnetic modeling of complex integrated circuit structures, and provides a basis for the design optimization of packaged integrated circuits.

[0020] (1) The finite element method can make up for the shortcomings of the moment method in dealing with inhomogeneous and anisotropic media, and give full play to its advantages in dealing with complex media structures;

[0021] (2) The method of moments is applied to the interior of the calculation area. Compared with the finite element method, it has higher computational efficiency and requires less computing resources when dealing with fine structures or high conductivity media.

[0022] (3) By selecting appropriate basis functions in different regions, boundary-free coupling is achieved at the interface between regions. Without the need to artificially introduce other additional boundary conditions, the physical model can be equivalent to obtain accurate near-field and far-field results, thereby improving the computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 is a flow chart of a frequency domain electromagnetic calculation method provided by an embodiment of the present invention;

[0024] Figure 2 It is a schematic diagram of the decomposition of the inner and outer regions provided by an embodiment of the present invention;

[0025] Figure 3 is a schematic diagram of an equivalent electromagnetic structure of an outer region provided by an embodiment of the present invention;

[0026] Figure 4 It is a flow chart of a frequency domain electromagnetic calculation method using conformal grid division provided by an embodiment of the present invention;

[0027] Figure 5 It is a flow chart of a frequency domain electromagnetic calculation method using non-conformal grid division provided by an embodiment of the present invention;

[0028] Figure 6 It is a schematic diagram of a 2×3 cable array placed on a layered medium provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0029] To make the purpose, technical solution and advantages of the embodiments of the present invention more clear, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0030] The present invention proposes an accurate and efficient frequency domain electromagnetic calculation method based on the finite element method-moment method for complex electromagnetic structures containing inhomogeneous media, high conductivity media, etc. When conformal meshing is used for implementation, by selecting appropriate basis functions in different regions, boundary-free coupling is achieved at the interface between regions. Without the need to artificially introduce additional boundary conditions, the physical model equivalence can be achieved and accurate near-field and far-field results can be obtained. In order to overcome the restriction that the discrete results on the boundary must coincide with each other and the discrete forms must be exactly the same when conformal meshing is implemented, the present invention proposes a specific implementation method using non-conformal meshing. The field value transfer and coupling between the non-conformal grids of the two regions are realized by the transfer matrix, and the calculation efficiency is further improved while ensuring the accuracy.

[0031] The specific plan is as follows:

[0032] 1. Frequency domain electromagnetic calculation method based on finite element method-moment method

[0033] like Figure 1 As shown, the overall implementation idea of ​​the frequency-domain efficient new electromagnetic calculation method based on the finite element method-moment method provided by the present invention is: first, the calculated area is divided into two parts, the outer area and the inner area, according to the material distribution characteristics of the complex integrated circuit structure. Then the two-dimensional finite element method is used to analyze and solve the outer area, and the two-dimensional moment method is used to analyze and solve the inner area. The solution process of the inner and outer areas is not completely independent, but a joint solution process. When solving the outer area, it is necessary to add the electromagnetic influence of the inner area on it to the equivalent current term containing the surface admittance operator to the electric field partial differential equation obtained by the non-homogeneous Helmholtz equation in the outer area, so as to obtain the electric field equation that it finally satisfies; when solving the inner area, it is necessary to first calculate the field distribution of the outer area from the electric field equation that is finally satisfied by the outer area, and then obtain the field distribution of the boundary of the inner area from the condition of the continuity of the longitudinal electric field, and finally obtain the field distribution of the inner area by combining the electric field integral equation.

[0034] Specifically, inhomogeneous media, background media, anisotropic media, etc. in complex integrated circuit structures are generally classified as outer regions, while high conductivity media and other complex structures are generally classified as inner regions.

[0035] by Figure 2 Taking the structure shown in the figure as an example, the boundary between the inner region and the outer region is denoted as Γ, and n in the figure represents the unit normal vector of each point on the contour. The longitudinal electric field in the inner region is denoted as The longitudinal electric field on the inner region boundary, i.e., the contour Γ, is The longitudinal electric field in the outer region is The longitudinal electric field on the outer region boundary, i.e., the contour Γ, is

[0036] To solve the field distribution in the outer region, we first use the non-homogeneous Helmholtz equation to obtain the equation for the outer region and boundary that satisfies the standard form of two-dimensional finite element analysis:

[0037]

[0038] in φ represents the electric field in the outer region is the external current source that may exist in the outer region and the boundary, k0 represents the free space wave number, k2 represents the outer region wave number, Z2 represents the outer region wave impedance, ε r2 , μ r2Represents the relative dielectric constant and relative magnetic permeability of the outer region medium. This equation does not take into account the influence of the inner region on the electromagnetic field distribution of the outer region, so it is not the electric field equation actually satisfied by the outer region. It is necessary to add the equivalent current term J containing the surface admittance operator equivalent to the electromagnetic influence of the inner region on the outer region. eq After correction, the final equation of the two-dimensional finite element analysis satisfied by the outer region and the boundary is as follows:

[0039]

[0040] in is the longitudinal electric field on the boundary of the inner region, i.e., the contour Γ. κ is the differential surface admittance operator, which is also the key to the coupling of the electromagnetic influence of the inner region to the external region. The derivation process of the differential surface admittance operator κ is introduced below.

[0041] Figure 2 The Helmholtz equation satisfied by the electric field in the inner region is: Where k1 represents the wave number in the inner region. The boundary condition of the longitudinal electric field continuity is met at the boundary. According to Maxwell's equations, the expression of the boundary tangential magnetic field can be calculated from the boundary longitudinal electric field:

[0042]

[0043] in, represents the normal derivative of the electric field at the inner region boundary, represents the tangential magnetic field on the boundary of the inner region, and μ1 represents the magnetic permeability of the inner region.

[0044] In order to derive the surface admittance operator κ, the present invention requires Figure 2 The structure shown in the figure is equivalent to the structure shown in the figure. The specific steps are as follows: replace the inner region medium with the outer region medium to obtain the equivalent electromagnetic structure for the outer region; secondly, introduce an equivalent current or an equivalent magnetic current to ensure that the relevant electromagnetic field remains unchanged before and after the equivalent. Here, it is necessary to ensure that the electromagnetic field outside the contour Γ is the same as before the equivalent. Due to the boundary condition that the longitudinal electric field remains unchanged before and after the equivalent, it is only necessary to introduce an equivalent current. No magnetic current is required. Figure 2 The equivalent electromagnetic structure outside the structure shown is as follows Figure 3 shown.

[0045] The electric field in the equivalent electromagnetic structure is Can get The Helmholtz equation satisfied is: Where k2 is the wave number in the outer region of the equivalent electromagnetic structure. The boundary condition of longitudinal electric field continuity is also satisfied at the boundary, and the tangential magnetic field on the boundary of the equivalent structure is obtained again through Maxwell equations:

[0046]

[0047] The equivalent current can be obtained from the equivalent front and rear boundary tangential magnetic fields: Will Substitution get:

[0048]

[0049] The above formula reflects the relationship between the longitudinal electric field and the equivalent current on the boundary, from which the surface admittance operator κ can be defined. The equivalent current obtained by the boundary longitudinal electric field under the action of the operator κ is expressed as:

[0050] Through the contour integral method, the contour integral equation is established, and finally the detailed expression of the differential surface admittance operator is obtained. From the above derivation process, it can be seen that, unlike the existing calculus equation hybrid method, the method of the present invention realizes coupling between the inner region and the outer region through the surface equivalent current, and does not need to artificially introduce other additional boundary conditions, which can be completely equivalent to the original physical model.

[0051] By solving the final equation of the outer region by substituting the surface equivalent current term, the electric field distribution in the outer region can be accurately calculated. On this basis, other required parameters can be further calculated, such as parasitic electromagnetic parameters, radar scattering cross section, etc.

[0052] Regarding the solution of the field distribution in the inner region, it can be seen from the condition of the continuity of the longitudinal electric field that in the electric field obtained by solving the final equation of the outer region, not only the electric field in the outer region is accurate, but the electric field on the boundary of the inner region is also real and reliable. Therefore, the internal electric field can be solved based on the electric field on the boundary of the inner region.

[0053] First, the final equation of the two-dimensional finite element analysis is solved to obtain the longitudinal electric field in the outer region including its boundary. According to the boundary condition of longitudinal electric field continuity, the longitudinal electric field at the inner region boundary is obtained According to the coupling relationship between electric field and magnetic field, The tangential magnetic field on the boundary of the inner region is calculated Finally, the electric field distribution in the inner area is calculated according to the electric field integral equation The calculation is as follows:

[0054] Where T = 1, represents the Green's function, represents the imaginary unit, represents the second-kind zero-order Hankel function, n′ represents the external normal vector at the source point, Represents an arbitrary vector at the defined source point.

[0055] 2. Implementation of the finite element method-moment method using conformal meshing

[0056] The specific steps of the frequency domain electromagnetic calculation method using conformal mesh generation are as follows: Figure 4 As shown in the figure: (1) Divide the region and mesh it to obtain a conformal network; (2) Establish the functional of the final equation of the outer region; (3) Substitute the discrete results of the outer region and its boundary into the functional of the sub-unit; (4) Accumulate the functionals of the sub-regions to obtain the total functional of the outer region; (5) Use the Galerkin method to obtain the differential surface admittance operator and substitute it into the total functional; (6) Let the total functional be differentiated to zero and obtain the matrix equation; (7) Solve the matrix equation to obtain the external field distribution; (8) The field distribution at the boundary of the inner region is obtained from the longitudinal electric field; (9) The field distribution of the inner region is obtained from the electric field integral equation.

[0057] Specifically, the calculation area can be divided into an inner area and an outer area according to its material properties and structural properties. The same grid is used to divide the inner and outer areas, so that the discrete results on the boundaries of the two areas overlap with each other, and the discrete forms are exactly the same, that is, the boundary discrete nodes completely overlap. In this way, the coefficients corresponding to the nodes on the boundaries of the inner and outer areas can be completely shared, so that the field values ​​at the boundaries of the two areas can be directly transferred, which facilitates the processing of electromagnetic influences between areas.

[0058] The computer solution of the field distribution in the outer region uses the two-dimensional finite element method. First, the Rayleigh-Ritz method is used to establish the functional of the final equation of the outer region. The functional of the final equation is as follows:

[0059]

[0060] Discretize the outer region to obtain M triangular units, is the i-th two-dimensional linear basis function defined on the triangular element e obtained by discretization in the outer region, When the field value of any point on the triangle element e is represented by the two-dimensional linear basis function mentioned above The coefficient of e is the electric field at any point on the triangular element e:

[0061]

[0062] Discretize its boundary to get M s Segment units, where is the i-th one-dimensional linear basis function defined on the segment unit s obtained by discretizing the boundary of the outer region, When the field value at any point on the segment unit s is represented by the aforementioned one-dimensional linear basis function The coefficient ofs is the electric field at any point on the segment unit s, When the equivalent current at any point on the segment unit s is represented by the one-dimensional linear basis function mentioned above The coefficient of is the equivalent current at any point on the segment unit s:

[0063]

[0064] φ e ,φ s as well as Substituting the functional of the final equation, we get the functional expression of each triangular element e in the outer region:

[0065]

[0066] The expressions of each matrix element are as follows:

[0067]

[0068] By summing up the functionals of each triangular unit, the total functional expression of the outer area is:

[0069]

[0070] Among them, [B~] is the N×N matrix expanded on the basis of accumulation, and N is the total number of nodes in the external area and the boundary;

[0071] Substitute the differential surface admittance operator matrix The total functional expression of the outer region is:

[0072]

[0073] in, By M s ×M s The N×N matrix is ​​obtained by expanding the differential surface admittance operator [κ], {φ1} is the boundary electric field vector of the inner region. Since conformal meshing is used and the boundary conditions of longitudinal electric field continuity are met, the boundary electric field vector of the inner region is equal to the boundary electric field vector of the outer region, which can be expanded to the electric field vector of the outer region {φ}. The total functional expression of the outer region is:

[0074]

[0075] in,

[0076] The solution method of [κ] is: use piecewise triangular basis functions to discretize the longitudinal electric field, tangential magnetic field, equivalent current, etc. on the boundary of the inner region, then substitute the discretization results into the contour integral equation, and then perform the Galerkin test to establish a matrix equation for the electric field and magnetic field coefficients, so that the boundary longitudinal electric field is used to express the solution of the boundary tangential magnetic field. Then, by the same method, the solution of the equivalent tangential magnetic field of the boundary is obtained, and [κ] can be obtained by the definition of the equivalent current. The present invention uses piecewise triangular basis functions when discretizing the boundary of the inner region to ensure that the discretization results of the outer region and the inner region are the same on the interface. Specifically, each piecewise triangular basis function on the boundary of the inner region is composed of a one-dimensional linear basis function defined on two adjacent segments on the boundary of the outer region. Different from the one-dimensional linear basis function on the boundary of the outer region, the piecewise triangular basis function on the boundary of the inner region is defined on the node, not on the segment. Based on this, the outer region and the inner region can obtain the same discrete electric field and discrete current on the boundary, thereby realizing boundary-free coupling.

[0077] In the total functional expression of the outer region, the expression after substituting each matrix element into the expansion function is:

[0078]

[0079] b i e 、c i e are the known coefficients before x and y corresponding to the ith basis function on the triangular unit e obtained by discretization in the outer region, Δ e is the area of ​​the triangular unit e obtained by discretizing the outer region. β e =-k0 2 ε r2 e . s represents the length of the sth segment, f e Indicates the possible external sources on the triangular element e obtained by discretizing the outer region. ij is defined as: when i=j, its value is 1, otherwise its value is 0.

[0080] Depend on We can get: {[K]+2[B]}{φ}=[b], let [A]=[K]+2[B], the above formula can be transformed into: [A]{φ}=[b]. The computer code solves the above non-homogeneous matrix equation to get the electric field distribution of the outer region including the boundary.

[0081] For the computer numerical solution of the inner region field distribution, firstly, according to the boundary condition of longitudinal electric field continuity, the inner region boundary electric field is obtained from the outer region boundary electric field obtained in the previous text, and the inner region boundary magnetic field distribution can be further obtained according to the coupling relationship between the electric field and the magnetic field.

[0082] Electric field distribution in the inner area To solve the problem, we use the piecewise triangular basis function to discretize the inner region boundary, and we get:

[0083]

[0084] In the formula Represents the position vector of each discrete point to be found in the inner region, represents the position vector of each discrete point on the inner region boundary, m represents the number of discrete segments of the inner region boundary, The nth piecewise triangular basis function representing the inner region boundary, h n 、e n The distribution represents the magnetic and electric fields at the first endpoint of the nth discrete segment on the boundary of the inner region, both of which are known quantities.

[0085] From the above formula, we can see that the electric field of each discrete point in the inner region can be expressed as the weighted sum of the known electric and magnetic fields at the discrete points on the boundary of the inner region. Therefore, by calculating the loop setting, The electric fields of all the points in the inner area can be obtained by going through the positions of all the points in the inner area one by one. According to the coupling relationship between the electric field and the magnetic field, all the magnetic fields and other electromagnetic parameters of interest of the points in the inner area can be further obtained.

[0086] 3. Specific implementation of the finite element method-moment method using non-conformal meshing

[0087] When using conformal meshing to implement the proposed finite element method-moment method frequency domain electromagnetic calculation method, it is impossible to flexibly adjust the regional meshing size according to the different properties of the materials and structures in the two regions. While losing the flexibility of the algorithm, it also sacrifices a certain degree of calculation efficiency. In order to overcome this limitation, the specific implementation of using non-conformal meshing will be introduced below. Figure 5 As shown in FIG, the steps of the method are: (1) dividing the region and meshing to obtain a non-conformal network; (2) determining the outer region unit where the inner region boundary node is located; (3) establishing a mapping relationship between the local and global numbers of the outer region unit nodes; (4) constructing a transfer matrix of the inner and outer boundary nodes; (5) using the finite element method to solve the outer region and bringing in the transfer matrix; (6) using the moment method to solve the inner region.

[0088] The calculation area can be divided into inner and outer areas according to its material properties and structural properties, and a more appropriate mesh size is selected to obtain a non-conformal mesh. Specifically, since the inner area may contain high conductivity media, fine structures, etc., a relatively fine mesh is used, while the outer area only contains background media such as air, so a sparser mesh is used.

[0089] Because of the introduction of non-conformal grids, the boundary discrete nodes no longer overlap, so the coefficients corresponding to the nodes cannot be shared. In order to realize the numerical transfer of different grids between regions, the boundary nodes of the inner region can be transformed to the outer region by interpolation. Assume that there are m discrete nodes on the boundary of the inner region, and the m nodes transferred to the outer region by interpolation s This interpolation transformation can be achieved by introducing a transfer matrix to describe.

[0090] In the formula, and Respectively represent the column vectors composed of the longitudinal electric field expansion coefficients before and after interpolation, that is:

[0091]

[0092] Among them, e i (i=1,2,...,m) and Respectively and The expansion coefficient of .

[0093] When non-conformal meshing is used to implement the proposed finite element method-moment method frequency domain calculation method, the electric field, current and other physical quantities are expanded using linear basis functions, which are essentially linear interpolation functions, so they can be used to construct the transfer matrix The linear interpolation function expression defined in the triangular mesh unit e is as follows:

[0094]

[0095] In the formula, a j e 、b j e 、c j e are the constant term corresponding to the jth basis function of the triangular mesh unit e and the known coefficients before x and y, Δ e Represents the area of ​​unit e. Δ j (j=1,2,3) is defined inside the triangular element e.

[0096] To construct the transfer matrix There are three basic steps required:

[0097] (1) Determine the triangular unit where the boundary node of the inner region is located in the outer region by defining the relationship between the linear interpolation functions on the unit e. The implementation method is as follows:

[0098] Among them, x0 and y0 represent the x and y coordinates of the boundary nodes of the inner region respectively; if the triangle unit satisfies this implementation method, the unit where the boundary node is located can be determined, and To represent the values ​​of the three interpolation coefficients;

[0099] (2) When the outer region triangular unit corresponding to each inner region boundary node is determined, the global number of the nodes constituting the triangular unit in the outer region is obtained. s Then establish a local number k (k = 1, ..., m s ) to its global number M(k), where M(·) represents the corresponding mapping function;

[0100] (3) Based on this mapping relationship, a transfer matrix is ​​constructed. For the i-th (i=1, 2, ..., m) node on the boundary of the inner region, it is located in the e-th node in the outer region. i triangle units, and the corresponding triangle vertices are globally numbered and The corresponding transfer matrix The element expression in is:

[0101] In this invention, since the inner region is meshed relatively finely, the nodes on its boundary will not all fall inside the triangle units in the outer region, but may fall on the edges of the triangle units. The only change in the derivation process is that the interpolation will no longer be based on the triangle unit interpolation, but on the edge unit interpolation. Therefore, the linear interpolation function will be reduced from two dimensions to one dimension, which contains two interpolation coefficients, denoted as The calculation expression is as follows:

[0102]

[0103] In the formula, represents the position vector of the two end points of the sth segment, l s Indicates the length of the segment. Respectively represent the distance from the interpolation point to the two end points.

[0104] In the transfer matrix After that, the finite element method can be used to solve the field distribution in the outer region. The specific steps are similar to those in conformal meshing, which include the following sub-steps: (1) Establish the functional of the final equation of the outer region; (2) Substitute the discretization results of the outer region and its boundary into the functional of the sub-unit; (3) Accumulate the functionals of the sub-regions to obtain the total functional of the outer region; (4) Use the Galerkin method to obtain the differential surface admittance operator and substitute it into the total functional; (5) Differentiate the total functional to zero and obtain the matrix equation; (6) Solve the matrix equation to obtain the outer field distribution.

[0105] Different from conformal meshing, in this case, after the differential surface admittance operator matrix {J}=[κ]{φ1} is introduced into the total functional in step (4), the total functional expression of the outer region is:

[0106]

[0107] However, due to the use of non-conformal meshing, the boundary discrete nodes no longer overlap, so the coefficients corresponding to the nodes cannot be shared. Therefore, the inner region boundary electric field vector {φ1} is not equal to the outer region boundary electric field vector. It is necessary to use the transfer matrix obtained above to represent the inner region boundary electric field vector with the outer region boundary electric field vector. Specifically, the inner region boundary electric field vector {φ1} is equal to The product of the electric field vector of the outer region boundary is then expanded to the electric field vector of the outer region {φ}. The total functional expression of the outer region is:

[0108]

[0109] in,

[0110] The method of moments is used to solve the field distribution in the inner region. The specific steps are similar to those in conformal meshing, including the following two sub-steps: (1) Using the transfer matrix to calculate the electric field distribution of the outer region boundary The electric field distribution at the inner region boundary is obtained (2) The field distribution in the inner region is obtained from the electric field integral equation.

[0111] Different from conformal decomposition, the electric field distribution at the boundary of the inner region is obtained from the electric field distribution at the boundary of the outer region not by using the boundary condition of the longitudinal electric field continuity, but by directly converting it using the transfer matrix obtained previously.

[0112] In order to improve the efficiency of electromagnetic modeling of complex integrated circuit structures, the present invention not only inherits the advantage of the moment method that only needs to discretize the conductor surface to reduce computing resources, but also combines the advantage of the finite element method that can handle non-uniform and anisotropic media well. Therefore, the calculation method based on the finite element method-moment method of the present invention has high efficiency and flexibility while ensuring the calculation accuracy, and has good performance in the frequency domain electromagnetic modeling calculation of complex integrated circuit structures.

[0113] The skin current of the cable array placed on the layered uniform medium is solved, and the corresponding structural parameters are as follows: Figure 6As shown in the figure, the entire cable array contains 6 cable structures, each of which contains three copper conductor columns with a circular cross-section and a radius of 0.25 mm, which are buried in a uniform dielectric substrate with a radius of 1 mm and a relative dielectric constant of 2.3. The entire cable array is placed 0.5 mm above the layered dielectric, and the relative dielectric constants of the layered dielectric from top to bottom are 2, 2.5, and 3. A TM polarized uniform plane wave with a frequency of 30 MHz is incident along the y-axis.

[0114] Table 1 Comparison of cable array FEM-Moment calculation methods and FEM calculation resources

[0115]

[0116] Table 1 shows the usage of computing resources when the finite element method and the finite element method-moment method calculation method of the present invention obtain the skin current density with the same accuracy. The codes are all implemented using Matlab, the computer CPU is Inteli7-7700 3.6GHz, and the memory is 64G. In addition, in order to ensure the fairness of the comparison, only single-threaded operation is used without any parallel processing. In the table, the fourth column ratio is defined as the ratio of the relevant quantity of the finite element method-moment method calculation method to the corresponding quantity of the finite element method. It can be seen from the second row in Table 1 that the peak skin current density obtained by the finite element method-moment method calculation method is 123964.9A / m 2 The result obtained by the finite element method is 123985.3A / m 2 , the ratio is 0.9997, so to a certain extent, it can be considered that the two methods obtain the skin current density with the same accuracy. On this basis, the computational resource consumption of the two methods is compared.

[0117] The third row in Table 1 shows the unknowns used by the two methods when solving the skin current density of the conductor inside the cable. The number of unknowns in the finite element method is 447485, while the finite element method-moment method calculation method only uses 159364, which is 64% less than the finite element method. In the actual solution process, for the copper conductor column inside the cable, the finite element method's segmentation size is 0.005mm, while the finite element method-moment method calculation method's segmentation size is only 0.02mm. This is because the finite element method-moment method calculation method combines the advantages of the moment method in dealing with conductors. It only needs to discretize the conductor surface, and does not need to use a very fine grid segmentation inside the conductor, thereby greatly reducing the number of unknowns required.

[0118] In terms of memory consumption, the calculation method of the finite element method-moment method has also been significantly reduced compared to the finite element method, as shown in the fourth row of Table 1, accounting for only 41%. From the fifth to the eighth row of Table 1, the calculation method of the finite element method-moment method and the calculation time of the finite element method for solving the skin current density are given in detail, including the total time, DSAO construction time, matrix filling time and solution time. As shown in the fifth row of Table 1, the total time of the finite element method is 10265.9s, and the total time of the finite element method-moment method is 976.5s, accounting for only 9%, which is a significant improvement. Specifically, although the calculation method of the finite element method-moment method requires the construction of a differential surface admittance operator, due to the significant reduction in unknown quantities, the time required for matrix filling and solving is greatly reduced, accounting for only 7% and 25% of the finite element method.

[0119] In summary, when obtaining the skin current density with the same accuracy, the calculation method based on the finite element method-moment method is more efficient than the traditional finite element method, with significant improvements in the number of unknowns, memory consumption and calculation time, which will become more significant as the inner area increases.

[0120] The above is only a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the form disclosed herein, and should not be regarded as excluding other embodiments, but can be used in various other combinations, modifications and environments, and can be modified within the scope of the concept described herein through the above teachings or the technology or knowledge of the relevant field. The changes and modifications made by those skilled in the art shall not deviate from the spirit and scope of the present invention, and shall be within the scope of protection of the claims attached to the present invention.

Claims

1. A frequency domain electromagnetic calculation method, characterized in that: The method comprises: S1. Dividing the calculated area into an outer area and an inner area according to the material distribution characteristics of the integrated circuit structure; S2. Use the two-dimensional finite element method to analyze and solve the outer region. When solving, the electromagnetic influence of the inner region on it is equivalent to the equivalent current term containing the surface admittance operator, which is added to the electric field partial differential equation obtained by the non-homogeneous Helmholtz equation in the outer region to obtain the electric field equation of the outer region. According to the electric field equation of the outer region, the field distribution of the outer region is obtained; S3. When solving the inner region, according to the field distribution in the outer region, the field distribution at the boundary of the inner region is obtained by the condition of longitudinal electric field continuity, and finally the field distribution in the inner region is obtained by combining the electric field integral equation; Step S2 specifically includes: To solve the field distribution in the outer region, we first use the non-homogeneous Helmholtz equation to obtain the equation of the outer region and boundary that satisfies the standard form of two-dimensional finite element analysis. The equation is as follows: in φ represents the longitudinal electric field in the outer region is the external current source that may exist in the outer region and the boundary, k0 represents the free space wave number, k2 represents the outer region wave number, Z2 represents the outer region wave impedance, ε r2 , μ r2 represents the relative permittivity and relative permeability of the medium in the outer region; The equivalent current term J containing the surface admittance operator is added to the equivalent electromagnetic influence of the inner region on the outer region eq After correction, the final equation of two-dimensional finite element analysis satisfied by the outer region and the boundary is obtained. The final equation is as follows: in is the longitudinal electric field on the boundary of the inner region, κ is the differential surface admittance operator; Step S3 specifically includes: The final equation of the two-dimensional finite element analysis is solved to obtain the longitudinal electric field of the outer region including its boundary. According to the boundary condition of longitudinal electric field continuity, the longitudinal electric field at the inner region boundary is obtained According to the coupling relationship between electric field and magnetic field, The tangential magnetic field on the boundary of the inner region is calculated Finally, the electric field distribution in the inner area is calculated according to the electric field integral equation The calculation is as follows: Where T = 1, represents the Green's function, represents the imaginary unit, represents the second-kind zero-order Hankel function, n′ represents the external normal vector at the source point, Represents an arbitrary vector at the defined source point.

2. The frequency domain electromagnetic calculation method according to claim 1, characterized in that: S1 specifically includes dividing the non-uniform medium, background medium and anisotropic medium in the integrated circuit structure into an outer region, and dividing the high conductivity medium and other structures into an inner region.

Images