Skin Effect Modeling in Conductors of Arbitrary Shape Through a
Surface Admittance Operator and the Contour Integral Method

arXiv:1509.08357v1 [] 28 Sep 2015

Utkarsh R. Patel and Piero Triverio

∗†

Submitted to IEEE Transactions on Microwave Theory and Techniques
on September 27, 2015

Abstract
An accurate modeling of skin effect inside conductors is of capital importance to solve transmission
line and scattering problems. This paper presents a surface-based formulation to model skin effect in
conductors of arbitrary cross section, and compute the per-unit-length impedance of a multiconductor
transmission line. The proposed formulation is based on the Dirichlet-Neumann operator that relates
the longitudinal electric field to the tangential magnetic field on the boundary of a conductor. We
demonstrate how the surface operator can be obtained through the contour integral method for conductors
of arbitrary shape. The proposed algorithm is simple to implement, efficient, and can handle arbitrary
cross-sections, which is a main advantage over the existing approach based on eigenfunctions, which is
available only for canonical conductor’s shapes. The versatility of the method is illustrated through
a diverse set of examples, which includes transmission lines with trapezoidal, curved, and V-shaped
conductors. Numerical results demonstrate the accuracy, versatility, and efficiency of the proposed
technique.

1

Introduction

Transmission line modeling is a crucial part of the computer aided design of a variety of systems, including
high-speed electronic boards [1, 2], integrated circuits [3], microwave systems [4], metamaterials [5], and
power grids [6]. In order to create a transmission line model, one must first obtain the per-unit-length
(p.u.l.) impedance and admittance of the line [7]. In several applications, such as the investigation of signal
integrity issues in high-speed electronic systems, such parameters are required over a wide band, extending
from DC up to tens of gigahertz [1]. Across such band, the development of skin effect significantly changes
the line behaviour [2], and must be accurately described. Skin effect modeling is also important in scattering
problems [8].
In order to obtain the p.u.l. impedance of a transmission line, one must describe the electromagnetic
fields both inside the conductors (interior problem) and outside (exterior problem). For the exterior problem,
integral equations [9–11] are commonly utilized. By requiring only a discretization of the conductors’ surface,
integral equations are typically more efficient that volumetric approaches such as the finite difference or finite
element method (FEM) [12, 13], which require the discretization of a large area around the conductors.
For the interior problem, which is responsible for capturing skin effect, both volumetric and surface
methods have been used. Volumetric approaches include the FEM [13, 14], conductor partitioning [15] and
volumetric integral equations [16]. Unfortunately, as frequency increases, these methods become computationally inefficient, since a very fine mesh is needed to capture the pronounced skin effect.
Surface methods solve the interior problem using only a discretization of the boundary of the conductors.
This is achieved by describing the electromagnetic behaviour of the conductor through a surface operator that
∗ This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (Discovery Grant
program) and in part by the Canada Research Chairs program.
† U. R. Patel and P. Triverio are with the Edward S. Rogers Sr. Department of Electrical and Computer Engineering,
University of Toronto, Toronto, M5S 3G4 Canada (email: utkarsh.patel@mail.utoronto.ca, piero.triverio@utoronto.ca).

1

relates the electric and magnetic fields on the boundary. This operator is responsible for modeling skin and
proximity effects [2]. The surface operator can be obtained analytically or numerically. In [17, 18], a surface
operator is derived using the so-called surface-impedance boundary conditions [19], under the approximation
of small curvature and well-developed skin effect [20]. While high-order boundary conditions have been
proposed to increase accuracy [21], this approach is mostly suitable for conductors with smooth boundaries.
Numerically, the surface operator can be obtained using finite differences [22], finite elements [23], the electric
field integral equation [24, 25] or the magnetic field integral equation [20]. All these approaches require, to
calculate the surface operator, a volumetric discretization of the interior problem and the calculation of
several kernel matrices whose size depends on mesh size. These features increase their computational cost
and complexity.
An interesting technique to derive the surface operator analytically was introduced by Knockaert and
De Zutter in [26]. Using the eigenfunctions of the Helmholtz equation, the surface admittance operator
is obtained analytically, avoiding a discretization of the interior problem. This key idea leads to a simple
and efficient formulation. However, since eigenfunctions can realistically be obtained only for canonical
geometries, this approach has been restricted so far to circular [26, 27], tubular [28], rectangular [26], and
triangular conductors [29].
In this paper, we show how the surface admittance operator for conductors of arbitrary shape can be
efficiently obtained from a simple contour integral, through the so-called contour integral method [30]. While
the possibility to extract the surface admittance operator from the contour integral is already known [31], it
has been applied only to simple scattering problems [31]. All works that followed adopted the eigenfunctions
approach [26, 29, 32–34] which is feasible for canonical shapes only. While arbitrary shapes can be decomposed into triangles [29], this significantly adds to the algorithm complexity. This paper provides a simple,
unified way to handle arbitrary cross-sections, and extract accurate and broadband p.u.l. parameters for a
variety of transmission lines. Numerical tests show that the proposed method is robust and computationally
efficient, even when compared against the analytical approach based on eigenfunctions [26]. In this paper, the
proposed method is applied to calculate transmission line parameters, but it can also be used for scattering
problems [35].
The paper is organized as follows. In Sec. 2, we define the problem, and in Sec. 3 we discuss how the
surface admittance operator is related to the contour integral method from a theoretical standpoint. In
Sec. 4, we discuss the numerical implementation of the proposed approach, and in Sec. 5 how it can be used
to compute the p.u.l. impedance of arbitrary transmission lines. Finally, in Sec. 6, we demonstrate the
accuracy, robustness, and computational efficiency of the proposed method through a comprehensive set of
examples.

