As shown in Figure 2, we consider a general circuit structure. Ω denotes the volume of the computational domain, and is bounded by two types of truncation surfaces, including the ports T _i (i = 1, 2, … , p) and the side wall T ₀. In the structure, Tc denotes the surface of a conductor and Td denotes the interface between different dielectric materials. In the computational domain, the electric field E in Ω must satisfy the vector wave equation ∇ × μ r − 1 ∇ × E − k 0 2 ε r E = − j k 0 Z 0 J i (1) where k ₀ and Z ₀ are the wavenumber and the intrinsic impedance of free space; and J _i is an impressed current density as an excitation source.

The boundary conditions applied in the electromagnetic simulation include the perfect electric conductor (PEC) boundary, the impedance boundary condition (IBC), the absorbing boundary condition (ABC) and port boundary condition (PBC). Implementation of these boundaries in the FEM is described in [12].

We adopt the Galerkin method, its scalar multiplication of Eq. 1 by a weighting function w and integration over the computational domain Ω, and obtain the weak form as follows ∇ × w , μ r − 1 ∇ × E Ω − j ω μ 0 n ^ × H , w ∂ Ω − k 0 2 w , ε r E Ω = − j k 0 Z 0 w , J i Ω (2) where the surface integral of two complex-valued vector functions will be shorted conveniently as u , v Γ = ∬ Γ u ⋅ v d S (3) and the volume integral of two complex-valued functions in a domain Ω is denoted by u , v Ω = ∭ Ω u ⋅ v d V (4)

After the finite element analysis of the high-speed interconnection system, the network parameters of the interconnection structure are usually extracted; and the electromagnetic field analysis is converted into a circuit analysis by applying the microwave network theory [18]. In general, the characteristics of a network can be analytically described by admittance parameters (Y-parameters), impedance parameters (Z-parameters), or scattering parameters (S-parameters) [19]. Among them, the S-parameter has been widely used in microwave network and microwave circuit analysis. Because its physical meaning is clear and it can be transformed into other network parameters, it is most suitable to describe the external characteristics of the network.

In order to obtain the modal scattering matrix of the port accurately, we use the TFEM. For the convenience of discussion and without loss of generality, we adopt a two-port structure and the first-order vector basis functions. The main idea of TFEM is to combine the weak form of the Helmholtz equation with the expansion of the mode function in the regular truncated region. When the modal impedance and mode function have been obtained, we analyze the mode function set, ( E _ij , H _ij ) (i = 1 or 2; j = 1, 2, … , N _m ) of the electric field and magnetic field of the port. It is assumed that N _m is sufficient to provide an accurate approximation of the electric field on the port boundaries T ₁ and T ₂, and their electric field mode functions are expressed as follows [20]. E = ∑ j = 1 N m α i j E i j , i = 1   o r   2 (5)

According to the modal orthogonality relationship [21], the transverse electric and magnetic fields are normalized to obtain a set of normalized mode functions ( e _ij , h _ij ) (i = 1 or 2; j = 1, 2, … , N _m ). If the computational domain uses tetrahedral meshes for the region discretization, there are triangular elements on the boundaries. Then the total electric and magnetic field of the entire computational domain can be divided into the following form: E = ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i + e 1 k H = ∑ i = 1 N m α 1 i h 1 i + ∑ i = 1 N m α 2 i h 2 i − h 1 k (6) where N _e is the number of edges in the mesh excluding those on the PEC surfaces, N _pe is the number of edges on the ports, and α _1i and α _2i are the coefficients of the i-th mode. e _1k denotes the excitation on T ₁ for the problem, which is the k-th mode with amplitude 1. The negative sign in front of h _1k is due to the propagation of the excitation wave into the computational domain, and our definition of the propagation directions of the normalized modal electric and magnetic fields is based on the outgoing unit normal along the waveguide axis on the port boundary.

A simplified alternative form of Eq. 6 is obtained by absorbing the incident electric field mode into the modal sum, and the two modified expansions become E = ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i H = ∑ i = 1 N m α 1 i h 1 i + ∑ i = 1 N m α 2 i h 2 i − 2 h 1 k (7)

Without considering an impressed current density, substitution of Eq. 7 into Eq. 2, then imposing the first-order ABC on T ₀ and applying the Galerkin process yields the following finite element equation for j = 1, 2, … , N _e − N _pe : ∇ × w j , μ r − 1 ∇ × ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i Ω + j ω μ 0 Z 0 n ^ × w j , n ^ × ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i T 0 − j ω μ 0 n ^ × ∑ i = 1 N m α 1 i h 1 i + ∑ i = 1 N m α 2 i h 2 i − 2 h 1 k , w j T 1 + T 2 − ω 2 μ 0 ε 0 w j , ε r ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i Ω = 0 (8)

