Resistance computations for multilayer packaging structures by ...

87

IEEE TRANSACTIONS ON COMPONENTS, HYBRIDS, AND MANUFACTURING TECHNOLOGY, VOL. 15, NO. 1, JANUARY 1992

Resistance Computations for Multilayer Packaging Structures by Applying: the Boundary Element Method U

Ruey-Beei Wu

Abstract-A boundary element formulation together with the moment method analysis is proposed to find the equivalent multiport dc resistances for packaging structures consisting of uniform conducting plates with arbitrarily shaped boundaries. Incorporated with a new equiangle division scheme and exact integration formulas, this approach is very efficient in the usage of computer storage and computation time such that many practical structures can be dealt with even by a small personal computer. The method is employed to investigate the effect on the equivalent resistances due to the presence of the interior apertures and the exterior boundary. Also presented are the resistance modeling technique for the multilayer packaging structures, and a novel cut and fill technique to estimate the resistances of complicated structures.

I. INTRODUCTION

T

OGETHER with the evolvement of today’s integrated circuit technology is the ever increasing device density, switching speed, and system performance. As the advent of IBM’s multilayer ceramic technology [ 11, [2], the multichip, multilayer module has become the essential trend of the present packaging technology to increase the interconnection density and decrease the transmission delay due to the chip packaging 131. Salient examples are the SX super computer packaging disclosed by NEC [4] and the AVP technology invented by AT&T [5]. Today, accurate electrical evaluation of any proposed circuit packaging scheme, at any level of packaging, is essential to the success of a given integrated circuit system. As to the electrical design and analysis of packaging structures, the signal and power distributions are the two major concerns [6]. In the power distribution design, one of the important issues is the dc voltage drop due to the packaging interconnections between the IC chips and the power supply. This issue becomes more and more critical, since the voltage immunity becomes lower due to the reduced device voltage level and the resistance per unit length becomes larger due to the miniaturized packaging wiring and redistribution plates. For multichip multilayer modules, the complexity of geometries of the conducting plates and the existence of numerous contact points or ports representing different boundary conManuscript received February 14, 1991; revised July 30, 1991. This work was supported in part by the National Science Council, Republic of China, under Grant NSC 79-0404-E002-017. The author is with the Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan. IEEE Log Number 9104489.

ditions make an analytic solution all but impossible. The simplifying assumptions that have to be made quite often result in a model which bears little resemblance to teality. For this sake, Sakkas applied the finite element method (FEM) to solve the potential distribution for arbitrary geometry and produced an equivalent resistance network as a model of the geometry [7]. Given the network for every conducting plate, arbitrary three-dimensional multilayer packaging structures can then be analyzed by general circuit analysis programs such as ASTAP [8] or SPICE. FEM calls for a mesh to discretize the conducting region of the plate into a lot of small elements, inside which basis functions with unknown coefficients are employed to approximate the desired potential distribution [7]. The method is well suited for arbitrary geometries and regions of nonuniform electrical properties. However, the discretization of the continuum requires a lot of unknowns, especially in the region near the corners of the boundaries where the potential varies rapidly. Also, a large fraction of the unknowns are usually spent in modeling the potential near the plate boundary which is of large extent, but has relatively small effects on the equivalent resistances between the ports in the center region of plate. As a result, the solution of the potential distribution for common structures usually requires larger computer storage or computational load than a small computer like IBM/PC can afford. In practical structures, the conducting plate is usually of uniform electrical properties but the shape and position of the ports may be rather arbitrary. The boundary element method (BEM), which requires the solution for the potential distribution along the boundaries only [9], is proposed in this paper to deal with the resistance computations. Usually, the discretization is much easier and the number of equations is one dimension lower. Hence, BEM can be numerically superior to FEM such that the resistance computations for practical problems can be executed on a small computer. This makes the method rather attractive for industry applications. 11. GOVERNING EQUATIONS FOR POTENTIAL DISTRIBUTION Fig. 1 shows a thin conducting planar structure which is uniform with thickness t and conductivity U . The plate is bounded by an exterior boundary and may include several interior apertures; both are represented by the boundaries r N . Inside the plate, there are several (say M ) ports serving as

0148-6411/92/$03.00 0 1992 IEEE

Authorized licensed use limited to: National Taiwan University. Downloaded on February 25, 2009 at 00:39 from IEEE Xplore. Restrictions apply.

IEEE TRANSACTIONS ON COMPONENTS, HYBRIDS, AND MANUFACTURING TECHNOLOGY, VOL. 15, NO. 1, FEBRUARY 1992

88

Y

ith port is given by evaluating the line integral of the current density around the port boundary, i.e.;

Exter 101 ,/boundary

A

I, =

f s ( z , y ) . ii dC =

-00

TD,

Fig. 1. A typical thin conducting plate with ports, apertures, and exterior boundary.

1

-

(2)

The boundary conditions to (1) are such that the potential should equal the impressed voltage, i.e.:

4 ( x , v ) = V,

(5)

TD,

#j

and

1 1 ~ ~(7)~ .

gii = j#i

Since the conducting plate always dissipates power, the matrix [GI is symmetric and semipositive definite. The resistance problem in a planar structure is in close analogy to the capacitance problem in a transmission line [lo]. The potential distributions in both problems satisfy the same Laplace equation (1). The current density ii . J(xI y) outward to the boundaries in the resistance problem corresponds to the surface charge density p ( z , y) in the capacitance problem. Hence, the two problems become mathematically dual with the dual quantities defined in the following:

(1)

in the conducting region A . Since the electric field is defined by the gradient of the potential distribution, the current distribution in the conducting region A is given by

Y) = -ffoV4(x, Y).

dC.

The equivalent resistances can be found from the matrix elements of [GI by the relations gzj = - I / R ~ . ~ for , i

input or output terminals whose boundaries are denoted by ro. In the common cases where the distance between any two boundaries is much larger than the thickness t, the current flow can be assumed to be planar. Then, the conducting property of the plate can be denoted by the sheet conductance gm which equals o t at dc or very low frequencies. Assume that each of the M ports has a terminal voltage imposed, say V , at the ith port (i = 1, 2 , . . . , M ) . Current will flow in the conducting plate due to the voltage drops between the terminals. In other words, there will be a potential distribution $(IC, y) which satisfies the Laplace equation

4.1

2

It is not difficult to verify from (1) the neutrality condition that the sum of total current over all the ports is zero. The relationship between the impressed voltages and port currents defines the equivalent resistance R,, = R,, between any two ports i and j in the conducting plate. More elegantly, let [VI and [ I ]be two 1 x M column vectors denoting the impressed voltages and the outward currents, respectively, at the M ports. These two vectors can be related by an M x M conductance matrix [GI, i.e.;

\ Interror

V24(X,Y) = 0

f

Resistance in a plate

Capacitance in a transmission line

potential ~ ( ry ) . potential d(s,y ) current density zi J(s.y) charge density d x ,Y) sheet conductivity UD permittivity e conductance matrix [GI shortcircuit capacitance matrix [C] equivalent conductance ~ / R L J two-terminal capacitance -CZJ( # 3 )

(3)

However, conventional capacitance analysis cannot be applied directly to the present resistance problem. In capacitance at the ithe port boundary ro, for i = 1, 2 , . . . , M , while problems, the potential field usually extends to the whole space along the exterior boundary or the boundaries formed by the and faces no exterior boundary. As a result, the equivalent apertures, the normal component of the current density should capacitance will be affected by the presence of a substrate be zero, i.e., ii . J ( x ,y) = 0, or or inhomogeneous material, if any. Different kinds of material inhomogeneity may warrant different analysis procedures [ 101. (4) In the resistance problems, the potential distribution of interest is confined inside the conducting plate which is bounded by at r N , where n is the unit vector outward normal to the an exterior boundary and may include some interior apertures. boundary. Usually, the plate is made of a particular metal with the same Associated with the Dirichlet boundary condition (3) and thickness such that the material can be assumed homogeneous. the Neumann boundary condition (4), a unique potential However, the presence of the exterior boundary and interior distribution $(x,y) can be solved from the partial differential apertures will increase the equivalent resistance and must be equation (PDE) (1).Then, the total current flowing out of the taken into account.


89

WU: RESISTANCE COMPUTATIONS FOR MULTILAYER PACKAGING STRUCTURES

111. BOUNDARY ELEMENTFORMULATION (BEM)

BEM, or boundary integral equation method, is a powerful approach to solve a PDE. It can reduce the two-dimensional PDE into one-dimensional integral equations along the boundaries [ l l ] . The number of unknowns is much smaller than its counterpart in FEM such that the numerical solution can be more efficient. Consider the geometry in Fig. 1 where the potential distribution satisfies the PDE (1). Without loss of generality, the sheet conductance is assumed to be unity in the following analysis. By applying BEM, we select suitable weighting functions w(T) = w(2, y) and require the residual of (1) to be orthogonal to them, i.e.:

By imposing the boundary conditions (3) and (4), (11) gives the desired surface integral equation for the variable 4 and its derivative 84 along the boundaries. Since the potential 4 is a given constant along the boundary of each port, (1 1) can be reduced to

v,=

f G(T;P’) . p ( ~d!) + f dG(r; T ’ ) . 4 ( ~d!) r D

(12)

rN

for T’ lying along the port boundary r D , (i = 1, 2 , . . . ,M ) . Here, we have applied the identity that

P

1

dG(r; T’) de

TD,

-1, -112, 0,

if 5;’ inside the region bounded by r D , if T‘ on the boundary r D , if T‘ outside the region bounded by r D , .

After integration by parts twice, (8) becomes

(13)

On the other hand, for T‘ along the exterior boundary or the interior aperture boundaries r N , the integral equation reads where dw is a short hand notation for the normal derivative d w l d n = ii . Vw and similarly for d4. The boundary r = U {r D , r N } is the union of the exterior boundary and the boundaries of all the terminal ports and the interior apertures. It should be noted that n is defined as the outward normal to the conducting region A. Along the boundaries of the ports, it actually points inward to the patch of the port terminal. Therefore, we will finally replace -84 along the port boundary by its dual function p(x, y) for the sake of convenience. The area integral in (9) can be greatly simplified if the weighting function is chosen to be the Green’s function. Then, we obtain the formula

f (4. dG - G

d! =

{

t(T’),

r

for T’ inside A for T‘ outside A (10)

1 2T

G(T; T ’ ) = - lnlF - 5;’1 and 1 cos(?& 5; - T’). 2TlT - T’l

It is noted that the potential distribution 4(x,y) at any point inside the plate boundary can be calculated from $I and its normal derivative 84 along the boundaries r by (10). When the point T’ lies on r, the singularity of function dG(T; T ’ ) at T = 5;’ makes the first term in the left-hand side (LHS) of (10) not integrable. Nonetheless, we can take the average of the two equations in (10) and obtain the formula

P

[

4.dGdl-

f

1 2

G . d $ d e = -4(T’)

f

p(T) d l = 0.

(15)

TD

IV. MOMENTMETHODSOLUTION The method of moments [13] is applied to solve the integral equations (12) and (14). In the discretization, the boundary of the nth port, r D , , is partitioned into Nn segments (n = 1, 2 , . . . ,M ) , while the exterior and aperture boundaries are partitioned into N M + segments. ~ Over each segment, the basis function is chosen to be the pulse function. Along the exterior or aperture boundary, the unknowns are the average potential over each segment

where

sdG(T; T’)=

In the case where the exterior boundary of the conducting plate is too far to affect the equivalent resistance, the problem can be modeled as of an open boundary type. Equations (12) and (14) should be modified by adding to the right-hand side (RHS) the potential at infinity qhm. The additional unknown dm is counteracted by the additional condition that the sum of all the port currents is zero, or

(11)

where r N , is the kth segment of r N and is its arc length. Along the port boundary, it is beneficial to choose the unknown along the segment r D n , , the jth segment of the boundary of the nth port, as qn,

=

J p(~)de

(17)

rDn,

where P denotes the Cauchy principal value [12].

which represents the total current flowing out of that segment.


~

90


be of arbitrary shapes, but are modeled as polygons. In the numerical computation, the first step is to divide all the polygon sides into the union of small line segments. The number of line segments decides the number of unknowns. Naturally, a finer division with more line segments can achieve better accuracy, but inevitably increases the required computer storage and the computation time. The conventional division for i = 1, 2 , . . . , N , and m = 1, 2 , .. . ,M . The matrix scheme based on the equilength criterion is not efficient for elements are defined by this problem, since the exterior boundary with large periphery may exhaust the available computer storage in many cases. Actually, the potential field satisfying Laplace equation is invariant under the conformal mapping [ 14, p. 2771. This implies that the suspension angle rather than the length of and the line segment is the critical factor in deciding the accuracy of the results. Hence, the program incorporates an automatic division scheme based on the equiangle criterion. For each polygon side of a certain boundary, the program first calculates its distance to each side of all the other boundaries. Then, the A similar procedure with (14) leads to another set of equanumber of divisions is chosen proportional to the ratio of the tions maximal and minimal distances. In this way, the far away M N, NM+I boundary which has a small influence on the field distribution 1 &,n, . qn, + ( ~ ~ -, Thplc) k . 4 k = 0 (21) is reasonably modeled with a few segments. Results of good n = l j=1 k=l accuracy can be achieved with a few unknowns by this for p = 1, 2 , . . . , N M + ~Here, . the definition of Pp,n, and improved division scheme. The second step is to calculate the matrix elements by Qp,k are similar to those in (19) and (20) except that Q p . k is zero if k = p . Also, the term h,, is unity when k = p and evaluating (19) and (20). Since the involved operation is a double integral, it may become quite time consuming if is zero if IC # p . Equations (18) and (21) form a system of equations for numerical methods such as Gaussian quadrature or Simpson the unknowns qn, and 4 k where j = 1, 2,...,Nn for integration are used. Here, we derive closed form formulas n = 1, 2 , . . . , M and k = 1 , 2 , . . . ,Nhf+l. After speci- for (19) and (20) in Appendix A. Though the expressions fying a set of suitable voltages Vm(m= 1, 2 , . . ., M ) , the in (A.7) and (A.ll) are exact, the argument of the complex boundary unknowns can be computed. Thus the total cur- logarithm function, In ( z - z’), should not be across the branch rent I n ( n = l , 2, . . . , M ) flowing out of each port is finally cut, i.e., the negative real axis, to circumvent the multiplicity obtained. In case the exterior boundary is neglected, the of In ( 2 - 2 ’ ) . A simple way is to replace ln(z - z’) with additional unknown qtm should be included on the LHS of (18) ln(z’ - 2 ) in (A.7) and (A.ll) when the center point of is and (21). Then, an additional equation due to the neutrality to the left of that of I”. By specifying the impressed port voltages as the boundary condition (15) is required, i.e.: conditions, we can solve the matrix equation and obtain the M hl N . outward port currents. To find all the equivalent resistances, the matrix equation must be solved several times correspondn=l n = l j=1 ing to different boundary conditions [7]. Here, the matrix is LU-factored and the multiple solutions are obtained by substitutions [15]. The direct method is advantageous here V. SOLUTION PROCEDURES since the involved factorization, which represents most of A number of electrically insulated plates, connected to one the work required to obtain a solution, is done only once. another through vias, can represent many present day packag- However, the matrix will lose its significance if the rows ing designs in IC technology. To analyze the dc properties of corresponding to the exterior boundary are factored before these multilayer structures, it is crucial to find the equivalent those corresponding to port boundaries. Hence, in the program, resistance network which is made of the via resistances and the unknowns are numbered first for the ports, then for the the equivalent resistances between each two ports in the same interior apertures, and last for the exterior boundary. It has plate [7]. The calculation for via resistance is obvious, given been found that the direct method is very efficient and always the cross section area and the length. However, some suitable gives satisfactory results. computer-aided analysis is necessary to obtain an accurate Merging the equivalent plate resistance models with the resistance model of the plate. via resistances can give an equivalent resistance network. Based on the aforementioned method, a Fortran program Since a lot of ports are internal to the packaging structure, R3D has been tailored to calculate the equivalent resistances it is advantageous to extract a macromodel consisting of the between multiple ports in a conducting plate of uniform equivalent resistances between the external input/output (I/O) conductivity. The ports, apertures, and exterior boundary may ports only. By applying the general network theorem, the Using Galerkin’s method [13], (12) is weighted by pulse functions such that its average for F’ lying on each segment I‘D,, is satisfied. Hence, we obtain a set of equations

C

r


91


macromodel can be obtained from the equivalent resistance network. The algorithm is described briefly as follows. All the ports which can be connected together by some physical paths construct a group. The original resistance network may include several groups, but the ports of the same group need to be considered at once. Let the voltages and currents of those ports be denoted by column vectors [V,] and [I,], respectively. They will satisfy a linear relation

[GI. [ K ]= [In]

Req 1.5

1.0

(23)

where the conductance matrix [G,] can be found from the resistance network easily. If there is a branch with conductance gij between ports i and j , [G,] is formed by adding gij to the elements G,(i,i) and G,(j,j) while adding -gij to the elements G,(i,j) and G,(j, i ) . However, the matrix formed by this way is singular. To solve the difficulty, we may choose one of the 1 / 0 ports in the group as the voltage reference. The diagonal term corresponding to the reference port is replaced by a very large number, say lolo, such that the matrix becomes nonsingular and the voltage at the reference port is properly confined to be almost zero. It is noted that the matrix [G,] is symmetric and sparse. Here, we employ a skyline representation [15] for the matrix to save the required computer memory and computation time in the following matrix solution. As far as the macromodel for the 1/0 ports is concerned, (23) has the boundary conditions that the net current of each interior port is zero. We may thus impose in the RHS of (23) that the current of theith 1 / 0 port be 1, the current of the reference port be -1, while all the others are zero. After solving the matrix equation (23) for the given RHS vector, the voltages at all the 1/0 ports can be extracted. Actually, if the relation between the voltages [V,] and currents [I0]of the 1 / 0 ports other than the reference one is written by

[VOI = P o l . [I01

2.0

(24)

then these extracted 1/0 port voltages just form the ith column of the matrix [R,].If the boundary conditions are specified for each i sequentially, the matrix [R,] in (24) can be obtained accordingly. Again, the matrix solutions involved here are done by the LU-factorization, but with skyline storage form [15] to take advantage of the sparsity of [G,]. The inverse of [R,], called [Go], gives the desired equivalent resistance network of the 1 / 0 ports. The negative of the off diagonal matrix element, say -Go(i, j ) , represents the equivalent conductance between the ith and jth 1/0 ports. The diagonal matrix element G,(i, i) represents the sum of the equivalent conductances between the ith 1 / 0 port and all other 1 / 0 ports including the reference port. Hence, the equivalent conductance between the ith 1/0 port and the reference port can be obtained from the sum of the matrix elements in [Go] over the ith row. Finally, the reciprocals of these equivalent conductances give the desired equivalent resistance network for the macromodel. VI. NUMERICAL RESULTS

Some examples are given in this section to investigate the finite boundary effects on equivalent resistances and to

0.5

0.0

f

I

I

1

10

100

d/r Fig. 2. Equivalent resistance between two circular ports in a circular plate.

discuss some resistance modeling techniques. For simplicity, all the conducting plates considered here are made of unit sheet conductance, if not specified explicitly. The equivalent resistances presented here are then the normalized results with respect to the sheet resistance. Let us first consider the equivalent resistance between two circular ports in a circular conducting plate with radius w (see the structure in Fig. 2). The two ports are of radius r and situated at a center-to-center distance of 2d. When w is very large, the exterior boundary has negligible effect. By applying suitable conformal mapping [14, p. 3291, the outward current density along the boundary of the left port under unit voltage drop can be analytically written as [16]

Jr(6') =

on

___

cosh-la

.

-

2r . ( a - cos 6') '

where 6' is the angle measured from the exact equivalent resistance is

a =d/r 2

(25)

axis. Also, the

It is noted that the equivalent resistance depends on the ratio a = d / r . That is, the normalized planar resistance remains invariant under scaling. The calculated resistance versus the ratio d / r is shown as the solid line in Fig. 2. The equivalent resistance equals the sheet resistance when d is about ten times of r. A conventional intuitive approach estimates the equivalent resistance by counting the number of squares between the two ports. This example shows that the intuitive approach can hardly be applied here. Also, the analytic results by (26) are shown by the triangular marks in the same figure. It is found that the calculated results are satisfactory over the very wide range of d / r from 1.1 to 100. By imposing a unit voltage on the left port while the right port is grounded, we calculate the outward current density along the boundary of the left port. The solid lines in Fig. 3


92


without boundary ----with boundary 1.5

10

05

,

0

,

,

.

.

. . . . .

. . . . .

,

610'

1

ioo

0.0 1 boo

e Fig. 3.

!

I

1

10

1

100

d/ r Fig. 4. Equivalent resistance between two square ports in a square plate.

Current distribution along one of the two circular ports in a circular plate.

show the numerical results with r = 1 and the ratio d / r = 1.1, 2, and 10 as the parameter. The curves are in excellent agreement with the analytic solutions given by (25), though not plotted here. For far away ports, or a large d / r ratio, say d / r = 10, the current density is nearly uniform along the port boundary. The current flow in the region distant from both ports is not negligible. As the d / r ratio decreases, the current density increases significantly at the side near the other port, but slightly at the far side. The current flow actually concentrates on the region between two ports. When the radius of the plate is not very large, the effect due to finite exterior boundary should be taken into account. As an example, we assume w = 2d to calculate the equivalent resistance and current density. The numerical results are shown by the dashed lines in Figs. 2 and 3, respectively. It is found that the presence of the exterior boundary makes the equivalent resistance higher since it will confine the current flow in the distant region. While the current density in the near side is hardly affected by the exterior boundary for all d / r ratios, the current density in the far side changes significantly. Especially for small d / r , e.g., d / r = 1.1, it decreases rapidly and becomes negligible for 0 larger than 150'. Fig. 4 shows the calculated resistances for the same structure except that the shape of the ports and the exterior boundary are square. Now, the parameters r and w denote the half side lengths of the ports and the exterior boundary, respectively. Though no analytic formula such as (26) exists for square ports, the equivalent resistance in Fig. 4 shows a similar behavior to that in Fig. 2. It is noted that its value is smaller at the same d / r ratio since the equivalent path between two ports is smaller for square ports than for circular ports. Now, we apply the R3D program to consider the resistance modeling for the mesh plane in the MLC structures [l]. Fig. 5(a) shows a typical mesh plane which basically is periodically perforated with holes of size a and pitch b. To connect to other planes through vias, some of the holes are

patched and they serve as the I/O ports. It is important to evaluate the equivalent resistance between any two ports in the mesh plane environment. However, thousands of apertures in the mesh plane make the numerical computation all but impossible. Instead, we may employ a suitable modeling technique to find a good approximation. Consider, for instance, the equivalent resistance between the two cross hatched ports in Fig. 5(a). We may cut the mesh plane at some boundary to form a new structure as shown in Fig. 5(b) which includes 3 . 5 = 15 holes only and is much easier for a numerical solution by R3D. The current flow in the distant region, albeit small, is artificially neglected such that the equivalent resistance calculated from Fig. 5(b) is an upper bound of the exact solution. On the other hand, we may fill the far away apertures and form another structure easily for the numerical solution as shown in Fig. 5(c). This structure offers some extra current flow paths and results in a lower bound of the exact equivalent resistance. The two bounds will become closer to each other if more holes are included by the present cut and fill technique. The exact resistance can thus be determined with good accuracy, if the computational labor can be afforded. Let R,, denote the resistance between two ports which are m and n holes apart in the x and y directions, respectively. Fig. 6 shows the resultant upper and lower bounds of the equivalent resistances R I OR11, , R20, and R21. Here, the lower bounds are obtained by considering open structures with 6 . 7, 6 . 8 , 7 . 7, and 7 . 8 holes, respectively. Also, the upper bounds are obtained by including 4 . 5. 4 . 6, 5 . 5, and 5 ' 6 holes together with the boundary of cut. When the hole size is small or when the two ports of interest are not close, a large extent of the current flow is in the distant region. The present cut and filZ model, which either neglects or overcounts the distant region, may result in a large discrepancy between the two bounds. Also, it is found that the adjacent resistance RIObecomes smaller as the hole size increases. However, all the other R,,'s become larger since the adjacent apertures



0

93

analysis of multilayer packaging systems, we consider an imaginary structure of three power distribution layers as shown in Fig. 7(a). The voltage supply comes in from the bar, i.e., node 10, in the middle plate, through the vias, upper plate, and pins to the chip pads, i.e., nodes 31 and 33. The return current flows from chip pad 32, passing through two plates to the lower plane, and to the ground bar, i.e., node 1. The dimensions of the plate thickness and via diameter are typical for the present MLC technology. For simplicity, the plate is assumed to be made of copper with conductivity 5.7 x lo7 mho/m. After executing the program R3D, we can obtain the equivalent resistance network as shown in Fig. 7(b). The via resistances are Ro and R I , while the resistance models for upper, middle, and lower plates are represented by the matrices [Ru],[ R M ]and , [ R L ]respectively. , The network is a detailed micromodel which includes a lot of nodes internal to the structure only. As far as the external circuits are concerned, only the equivalent resistances between 1/0 ports (nodes 1, 10, 31,32, and 33) are relevant. In this respect, we may employ the general network theorem to obtain the macromodel as shown in Fig. 7(c). Given the resistance model, we can investigate the voltage drops at the chip pads due to a current flow, say 10 A. If the chip pad with node 31 is ON while the node 33 is left floating, the incoming voltage will suffer from a voltage drop of - 1 0 . [4.811//(4.762 16.72)] = -39.31 mV at node 31 and -39.31 . 4.762/(4.762 16.72) = -8.71 mV at node 33 with respect to the power supply. Meanwhile, the ground pad 32 will present a residual voltage of +80.24 mV with respect to the package ground. In other words, the voltage supply to the chip circuits is decreased by 39.31 80.24 120 mV. Similarly, the voltage drops will be -39.00, -8.71, and +80.24 mV at the chip pads 33, 31, and 32, respectively, when the node 33 is ON and the node is 31 is OFF.

+

0 1 0 O D O (c) Fig. 5. Resistance modeling for a mesh plane by the cut and $11 technique. (a) Original structure. (b) Upperbound modeling. (c) Lowerbound modeling.

Re,

4-

- without

--- w i t h

boundary

VII.

2-

0.1 -

RlO l

06

'

'

+

boundary

2-

0.5

+

'

'

l

0.7

'

'

"

l

~

'

'

l

i

0.9

0.8

a/b Fig. 6. Equivalent resistances, R I O .R11, R20. and R21 versus the a / b ratio. The dashed and solid curves depict the resultant upper and lower bounds of the exact solution, respectively.

are now large enough to block most of the current flow. To illustrate the applications of R3D for the dc electrical

CONCLUSIONS AND

DISCUSSION

A computer-aided analysis tool R3D has been established to compute the equivalent resistances for packaging structures with multiple conducting plates. Based on the BEM together with the method of moments, this approach calls for the unknown current density along the port boundaries and the potential distribution along the interior apertures and exterior boundary. The equivalent resistance model for each plate can be found from these unknowns. By assembling the equivalent models for all the conducting plates with the via resistances, the equivalent resistance network of a multilayer packaging structure can be established. As far as the external circuits are concerned, the macromodel consisting only of the equivalent resistances between the I/O ports can be obtained by applying the general network theorem. Given the macromodel, the dc electrical properties of the power distribution in the packaging structures can be investigated. In the numerical solution, the program incorporates two novel features. One is the automatic division scheme based on an equiangle criterion, the other is the closed form integral formulas for the evaluation of matrix elements associated with


IEEE TRANSACIlONS ON COMPONENTS, HYBRIDS, AND MANUFACTURING TECHNOLOGY, VOL. 15, NO. 1, FEBRUARY 1992

94

Unit = m f l

hole rodbus=& via radius= 2 plate = 2 O O x 200x1

31

TA

32

n

unil

E

mil

3:

33

I

I

32

33

n

l -

node

22

21

2 343

22

23

3 770 6 463

23

24

25

1 355

3 770

5 947

1 355

2 343

1 355

5 947

24

5 947

25

5 946

1

355

9 G73

1

[RMI

26

6 462

nodel

::

12

13

14

10

2 157

1878

6688

2 246

6 693

8 751

5371

2 757

2 343

13 Power s u p p l y

5374

-31

32

33

[U,]

[Rpl

--

10

I

P P I

q 5 33

1 762

u' G1*

(4 Fig. 7. An example to illustrate resistance modeling for a multilayer structure. (a) Original structure. (b) Equivalent network. (c) Macromodel.

the Laplace problems. As a result, the approach is very efficient in the usage of computer storage and computation time such that a lot of practical structures can be handled by a small computer like IBM/PC-AT within a few minutes. For example, the numerical analysis for the structure shown in Fig. 7(a) calls for about 80 unknowns, or 80 .80 = 6.4K computer storage. The required runtime is about 40 s by a 28-MHz 80386-based PC coped with a 80387 coprocessor. Several numerical examples are included to demonstrate the versatility of the proposed method. The effects on the equivalent resistance and the current distribution due to the presence of the exterior boundary are discussed. For more complicated structures such as a mesh plane in the MLC structure, we propose a novel cut and $11 technique to find upper and lower bounds for the desired resistance. This technique, if properly applied, can give satisfactory results without too much computational effort.

INTEGRATION

APPENDIXA FORMULASFOR LAPLACE PROBLEMS

This appendix presents the analytic formulas for the integrations in (19) and (20). As shown in Fig. 8, we first consider the integral

where F = (5, y) and F' = (XI, y') denote the coordinates in the two line segments r and r' which are of lengths e and e', respectively. By the complex variable technique, we may define the dummy complex variables z and z' corresponding to F and F', i.e.:

z = x + iy;

z' = x'


+ iy'

(A4

95


The integration in (20) considers

In terms of complex variables, it is not difficult to verify that cos(ii, r - r’)

IV - V’l

1 de = Irn{ -d z } . 2

- z‘

(A.9)

Therefore, it becomes obvious that the inner integral in (A.8) gives

J r

Fig. 8. Integration along two straight line segments to evaluate matrix elements.

respectively. Then, the integrand in (A.l) is just the real part of an analytic complex function In ( z - 2 ’ ) . It is thus reasonable to consider the complex integral 2:

10

=

If -

-

?=’I

TI)

ln(z - z’)lE:}.

(A.lO)

Note that the RHS in (A.10) is just the suspension angle between the two end points z1 and z2 with respect to the point z’. It is independent of the integration path such that (A.lO) remains valid even when I? is not a straight element. Since the term in the bracket of (A.lO) is an analytic function, the same procedure in deriving (A.7) can be applied here. As the result, we obtain

z2

J J ln(z - z‘)dz dz’ 2;

cos(ii, r

(A.3)

21

where z1 and z2 correspond to the coordinates of the two end points of I?, and similarly for z i , zh, and I?’. By an analytic integration, the closed form solution to (A.3) is

Due to the presence of (2; - 2 ; ) in the denominator, ( A . l l ) is exact only when the integration path I” in (A.8) is a straight segment.

REFERENCES On the other hand, the dummy complex variables can be written in terms of two real dummy variables s and s’ by

z = z1

+

(22 -

z1) . s / e and

z’ = 2:

+ (2; -

2:)

. s’/!’ (-4.5)

respectively. Then, the integral in (A.3) becomes

I o = (z2 - 2 1 ) . ( 2 ;

- 2:)

. 0

0

(A4

Equating (A.6)-(A.4), the integral in (A.l) which is just the real part of the bracket in (A.6) is given by

111 B. T. Clark and Y. M. Hill, “IBM multichip multilayer ceramic modules for LSI chips-Design for performance and density,” IEEE Trans. Comp., Hybrids, Manuf: Technol., vol. CHMT-3, pp. 89-93, Mar. 1980. [2] A. J. Blodgett, Jr., “Microelectronic packaging,” Scientif: Amer., vol. 249, pp. 86-96, July 1983. [3] T. Dixon, “Multilayer ceramics: The key to high density interconnections,’’ Electron. Pack. Production, pp. 76-82, Feb. 1983. [4] T. Watari and H. Murano, “Packaging technology for the NEC SX supercomputer,” IEEE Trans. Comp., Hybrids, Manu$ Technol., vol. CHMT8, pp. 462-467, Dec. 1985. [5] C. J. Bartlett, J. M. Segelken, and N. A. Teneketges, “Multichip packaging design for LSI-based system,” IEEE Trans. Comp., Hybrids, Manuf: Technol., vol. CHMT-10, pp. 647-653, Dec. 1987. [6] E. E. Davidson and G. A. Katopis, “Package electrical design,” in Microelectron. Packaging Handbook, R. R. Tummala and E. J. Rymaszewski ed., New York Van Nostrand Reinhold, 1989, ch. 3, pp. 111-166. [7] T. M. Sakkas, “Potential distribution and multi-terminal dc resistance computations for LSI technology,” IBM J. Res. Develop., vol. 23, pp. 640-651, NOV. 1979. [8] W. T. Weeks, A. J. Jimenez, G. W. Mahoney, D. Mehta, H. Qassemzadeh, and T. R. Scott, “Algorithms for ASTAP-A network-analysis program,” IEEE Trans. Circuit Theory, vol. CT-20, pp. 628-634, Nov. 1973. 191 C. A. Brebbia, The Boundary Element Method for Engineers, London, U.K.: Pentech, 1978, p. 181. [lo] W. T. Weeks, “Calculation of coefficients of capacitance of multiconductor transmission lines in the presence of a dielectric interface,” IEEE Trans. Microwave Theory Tech., vol. MlT-18, pp. 35-43, Jan. 1970. [ l l ] L. Lapidus and G. F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering. New Jersey: Academic, 1982, sec. 5.5, pp. 461-481. [12] J. Mathews and R. L. Walker, Mathematical Methods of Physics. New York: W. A. Benjamin (reprint), 1970, second ed., p. 480.


96


[13] R. F. Harrington, Field Computation by Moment Methods. New York: MacMillan, 1968. [14] R. V. Churchill and J. W. Brown, Complex Variables and Applications. New York: McGraw-Hill, 1984. [15] K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. Englewood Cliffs, NJ: Prentice Hall, 1976, ch. 7, pp. 239-258. [16] Y.-C. Lin and R.-B. Wu, “A hybrid approach for two dimensional electrostatic problems,” Bull. College Eng., National Taiwan Univ., no. 47, pp. 153-162, Oct. 1989.

Ruey-Beei Wu was born in Tainan, Taiwan, Republic of China, on October 27, 1957. He received the B.S.E.E. and Ph.D. degrees from National Taiwan University, Taipei, Taiwan, in 1979, and 1985, respectively. In 1982, he joined the faculty of the Department of Electrical Engineering, National Taiwan University, where he is now a Professor. In 1986, he was a Visiting Scientist for one year with the IBM General Technology Division Laboratory, East Fishkill Facility, Hopewell Junction, NY. He currently works primarily on the applications of numerical methods to electromagnetic field problems He has been engaged in research on dielectric waveguides, optical fibers, wave scattering of anisotropic objects, slot antennas, microstrip discontinuities, and interconnection modeling for computer packaging

-_