2

Problem Definition

We consider a system of P conductors of arbitrary shape, having conductivity σ, permittivity ε, and permeability µ. For the sake of generality, the conductors are assumed to be inside a stratified medium, where
each layer has permittivity εl and permeability µl . If a conductor extends into two layers, it is decomposed
into two parts, as discussed in [36]. Our final goal is to calculate the partial p.u.l. resistance R(ω) and
inductance L(ω) matrices in the Telegrapher’s equation [2]
∂V
= − [R(ω) + jωL(ω)] I .
(1)
∂z

T

T
In the equation above, V = V1 . . . VP collects the potential Vp of each conductor, while I = I1 . . . IP
contains the current Ip flowing in each conductor. The impedance parameters will be calculated with the
surface admittance approach of [26] using, however, a contour integral to obtain the surface admittance
operator.

2

Top layers

µl , εl

µl , εl

t̂0

σ, µ, ε

0

n̂

~r 0

O

d~

~r
n̂

y

Top layers
z

x

µl , εl

γ (p)

γ (p)
(p)

Js (~r)

t̂

Bottom layers

Bottom layers

Figure 1: Left panel: sample geometry of a conductor with arbitrary cross-section inside a stratified medium.
Right panel: equivalent configuration obtained after replacing the conductor with the surrounding medium,
(p)
and introducing an equivalent current density Js (~r) on contour γ (p) .

3
3.1

Surface Admittance Operator Through the Contour Integral
Method
Surface Admittance Formulation

We discuss the surface admittance formulation used to solve the interior problem by considering the p-th
conductor of a multiconductor transmission line. The conductor has arbitrary cross section, has a simplyconnected contour γ (p) , and is depicted in the left panel of Figure 1. The position vector of an arbitrary
point on γ (p) is denoted by ~r. Finally, the unit vectors normal and tangential to contour γ (p) are denoted
by n̂ and t̂, respectively, as shown in Fig. 1.
In order to derive the surface admittance operator, we initially follow [26], and consider a boundary value
(p)
problem. Suppose that the z-directed electric field on contour γ (p) is given by Ez (~r). Under the quasi(p)
(p)
can be obtained from the scalar Helmholtz
TM assumption [26], the electric field Ez inside contour γ
equation [37]
∇2 Ez(p) + k 2 Ez(p) = 0
(2)
subject to the boundary condition
Ez(p) (~r) = Ez(p) (~r)

~r ∈ γ (p) .

(3)

p

In (2), k = ωµ (ωε − jσ) is the wavenumber inside the conductor. Furthermore, the tangential mag(p)
netic field Ht (~r) along the contour γ (p) follows directly from Maxwell’s equations under the quasi-TM
assumption [37]
"
#
(p)
1
∂Ez (~r)
(p)
Ht (~r) =
,
(4)
jωµ
∂n
(p)
~
r ∈γ

where the derivative is taken along the direction normal to contour γ (p) , as shown in the left panel of Fig. 1.
Next, we replace the conductor by the material of the surrounding layer [26], as shown in the right panel
of Fig. 1. From here onwards, we will call this the equivalent configuration. Due to this modification, the
(p)
electric field inside the conductor changes to Eez (~r) which satisfies the Helmholtz equation
2 e(p)
∇2 Eez(p) + kout
Ez = 0

(5)

subject to the same Dirichlet boundary condition (3)
Eez(p) (~r) = Ez(p) (~r)

~r ∈ γ (p) .

√

(6)

In (5), kout = ω µl εl is the wavenumber inside the layer surrounding the conductor. In this new configuration, the tangential magnetic field along γ (p) is
"
#
ez(p) (~r)
1
∂
E
(p)
e t (~r) =
H
.
(7)
jωµl
∂n
(p)
~
r ∈γ

3

By replacing the conductor with the surrounding medium, we have modified the fields both inside and
outside γ (p) . Hence, to restore the original fields outside γ (p) , we invoke the equivalence principle [38] and
(p)
introduce an equivalent surface current density Js (~r) on γ (p) [26], as shown in the right panel of Fig. 1.
The equivalent current is directed along z, and is given by
(p)

(p)

e t (~r) .
Js(p) (~r) = Ht (~r) − H

(8)

It is important to note that an equivalent magnetic current on γ (p) is not required, since the electric field
on γ (p) in the original and equivalent configuration remain the same due to the Dirichlet conditions (3) and
(6). By substituting (4) and (7) into (8), we obtain [26]
"
#
(p)
(p)
1 ∂ Eez (~r)
1 1 ∂Ez (~r)
(p)
.
(9)
−
Js (~r) =
jω µ ∂n
µl ∂n
(p)
~
r ∈γ

(p)

Equation (9) defines a surface admittance operator Ys that relates the longitudinal electric field and
equivalent current on γ (p) as
Js(p) (~r) = Ys(p) Ez(p) (~r) .
(10)
(p)

An explicit expression for Ys can be written in terms of the eigenfunctions of the Helmhotz equations (2)
and (5), as shown in [26]. However, this approach is viable only for canonical conductor shapes, for which
eigenfunctions are known analytically [26–29]. For arbitrary shapes, eigenfunctions can only be computed
numerically. Since many eigenfunctions are needed to accurately model the operator, this approach can
be very time consuming, and is typically avoided. In the next section, we show that the contour integral
method [30] provides an efficient and robust way to compute such operator numerically for arbitrary shapes.