The surface integral on T ₁ and T ₂ is zero, because edge elements on T ₁ and T ₂ have been excluded from the expansion of the electric field vector, and thus n ^ × w _j = 0 on T ₁ and T ₂. Combining the above N _e −N _pe equations in matrix form, we have S + j ω Z − ω 2 T x E + P 1 x 1 α + P 2 x 2 α = 0 (9) where S i , j = ∇ × w i , μ r − 1 ∇ × w j Ω Z i , j = μ 0 Z 0 n ^ × w i , n ^ × w j T 0 T i , j = μ 0 ε 0 w i , ε r w j Ω (10) and the vector x _E contains the expansion coefficients for the N _e −N _pe edge elements; x _1α and x _2α are the vector containing the expansion coefficients of the modal field representation of the tangential electric field on the port boundaries T ₁ and T ₂. The elements of the matrix P ₁ and P ₂ are given by P l i j = ∇ × w i , μ r − 1 ∇ × e l j Ω + j ω μ 0 Z 0 n ^ × w i , n ^ × e l j T 0 − ω 2 μ 0 ε 0 w i , ε r e l j Ω , l = 1   o r   2 (11)

In a similar manner, multiplying the mode functions e _lj (l = 1 or 2; j = 1, 2, … , N _m ) on ports as the weighting function and integrating over the computational domain yields. ∇ × e l j , μ r − 1 ∇ × ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i Ω + j ω μ 0 Z 0 n ^ × e l j , n ^ × ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i T 0 − j ω μ 0 n ^ × ∑ i = 1 N m α 1 i h 1 i + ∑ i = 1 N m α 2 i h 2 i − 2 h 1 k , e l j T 1 + T 2 − ω 2 μ 0 ε 0 e l j , ε r ∑ i = 1 N e − N p e w i E i + ∑ i = 1 N m α 1 i e 1 i + ∑ i = 1 N m α 2 i e 2 i Ω = 0 (12)

Considering the surface integral over T ₁ and making use of the orthogonality of modes, only the term associated with h _1j is nonzero. In matrix form, Eq. 12 become P 1 T x E + R 1 + j ω μ 0 I x 1 α = f P 2 T x E + R 2 + j ω μ 0 I x 2 α = 0 (13) where I is an identity matrix, and the elements of R ₁ and R ₂ are given by R l i j = ∇ × e l i , μ r − 1 ∇ × e l j Ω + j ω μ 0 Z 0 n ^ × e l i , n ^ × e l j T 0 − ω 2 μ 0 ε 0 e l i , ε r e l j Ω , l = 1   o r   2 (14)

Regarding the excitation vector f, all elements are zero except the element in the row associated with the excitation mode k, which has the value f _k = 2jωμ ₀.

Finally, the finite element matrix statement of the problem is obtained by combining Eqs 9, 13, S + j ω Z − ω 2 T P 1 P 2 P 1 T R 1 + j ω μ 0 I 0 P 2 T 0 R 2 + j ω μ 0 I x E x 1 α x 2 α = 0 f 0 (15)

Next extending the tangential part of the modal electric fields according to the N _pe1 and N _pe2 edge elements on the port boundaries T ₁ and T ₂, we have e 1 i = ∑ n = N e − N p e + 1 N e − N p e 2 w n b n 1 i , i = 1,2 , ⋯ , N m e 2 i = ∑ n = N e − N p e 2 + 1 N e w n b n 2 i , i = 1,2 , ⋯ , N m (16) where b n 1 i and b n 2 i are the expansion coefficients in the edge-element expansion of the i-th modal field on T ₁ and T ₂, respectively. Since the transverse modal electric fields are known, these coefficients are easily calculated. Considering all N _m modes, the matrix form expressions of the above equation are e 11 , e 12 , ⋯ , e 1 N m = w N e − N p e + 1 , ⋯ , w N e − N p e 2 b 11 , b 12 , ⋯ , b 1 N m ⏟ B 1 e 21 , e 22 , ⋯ , e 2 N m = w N e − N p e 2 + 1 , ⋯ , w N e b 21 , b 22 , ⋯ , b 2 N m ⏟ B 2 (17) where the vector b ¹ⁱ and b ²ⁱ contain the expansion coefficients in the edge-element expansion of the i-th mode.

Thus the space W of all edge elements used in the finite element approximation is decomposed into three subsets as follows: W = W I ⊕ W P 1 ⊕ W P 2 W I = w 1 , ⋯ , w N e − N p e 1 − N p e 2 W P 1 = w N e − N p e + 1 , ⋯ , w N e − N p e 2 W P 2 = w N e − N p e 2 + 1 , ⋯ , w N e (18) where W _I denotes the subset of the edge elements excluding those on PEC surfaces and ports, and W _P1 and W _P2 denote the subset of the edge elements on T ₁ and T ₂. We define the notation W _q , q∈{I, P ₁, P ₂}, as the row vector containing the elements of subset, and W q T , W q = S W q T , W q + j ω Z W q T , W q − ω 2 T W q T , W q (19)