3.2

Contour Integral Method

Equating (9) and (10), we see that, in order to derive an explicit expression for the surface admittance
(p)
operator, we need a relation between the electric field Ez on the boundary and its normal derivative
(p)
∂Ez /∂n. With this goal in mind, we reconsider the way we solve (2). The contour integral method [30]
gives the solution of the Helmholtz equation in terms of the electric field and its normal derivative on γ (p)
"
‰
∂G(~r, ~r 0 ) (p) 0
j
(p)
Ez (~r )
Ez (~r) =
2 γ (p)
∂n0
#
(p)
0
(~
r
)
∂E
z
− G(~r, ~r 0 )
dr0
(11)
∂n0
where ~r and ~r 0 are both on γ (p) , and n0 is the unit vector normal to the contour at the point ~r 0 . The Green’s
function is
G(~r, ~r 0 ) = C0 J0 (kd) − jY0 (kd)
(12)
where J0 (.) and Y0 (.) are the zero-th order Bessel and Neumann functions [39], respectively. As shown in
the left panel of Fig. 1, the distance between points ~r and ~r 0 is denoted as
d~ = ~r 0 − ~r ,

(13)

and d = d~ . Equation (11) stems from the planar Green’s theorem, and shows that the electric field along
the contour γ (p) can be interpreted as the superposition of cylindrical waves that originate from points along
γ (p) . Constant C0 can be any complex number. If C0 = 1 the Green’s function is the Hankel function of the
second kind, which represents outgoing cylindrical waves. In Sec. 4.5, we will discuss how to choose C0 at
low and high frequency to achieve high numerical robustness. The contour integral equation (11) relates the
electric field on the boundary to its normal derivative. It can thus be used, after numerical discretization,
to derive an explicit expression for the surface admittance operator (10).
4

4
4.1

Numerical Formulation
Discretization of Electric Fields, Magnetic Fields, and Equivalent Current

We discretize the contour integral equation with the method of moments [10], using point matching to test
(p)
the resulting equation. We divide contour γ (p) into Np segments γn , and expand the longitudinal electric
field in terms of pulse basis functions
Ez(p) (~r) =

Np
X

en(p) Π(p)
r) ,
n (~

(14)

n=1
(p)

where Πn (~r) is the n-th pulse basis function which is one if ~r belongs to the n-th segment, and zero
otherwise. We also define ~rn to be the position vector of the midpoint of the n-th partition. Similarly, we
(p)
e t(p) (~r) that appear in (8) as
discretize the tangential magnetic fields Ht (~r) and H
(p)

Ht (~r) =

Np
X

(p)
h(p)
r) ,
n Πn (~

(15)

h̃n(p) Π(p)
r) .
n (~

(16)

n=1

and
e t(p) (~r) =
H

Np
X
n=1

For simplicity of notation, we cast all expansion coefficients into column vectors
h
iT
E(p) = e1(p) e2(p) . . . e(p)
,
Np
iT
h
(p)
,
H(p) = h1(p) h2(p) . . . hN
p
h
iT
e (p) = h̃(p) h̃(p) . . . h̃(p)
H
.
1
2
Np

(17)
(18)
(19)

(p)

Similarly, we discretize equivalent current Js (~r) along the contour γ (p) using pulse basis functions as
Js (~r) =

Np
X

jn(p) Π(p)
r) ,
n (~

(20)

n=1
(p)

with coefficients jn stored into vector
h
J(p) = j1(p)

(p)

j2

...

(p)

jNp

iT

.

(21)

From (8), we have the following relation between the coefficients of equivalent current and magnetic fields
e (p) .
J(p) = H(p) − H

(22)

e (p) to the electric field coefficients E(p) via the contour integral equation (11) to
Next, we relate H(p) and H
obtain the surface admittance operator.

4.2

Magnetic Field in the Original Configuration

After combining (4) and (11), we substitute the electric and magnetic field expansions (14) and (15) into
the contour equation to obtain
"
‰
Np
Np
X
∂G(~r, ~r 0 ) X (p) (p) 0
j
(p) (p)
em Πm (~r) =
en Πn (~r )
2 γ (p)
∂n0
m=1
n=1
#
Np
X
0
(p) (p)
0
−jωµG(~r, ~r )
hn Πn (~r ) dr0
(23)
n=1

5

Using point-matching [38], we test the equation above at the midpoints ~rm of all Np segments, obtaining
ˆ
Np
j X (p)
∂G(~rm , ~r 0 ) 0
dr
en
(p)
2 n=1
∂n0
γn

e(p)
m =

+

ˆ
Np
ωµ X (p)
hn
G(~rm , ~r 0 )dr0
(p)
2 n=1
γn

(24)
(p)

for m = 1, . . . , Np . Note that in (24), the integration is performed only over the n-th segment γn . All
relations (24) can be compactly written in matrix form as
UE(p) = PH(p) ,

(25)

where U and P are square matrices with dimension Np × Np . Element (m, n) of matrix U, if m 6= n, is given
by
[U]m,n

jk
=
2

ˆ
(p)

γn

d~m · n̂0
[C0 J1 (kdm ) − jY1 (kdm )] dr0 ,
dm

(26)

where d~m = ~r 0 − ~rm , dm = d~m , and we used the fact that the normal derivative of (12) can be written as
d~ · n̂0
∂G(~r, ~r 0 )
= −k
[C0 J1 (kd) − jY1 (kd)] .
0
∂n
d

(27)

The diagonal entries of U are given by
[U]m,m = 1 .
The (m, n)-th entry of P in (25) is given by
ˆ
ωµ
[P]m,n =
[C0 J0 (kdm ) − jY0 (kdm )] dr0 .
2 γn(p)

(28)

(29)

In all numerical examples of Sec. 6, the integrals in (26) and (29) were evaluated using a 5-point Gaussian
quadrature routine. Since the Neumann function approaches infinity for small arguments [39], the diagonal
entries of P must be evaluated analytically using the small-argument approximation of the Bessel functions,
as done in [30].
From (25), we can express the magnetic field H(p) in terms of the electric field on the same conductor as
H(p) = P−1 UE(p) .

4.3

(30)

Magnetic Field in the Equivalent Configuration

Next, we find the magnetic field in the equivalent configuration by repeating all the steps of Sec. 4.2, but with
the material parameters of the conductor replaced by the parameters of the surrounding medium. These
steps lead to
e (p) = P−1 Uout E(p) ,
H
(31)
out
which is analogous to (30). In the equation above, entries of Pout and Uout are calculated with formulas
(29), (26), and (28), with wavenumber k replaced by kout and permeability µ replaced by µl .

4.4

Surface Admittance Operator

Finally, by substituting (30) and (31) into (22), we obtain
J(p) = Y(p) E(p)
6

(32)

where
Y(p) = P−1 U − P−1
out Uout

(33)

is the discretized surface admittance operator of the p-th conductor. Therefore, we see that with the proposed
contour integral approach, it is sufficient to evaluate matrices P, Pout , U, and Uout to easily obtain the
surface admittance operator for a conductor of arbitrary shape.
In order to compact the notation, we collect the coefficients of the electric field and equivalent current of
all conductors into two column vectors E and J. The surface admittance relations (32) can be thus compactly
written as
J = Ys E ,
(34)
where

 (1)
Y


Ys = 



Y(2)
..



.


.

(35)

Y(P )

4.5

Choice of C0

In this section, we discuss how to set the constant C0 at low and high frequency in order to achieve a
well-conditioned algorithm.
4.5.1

Low Frequency

At low frequency, where skin effect has not yet developed, computing (29) requires the evaluation of the
Green’s function (12) for very small arguments. For small arguments, the Neumann function Y0 (.) dominates
the Bessel function J0 (.) [39]. Hence, to maintain a good numerical contrast between the Bessel and Neumann
functions, we must set C0 to a high value, as discussed in [40]. The value of 106 provided accurate results
for all numerical tests we performed, including the examples of Sec. 6.
4.5.2

High Frequency

At high frequency, we set C0 = 1. This makes the Green’s function (12) become the Hankel function of the
second kind, which is well-behaved at high frequency.
4.5.3

Switching Condition

As previously discussed, two different values of C0 are appropriate at low and high frequency. Numerical
tests demonstrated that there is a large range of intermediate frequencies where both values of C0 provide
accurate results. Therefore, the choice of the frequency where one should switch from one C0 value to the
other is not critical. In our implementation, we use the low-frequency value of C0 on the p-th conductor
when
∆p
≤ t,
(36)
δp
where ∆p and δp are, respectively, the minimum transversal dimension and the skin depth in the p-th
conductor. In all numerical test we performed, any t value between 0.2 and 0.5 gave good results, and t = 0.5
has been used in all numerical examples presented in Sec. 6. When frequency increases and inequality (36)
no longer holds, the high-frequency value of C0 is used.

5

Exterior Problem and Impedance Computation

In this section, we use the electric field integral equation to express the relation between the electric field
on the conductors’ boundary, and the equivalent currents dictated by the region outside the conductors [26].

7

Combined with the surface admittance operator (34), this will lead to the p.u.l. impedance of the transmission
line. On the contour of the p-th conductor, the electric field integral equation [37] reads
Ez(p) (~r) = jωµl

P ‰
X

Js(q) (~r 0 )G0 (~r, ~r 0 )ds0 −

γ (q)

q=1

∂Vp
∂z

(37)

where ~r ∈ γ (p) and the integration is performed over the contour γ (q) of each conductor (q = 1, . . . , P ).
The integral kernel G0 (~r, ~r 0 ) is the Green’s function of the surrounding medium. Since, through the surface
admittance formulation, all conductors have been replaced with the surrounding medium, G0 (~r, ~r 0 ) is simply
the Green’s function of a stratified medium [41]. Following [16,26], we substitute (1), (14) and (20) into (37),
obtaining
Np
X

(p)
e(p)
r) =jωµl
m Πm (~

m=1

P ‰
X
q=1

+

P
X

Nq
X

γ (q)

jn(q) Π(q)
r 0 )G0 (~r, ~r 0 )ds0q
n (~

n=1

[Rpq (ω) + jωLpq (ω)] Iq

(38)

q=1

for p = 1, . . . , P . Using point-matching [38], this integral equation can be converted into the system of
algebraic equations
E = jωµl G0 J + Q [R(ω) + jωL(ω)] I
(39)
where G0 is the Green’s matrix. In (39), Q is a block diagonal matrix


11


12


Q=

.
..


1P

(40)

where 1p is a vector of size Np × 1 whose entries are all ones. Furthermore, as shown in [26], the total current
inside each conductor is equal to the contour integral of the equivalent current (20), and hence
I = QT WJ

(41)

where W is a diagonal matrix in the form
 (1)
w


W=



w(2)
..






.

(42)

w(P )
where w(p) is a diagonal matrix of size of Np × Np , where entry (n, n) is the width of the n-th segment of
γ (p) .
As shown in [27], we can manipulate (39) using (34) and (41) to obtain the partial p.u.l. resistance
h
i−1
−1
R(ω) = QT W (1 − jωµl Ys G0 ) Q

(43)

h
i−1
−1
L(ω) = ω −1 QT W (1 − jωµl Ys G0 ) Q
.

(44)

and partial p.u.l. inductance

The p.u.l. impedance can be obtained from the partial p.u.l. impedance by taking one of the conductors as
reference.
8

3.5
1

2

2
0.2
0.5
0.2
2

Figure 2: Cross section of the two-conductor lines analyzed in Sec. 6.1. All dimensions are in millimeters.
Table 1: Example of Sec. 6.1: CPU time required to compute the p.u.l. impedance at one frequency
Test Case
Proposed MoM-SO [27] De Zutter [26]
Round conductors
0.04 s
0.0006 s
N/A
Rectangular conductors
0.16 s
N/A
0.12 s

6

Numerical Results

In this section, we demonstrate the efficiency and robustness of the proposed method, and compare it
against the state of the art. We illustrate the versatility of the technique by considering a comprehensive
set of transmission lines made by circular, rectangular, trapezoidal, V-shaped, and conformal conductors.
The method is validated against the eigenfunctions approach [26, 27], as well as a completely independent
technique (FEM) [42]. All computations were performed on a computer with 16 GB of memory and a
3.4 GHz processor. All techniques based on a surface admittance operator were implemented in MATLAB.

6.1
6.1.1

Two-conductor lines
Round conductors

We first consider a transmission line made by two round conductors with radii a1 = 1 mm, a2 = 2 mm,
and spacing d = 3.5 mm. The line cross section is shown in the left panel of Fig. 2. The conductivity of
both conductors is σ = 5.8 · 107 S/m. We calculated the p.u.l. impedance using the proposed technique
and the MoM-SO algorithm [27], which uses the eigenfunctions method to derive the surface admittance
operator. Figure 3 shows the p.u.l. resistance and inductance obtained with the proposed approach and
with MoM-SO. The two methods are in excellent agreement. For the proposed approach, the boundaries
enclosing the small and large conductors were discretized with N1 = 28 and N2 = 60 pulse basis functions,
respectively. As shown in Table 1, for each frequency point the proposed approach took only 0.04 s. In this
case, the eigenfunctions approach is faster since it can capture, with few Fourier basis functions, the field
distribution inside the circular conductors. The advantage of the proposed approach is generality: it can be
applied to arbitrary shapes, while MoM-SO is limited to round conductors, either solid [26,27] or hollow [28].
6.1.2

Rectangular conductors

We now consider the case of two rectangular conductors presented in [26]. Each conductor has conductivity
σ = 5.6 · 107 S/m and dimension 2 mm × 0.2 mm, as shown in the right panel of Fig. 2. Figure 4 shows the
p.u.l. resistance and inductance for this transmission line computed with the proposed method and with [26].
The latter method computes the surface admittance operator analytically using sinusoidal functions, which
are the eigenfunctions of the Helmholtz equation on a rectangular domain. In order to interface the surface
admittance operator to the electric field integral equation (38), the sinusoidal functions have to be mapped
onto pulse basis functions [26]. From Fig. 4, we see that the resistance and inductance values obtained
with both approaches match very well. For the proposed approach, each conductor was discretized with
N1 = N2 = 106 pulse basis functions. For the method of [26], the number of sinusoidal harmonics were
chosen to be M = 400 (see [26] for definition of M ), and each conductor was discretized with N1 = N2 = 280
pulse basis functions. The approach of [26] required more basis functions than the proposed approach to
get accurate results at very high frequency. The computational time taken by both techniques is given in
9

Resistance p.u.l. [Ω/m]

MoM−SO
Proposed

0

10

−1

10

−2

10

0

10

2

10

4

10
Frequency [Hz]

6

10

8

10

Inductance p.u.l. [µH/m]

0.5
0.45
0.4
0.35
0.3
0.25
0.2
0
10

MoM−SO
Proposed
2

10

4

10
Frequency [Hz]

6

10

8

10

Figure 3: P.u.l. resistance and inductance of the transmission line with two round conductors considered in
Sec. 6.1.1, obtained with the proposed method and MoM-SO [27].
Table 2: Examples of Sec. 6.2, 6.3 and 6.4: CPU time required to compute the p.u.l. impedance at one
frequency with the proposed technique and the FEM [42].
Section
Test Case
Proposed FEM
Sec. 6.2
Valley microstrip
0.52 s
73.8 s
Sec. 6.3
On-chip line
0.18 s
9.11 s
Sec. 6.4 Curved microstrips
0.51 s
40.3 s

Table 1. We see that the computational cost of the proposed approach is not far from the cost of [26],
although the proposed method is more general, as it can handle arbitrary shapes.

6.2

Valley Microstrip Line

We consider a valley microstrip line from [43], having the cross-section depicted in Fig. 5. This line type is
used in low-loss microwave integrated circuits, and has been considered in several previous publications [44,
45]. All conductors are made of copper (σ = 5.8 · 107 S/m). The top conductor is the signal line, while the
two lower conductors form the reference line.
Figure 6 shows the p.u.l. inductance and resistance of the system obtained with the FEM [42, 46] and
with the proposed approach. The results from the two methods agree well both at very low (1 MHz) and
very high (1 THz) frequency. This test further validates the proposed method and shows its numerical
robustness. In the FEM simulation, we had to use a fine mesh with 3,102 boundary elements and 188,107
triangular elements, in order to properly resolve the pronounced skin effect at high frequency. In the proposed
technique, the boundary of the signal conductor was discretized with N1 = 194 pulse basis functions, while
each ground conductor was discretized with N2 = N3 = 92 pulse basis functions. The CPU time taken by
both techniques is reported in Table 2 and shows the efficiency of the proposed method, which requires only

10

De Zutter et al.
Proposed

1

Resistance p.u.l. [Ω/m]

10

0

10

−1

10

2

4

Inductance p.u.l. [µH/m]

10

10

6

8

10
Frequency [Hz]

10

10

10

De Zutter et al.
Proposed

0.28
0.26
0.24
0.22
0.2
2
10

4

10

6

8

10
Frequency[Hz]

10

10

10

Figure 4: P.u.l. resistance and inductance of the transmission line with two rectangular conductors considered
in Sec. 6.1.2. Results obtained with [26] are labeled as “De Zutter et al”.
7.5

1

1

7.5

3

13
10

1

13

1

Figure 5: Valley microstrip line considered in Sec. 6.2. All dimensions are in micrometers.
0.52 s to extract the p.u.l. impedance at one frequency for this non-trivial line. While in such cases one
can resort to a general 2D FEM approach, this significantly increases the computational cost, as shown in
Table 2.

6.3

On-chip Transmission Line with Trapezoidal Conductors

Next, we consider the four-conductor transmission line [47] shown in Fig. 7. All conductors are made of
Al-oxide (σ = 3.57·107 S/m). This type of transmission line is typical for interconnects in integrated circuits.
The trapezoidal shape of the signal lines arises is cause, for example, by underetching or electrolytical growth.
As shown in [47], approximating signal lines with perfect rectangles results in a non-negligible error on the
p.u.l. resistance and inductance. Figure 8 shows various entries of the p.u.l. resistance and inductance
matrices of the transmission line, computed with the proposed technique and the FEM [42]. The proposed
technique correctly captures the non-trivial impedance behaviour over frequency. In the proposed method,
each trapezoidal conductor was discretized with 40 pulse basis functions, while 116 functions were used for
the lower conductor. In the FEM simulation, 27,596 triangular and 718 boundary elements were used to

11

Resistance p.u.l. [Ω/m]

FEM
Proposed

4

10

3

10

6

8

10

10

10

12

10

10

Frequency [Hz]

Inductance p.u.l. [µH/m]

0.28
FEM
Proposed

0.27
0.26
0.25
0.24
0.23
0.22
6
10

8

10

10

12

10

10

Frequency [Hz]

Figure 6: P.u.l. resistance and inductance of the valley microstrip line considered in Sec. 6.2, computed with
the proposed method and the FEM [42].

#3
ε0 , µ0

0.75
#1
1.5 0.5
#0

#2

1
1
1

13

Figure 7: On-chip interconnect of Sec. 6.3. All dimensions are in micrometers.
mesh the line’s cross section. Table 2 shows that the proposed method took just 0.18 s to calculate the whole
p.u.l. impedance matrix at a single frequency.

6.4

Curved Microstrip Lines

Flexible dielectrics allow for the creation of curved interconnects conformal to a cylindrical surface [48, 49].
We analyse the configuration shown in Fig. 9, which features four copper conductors (σ = 5.8 · 107 S/m).
Conductors were chosen to be very thin and wide to demonstrate the robustness of the proposed approach in
handling conductors with large aspect ratio. Figure 10 shows selected entries of the resistance and inductance
matrices, computed with the proposed method and the FEM [42]. A close match can be observed, which
validates the proposed technique. Because of the low skin depth at high frequency, the FEM undergoes
some numerical difficulties at high frequency in the computation of inductance, that we could not eliminate
by refining the mesh. This issue is indicative of the intrinsic challenge in meshing the region inside the
conductors at high frequency, where current is confined in a small region near the edges. Surface methods,
such as [26] and the one in this paper, completely avoid this issue. In the proposed method, a total of 373
pulse basis functions were used to discretize the conductors boundary. Table 2 shows that, with the proposed

12

Resistance p.u.l. [Ω/m]

FEM
Proposed

5

10

R11
4

10

R12
R23

3

10

7

Inductance p.u.l. [µH/m]

10

8

10

0.5

L22

0.4

L11

9

10

9

10

10
10
Frequency [Hz]

11

10

12

10

0.3
0.2
0.1
0
7
10

L12
L32
8

10

10
10
Frequency [Hz]

11

10

12

10

Figure 8: P.u.l. resistance and inductance of the on-chip transmission line with trapezoidal conductors
considered in Sec. 6.3.
method, the p.u.l. impedance of the conformal interconnect can be obtained in 0.51 s per frequency.

7

Conclusion

We presented an efficient approach to model skin effect in conductors of arbitrary cross section. The proposed
approach is based on a surface admittance operator. We show how the operator can be efficiently computed
from a contour integral. The proposed approach can handle conductors of arbitrary shape, unlike the
popular eigenfunctions method which is viable only for canonical shapes, where eigenfunctions are available
analytically. The novel method has been applied to the computation of the resistance and inductance
of transmission lines with rectangular, circular, trapezoidal, V-shaped and curved conductors. Numerical
results demonstrate the efficiency, robustness and accuracy of the proposed technique.

References
[1] B. Young, Digital signal integrity: modeling and simulation with interconects and packages.
Hall, 2001.
[2] C. R. Paul, Analysis of Multiconductor Transmission Lines, 2nd ed.

Prentice

Wiley, 2007.

[3] F. Caignet, S. D. Bendhia, and E. Sicard, “The challenge of signal integrity in deep-submicrometer
CMOS technology,” Proc. IEEE, vol. 89, no. 4, pp. 556–573, 2001.
[4] J. G. Hong and M. J. Lancaster, Microstrip filters for RF/microwave applications. John Wiley & Sons,
2004, vol. 167.
[5] C. Caloz and T. Itoh, Electromagnetic metamaterials: transmission line theory and microwave applications. John Wiley & Sons, 2005.
13

#1

h
#3

#2

h

#0

h

α4
α3
α2
α1

R2

R1

Figure 9: Cross-section of the curved microstrips of Sec. 6.4. Geometrical dimensions are: α1 = 10◦ , α = 20◦ ,
α = 30◦ , α = 35◦ , R1 = 80 µm, R2 = 100 µm, h = 2 µm.
4

Resistance p.u.l. [Ω/m]

10

FEM
Proposed
3

R22

10

R11
2

10

R12
6

10

7

8

10

9

10
10
Frequency [Hz]

10

10

11

10

Inductance p.u.l. [µH/m]

0.39
L11
0.38
0.37
0.36
0.35
FEM
Proposed

0.34
0.33
6
10

7

10

8

9

8

9

10
10
Frequency(Hz)

10

10

11

10

Inductance p.u.l. [µH/m]

0.75
L22
0.7

0.65
FEM
Proposed
0.6
6
10

7

10

10
10
Frequency(Hz)

10

10

11

10

Figure 10: Selected entries of the p.u.l. resistance and inductance matrices of the curved microstrip lines of
Sec. 6.4.

14

[6] A. Ametani, N. Nagaoka, Y. Baba, T. Ohno, Power System Transients: Theory and Applications. Boca
Raton, FL: CRC Press, 2013.
[7] C. R. Paul, Transmission lines in digital and analog electronic systems: signal integrity and crosstalk.
John Wiley & Sons, 2011.
[8] L. Tsang, J. A. Kong, K. Ding, Scattering of electromagnetic waves: theories and applications.
Wiley, 2000.

J.

[9] M. Kamon, M. J. Tsuk, and J. K. White, “FASTHENRY: A multipole-accelerated 3-D inductance
extraction program,” IEEE Trans. Microw. Theory Tech., vol. 42, no. 9, pp. 1750–1758, 1994.
[10] R. Harrington, Time-Harmonic Electromagnetic Fields.

McGraw-Hill, 1961.

[11] A. Ruehli, G. Antonini, and L. J. Jiang, “Skin-effect loss models for time- and frequency-domain peec
solver,” Proc. IEEE, vol. 101, no. 2, pp. 451 – 472, Feb. 2013.
[12] J. Jin, The finite element method in electromagnetics.

John Wiley & Sons, 2014.

[13] S. Cristina and M. Feliziani, “A finite element technique for multiconductor cable parameters calculation,” IEEE Trans. Magn., vol. 25, no. 4, pp. 2986–2988, 1989.
[14] G. Costache, “Finite element method applied to skin-effect problems in strip transmission lines,” IEEE
Trans. Microw. Theory Tech., vol. 35, pp. 1009–1013, 1987.
[15] E. Comellini, A. Invernizzi, G. Manzoni., “A computer program for determining electrical resistance
and reactance of any transmission line,” IEEE Trans. Power App. Syst., no. 1, pp. 308–314, 1973.
[16] K. Coperich, J. Morsey, V. Okhmatovski, A. Cangellaris, and A. Ruehli, “Systematic development
of transmission-line models for interconnects with frequency-dependent losses,” IEEE Trans. Microw.
Theory Tech., vol. 49, no. 10, pp. 1677–1685, Oct 2001.
[17] M. J. Tsuk and J. A. Kong, “Boundary integral equation formulation of skin effect problems in multiconductor transmission lines with arbitrary cross section,” IEEE Trans. Microw. Theory Tech., vol. 39,
no. 8, pp. 1336–1347, Aug 1991.
[18] M. Al-Qedra and V. Okhmatovski, “Full-periphery surface impedance for skin-effect approximation in
electric field integral equation,” Microwave and Wireless Components Letters, IEEE, vol. 19, no. 1, pp.
9–11, Jan 2009.
[19] S. V. Yuferev and N. Ida, Surface Impedance Boundary Conditions: A Comprehensive Approach. CRC
Press, 2009.
[20] Z. G. Qian, W. C. Chew, and R. Suaya, “Generalized impedance boundary condition for conductor
modeling in surface integral equation,” IEEE Trans. Microw. Theory Tech., vol. 55, no. 11, pp. 2354–
2364, Nov 2007.
[21] L. Di Rienzo, S. Yuferev, N. Ida, “Computation of the Impedance Matrix of Multiconductor Transmission Lines Using High-Order Surface Impedance Boundary Conditions,” IEEE Trans. Electromagn.
Compat., vol. 50, no. 4, pp. 974–984, Nov. 2008.
[22] K. Coperich, A. Ruehli, and A. Cangellaris, “Enhanced skin effect for partial-element equivalent-circuit
(PEEC) models,” IEEE Trans. Microw. Theory Tech., vol. 48, no. 9, pp. 1435 – 1442, Sep. 2000.
[23] A. Siripuram and V. Jandhyala, “Surface interaction matrices for boundary integral analysis of lossy
transmission lines,” IEEE Trans. Microw. Theory Tech., vol. 57, no. 2, pp. 347–382, Feb. 2009.
[24] A. Menshov and V. Okhmatovski, “New Single-Source Surface Integral Equations for Scattering on
Penetrable Cylinders and Current Flow Modeling in 2-D Conductors,” IEEE Trans. Microw. Theory
Tech., vol. 61, no. 1, pp. 341–350, Jan 2013.

15

[25] M. S. Tong, G. Z. Yin, R. P. Chen, and Y. J. Zhang, “Electromagnetic modeling of packaging structures with lossy interconnects based on two-region surface integral equations,” IEEE Transactions on
Components, Packaging and Manufacturing Technology, vol. 4, no. 12, pp. 1947–1955, Dec 2014.
[26] D. De Zutter, and L. Knockaert, “Skin Effect Modeling Based on a Differential Surface Admittance
Operator,” IEEE Trans. Microw. Theory Tech., vol. 53, no. 8, pp. 2526 – 2538, Aug. 2005.
[27] U. R. Patel, B. Gustavsen, and P. Triverio, “An Equivalent Surface Current Approach for the Computation of the Series Impedance of Power Cables with Inclusion of Skin and Proximity Effects,” IEEE
Trans. Power Del., vol. 28, pp. 2474–2482, 2013.
[28] ——, “Proximity-Aware Calculation of Cable Series Impedance for Systems of Solid and Hollow Conductors,” IEEE Trans. Power Delivery, vol. 29, no. 5, pp. 2101–2109, Oct. 2014.
[29] T. Demeester and D. De Zutter, “Construction of the Dirichlet to Neumann boundary operator for
triangles and applications in the analysis of polygonal conductors,” IEEE Trans. Microw. Theory Tech.,
vol. 58, no. 1, pp. 116–127, 2010.
[30] T. Okoshi, Planar circuits for microwaves and lightwaves.

Springer Berlin Heidelberg, 1985.

[31] L. Knockaert, D. De Zutter, G. Lippens, H. Rogier, “On the Schur complement form of the Dirichletto-Neumann operator,” Wave Motion, vol. 45, no. 3, pp. 309–324, Jan 2008.
[32] D. De Zutter, H. Rogier, L. Knockaert, and J. Sercu, “Surface current modelling of the skin effect for
on-chip interconnections,” IEEE Trans. Adv. Packag., vol. 30, no. 2, pp. 342–349, May 2007.
[33] T. Deemester and D. De Zutter, “Modeling the broadband inductive and resistive behavior of composite
conductors,” IEEE Microwave and Wireless Components Letters, vol. 18, no. 4, pp. 230–232, April 2008.
[34] T. Demeester, and D. De Zutter, “Fields at a finite conducting wedge and applications in interconnect
modeling,” IEEE Trans. Microw. Theory Tech., vol. 58, no. 8, pp. 2158–2165, Aug 2010.
[35] H. Rogier, D. De Zutter, L. Knockaert, “Two-dimensional transverse magnetic scattering using an exact
surface admittance operator,” Radio Science, vol. 42, no. 3, 2007.
[36] F. Olyslager, D. De. Zutter, and K. Blomme, “Rigorous Analysis of the Propagation Characteristics of
General Lossless and Lossy Multiconductor Transmission Lines in Multilayered Media,” IEEE Trans.
Microw. Theory Tech., vol. 41, no. 1, pp. 79–88, 1993.
[37] C. Balanis, Advanced Engineering Electromagnetics.

John Wiley & Sons, 1989.

[38] C. Balanis, Antenna Theory: Analysis and Design, 3rd ed.

Wiley, 2005.

[39] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and
Mathematical Tables. New York: Dover, 1964.
[40] Y. Ayasli, “Analysis of Wide-Band Stripline Circulators by Integral Equation Technique,” IEEE Trans.
Microw. Theory Tech., vol. 28, no. 3, pp. 200–209, Mar 1980.
[41] W. C. Chew, Waves and fields in inhomogeneous media.

IEEE Press, 1995.

[42] COMSOL Multiphysics. COMSOL, Inc. [Online]. Available: https://www.comsol.com/
[43] R. Garg, I. Bahl, and M. Bozzi, Microstrip Lines and Slotlines.

Artech House, 2013.

[44] A. Rong and Z. Sun, “A novel quasi-TEM analysis of valley microstrip lines with slit for use in multilayered MMICs,” in IEEE MTT-S International Microwave Symposium,, June 1993, pp. 955–958 vol.2.
[45] T. Hasegawa, S. Banba, and H. Ogawa, “Characteristics of valley microstrip lines for use in multilayer
MMIC’s,” IEEE Microw. Guided Wave Lett., vol. 1, no. 10, pp. 275–277, 1991.

16

[46] Y. Yin and H. W. Dommel, “Calculation of frequency-dependent impedances of underground power
cables with finite element method,” IEEE Trans. Magn., vol. 25, no. 4, pp. 3025–3027, 1989.
[47] T. Demeester and D. De Zutter, “Modeling the broadband resistive and inductive behavior of polygonal
conductors,” in International Conference on Electromagnetics in Advanced Applications (ICEAA ’09),
Sept 2009, pp. 210–213.
[48] L. Zheng and Y. Wang, “Accurate ssolution of elliptical and cylindrical striplines and microstrip lines,”
IEEE Trans. Microw. Theory Tech., vol. 34, pp. 259–265, 1986.
[49] C. H. Chan, and R. Mittra, “Analysis of a class of cylindrical multiconductor transmission lines using
an iterative approach,” IEEE Trans. Microw. Theory Tech., vol. 35, pp. 415–424, 1987.

17