Thus the expressions for the matrices P ₁, P ₂, R ₁ and R ₂ are given by P 1 = W I T , W P 1 B 1 P 2 = W I T , W P 2 B 2 R 1 = B 1 T W P 1 T , W P 1 B 1 R 2 = B 2 T W P 2 T , W P 2 B 2 (20)

For the case of truncating the wave ports using the TFEM, we obtain the following form of the finite element matrix W I T , W I W I T , W P 1 B 1 W I T , W P 2 B 2 B 1 T W P 1 T , W I B 1 T W P 1 T , W P 1 B 1 + j ω μ 0 I 0 B 2 T W P 1 T , W I 0 B 2 T W P 2 T , W P 2 B 2 + j ω μ 0 I x E x 1 α x 2 α = 0 f 0 (21)

The solution of the above matrix equation can directly obtain the electric field solution inside the electromagnetic device according to x _E and the modal S-parameters of the modal field. The i-th elements of vector x _1α and x _2α are equal to the reflection and transmission coefficient of the i-th mode. According to Eq. 7, the reflection coefficient of the excited mode is obtained by subtracting 1 from the k-th element of x _1α .

A multiconductor port is composed of coupled signal conductors, and a terminal is assigned on a conductor in contact with the port. The number of modes for a port is determined by the number of terminals that touch the port. If N + 1 different conductors touch the port, there are N terminals and one reference conductor commonly called ground.

When the multiconductor port is used for an excitation, a modal scattering matrix can be directly obtained by exciting each port with a number of modes equal to the number of terminals. Then the modal scattering matrix can be converted to a terminal impedance (admittance) matrix and a terminal scattering matrix referenced to either the actual line characteristic impedance matrix or a normalized impedance (for example, 50 Ω).

The modal port representations of the tangential electric and magnetic fields [22] are E t = ∑ n = 1 N a n + b n e t n H t = ∑ n = 1 N a n − b n h t n (22) where a and b are the complex modal coefficient vectors of the ingoing and outgoing modal fields, respectively, and a is related to b by the modal scattering matrix as b = S m a (23) where S _m denotes the modal scattering matrix.

According to the definition of voltage and Ampere’s law, the terminal voltages and currents in matrix notation are given by v = M v a + b a n d i = M i a − b (24) where M _v and M _i are the modal voltage and current matrices. Combining Eqs 23, 24 yields the following relationship i = Y t v (25) where the terminal admittance matrix Y _t is given by Y t = M i I − S m I + S m − 1 M v − 1 (26) and I is identity matrix.

In order to obtain the terminal scattering matrix S _t , we assume that α and β are ingoing and outgoing nodal waves, respectively. α is related to β by the terminal scattering matrix as β = S t α (27)

Using the following definitions

V = Z r α + β a n d   I = Z r − 1 α + β (28) and substituting Eq. 25 into Eq. 28, we obtain β = Z r − 1 Z t − Z r Z t + Z r − 1 Z r α (29) where Z _t = Y _t ⁻¹ is an impedance matrix and Z _r is the reference impedance matrix.

The terminal scattering matrix S _t is given by S t = Z r − 1 Z t − Z r Z t + Z r − 1 Z r (30)

Thus the modal S-parameters are converted into the terminal S-parameters, and other characteristic parameters, such as the terminal Y-parameters and the terminal Z-parameters, can also be obtained.

In theory, the combination of two single-ended conducting traces is a differential pair. The differential pair is an effective way to remove noise from a signal. As shown in Figure 3, a differential signal is transmitted through the differential pair, and the voltages v ₁ and v ₂ on the two traces are equal and opposite. The differential and common voltages v _d and v _c are defined by v d = v 1 − v 2 a n d   v c = v 1 + v 2 / 2 (31)

According to the law of conservation of energy, the differential and common currents i _d and i _c are defined by i d = i 1 − i 2 / 2 a n d   i c = i 1 + i 2 (32)

Eqs 31, 32 can be represented as follows: v = Q e   a n d   i = Q − T u (33) where v = [v ₁ v ₂] ^T , e = [v _d v _c ] ^T , i = [i ₁ i ₂] ^T , u = [i _d i _c ] ^T , and Q = 1 / 2 1 − 1 / 2 1 (34)

Substituting Eq. 33 into Eq. 25 yields u = Q T Y t Q e (35) and the differential admittance matrix Y d t = Q T Y t Q (36)

In a similar way, we can obtain the differential impedance matrix Z d t = Q Z t Q T (37) and the differential scattering matrix S d t = Z d r − 1 Z d t − Z d r Z d t + Z d r − 1 Z d r (38) where Z _dr is the differential reference impedance matrix.

The first example is a coaxial cable as shown in Figure 4 to demonstrate that the developed tool has basic terminal network parameter analysis capabilities. The copper conductor is surrounded by a homogeneous dielectric, whose relative permittivity is 4.4. The radii of the copper conductor and the dielectric are 1.5 mm and 3 mm. The length of the coaxial cable is 10 mm. The operation frequency is 7 GHz. The frequency scan starts at 6 GHz and ends at 8 GHz.

The calculation results of the proposed method (TFEM) and FEM are compared with HFSS. The result comparisons of terminal network parameters diagram of the coaxial cable are shown in Figure 5. From it we can see that the results calculated by the TFEM and FEM are consistent with HFSS, thereby verifying the accuracy of the developed tool.

In addition, the computational performance of the TFEM and FEM are investigated. As shown in Table 1, it can be seen that the calculation time and peak memory of the TFEM are less than those of the FEM, which proves the advantages of the proposed method.

The second example is a differential microstrip line as shown in Figure 6A, which is used to demonstrate that the developed tool has the basic terminal network parameter analysis capability of differential pairs. The differential microstrip line comprises of two copper traces and a dielectric substrate enclosed in an air-box. The relative permittivity of this substrate is 4.4. The thickness of the trace is 0.7 mil, and the width of the model is 2 mil. The remaining dimensions are indicated in Figure 6A in mils. The finite conductivity boundary is assigned on the bottom face of the model, and its parameters are the same as those of cooper. The radiation boundary is assigned on the top, left and right faces.

The two ports are assigned in the remaining faces and deembedded with 499 mil outwards from the structure. There are 4 terminals assigned as shown in Figure 6B, and their renormalizing impedance is 50 Ω. The terminals 1_a and 1_b form a differential pair named D1, the terminals 2_a and 2_b form another differential pair named D2. The differential mode impedance of the differential pairs is 100 Ω. The operation frequency is 3 GHz, and the frequency sweep is from 2 GHz to 4 GHz.

The computational results of the terminal network parameters of differential pairs are shown in Figure 7. It is seen that when the operating frequency is 3.4 GHz, the return loss of the differential signals is the lowest. It can be observed that the results calculated by the developed tool show good agreement with HFSS so that demonstrate the terminal network parameter analysis capabilities of differential pairs of the developed tool.

In order to demonstrate the multi-mode terminal network parameter analysis capability of the developed tool, the third example is a set of seven copper striplines as shown in Figure 8. The first, fourth and seventh stripline from the left are ground conductors. The relative permittivity of this substrate is 4.4. The thickness of the trace is 0.7 mil, and the width of the model is 2 mil. The remaining dimensions are indicated in Figure 8 in mils. The finite conductivity boundary, whose parameters are the same as those of cooper, is assigned on the bottom and top faces of the model.

We assign two ports on the front and back faces and deembed them with 499 mil outwards from the structure. There are eight terminals assigned as shown in Figure 9, and their renormalizing impedance is 50 Ω. The terminals form four differential pairs, and their differential mode impedance is 100 Ω. The operation frequency is 3 GHz, and the frequency sweep is from 2 GHz to 4 GHz.

As shown in Figure 10, when the operating frequency is 2.8 GHz, the return loss of the differential signals is the lowest. It is observed that the results of the terminal network parameters of differential pairs calculated by the developed tool are in good agreement with HFSS, which can demonstrate the correctness of the developed tool for multi-mode terminal network parameter analysis.

The above examples verify the accuracy of the developed tool, then a model with more complex structure is used to further verify the efficiency improvement of the developed tool. The final example is a via model as shown in Figure 11, which is widely used for signal connection between different layers in multilayer PCBs. It has a pair of microstrip lines that transition through the vias to a pair of striplines on a lower layer. The width of the structure is 128 mil. The conductors are copper and the relative permittivity of the dielectric is 4.4. The thicknesses of the conductors and dielectric are 1.2 mil and 5 mil, respectively. The remaining dimensions are indicated in Figure 11 in mils. The radiation boundary is applied to the air box.

The two ports are assigned as shown in Figure 11 and deembedded with 48 mil outwards from the structure. The two microstrip lines are each assigned a terminal in the coupled microstrip port. The same goes for the two striplines at the opposite end. Their renormalizing impedance is 200 Ω. The operation frequency is 4 GHz, and the frequency sweep is from 2 GHz to 6 GHz.

Figure 12 shows comparisons of the calculation results of the terminal parameters calculated by the developed tool and HFSS. It can be seen the results are consistent, which further demonstrates the accuracy of the developed tool.

We further analyze the computational efficiency of the via model. The comparisons of the developed tool and HFSS are given in Table 2. The developed tool has advantages in terms of CPU times and memory consumption with a similar number of tetrahedron grids in the via model simulation, thus proving the efficiency of the developed tool.