In this study, a multiscale interface element method (MIEM) is developed to evaluate the interfacial peel and shear stress distributions in multilayer packaging structures. A newly developed polygon interface element contained many sides that could easily connect with many other polygon elements of much smaller scale. Simple model and less computation can obtain accurate stress distributions. After verification of the validity of the developed model, MIEM is applied to analyze thermal stress distribution of the chip scale package (CSP). In general, a good qualitative agreement has been observed between the various results. Moreover, MIEM can directly obtain the stress value at the interface, which is also one of the advantages of this method in dealing with multilayer structures.

1. Introduction

Higher integration has become a trend in modern electronic packaging. One of the major challenges is to improve the thermomechanical reliability when the package is subjected to the thermal loads. However, the conjunct interfaces near the free edge always suffer high stress gradients and even can generate cracks, because the packaging components are fabricated with different thermal and mechanical properties [1, 2]. Suhir [3] predicted the thermal stresses in the layered structures. The method provided an easy-to-use estimation method for interfacial stresses between layers. Pao et al. [4, 5] developed Suhir’s model for extension to multilayer structures. The proposed model can estimate the total force produced by the solder joint from the equation of force balance. Oda and Sakamoto [6] proposed a special finite-element method for multiple laminated structures. From the results of the basic IC model, the influences of varying chip size and package thickness for the stress distributions are presented. Ye et al. [7] analyzed a simple insulated gate bipolar transistor (IGBT) model and obtained the stress and strain in the solder layer. From the results, the thermal cycling had significant effect on the solder layer reliability. Ghorbani and Spelt [8, 9] developed a two-dimensional analytical model to investigate the interfacial thermal stresses in the solder joints of a leadless chip resistor (LCR). And the results were compared with finite element analysis.

Finite element method (FEM) has been widely used to solve engineering problems in recent decades. But for displacement discontinuity problems, the stability of numerical results will not be ensured. Zhuo et al. [10] presented the interface stress element method (ISEM) for solving the problem of discontinuous media mechanics. The idea came from rigid body spring model (RBSM) presented by Kawai [11]. Zhuo and Zhang [12] applied ISEM to analyze engineering problems; the ISEM is more accurate and applicable compared with the traditional FEM. Zhang et al. [13] proposed an ISEM-FEM-IEM hybrid model, which solved the conflict between computational accuracy and efficiency. These examples demonstrate the applicability and adaptability of the model to the engineering. In addition, Zhang et al. [14] performed algorithm of ISEM for block elements of arbitrary shape using natural coordinate. The numerical results show that the method was feasible and accurate.

As mentioned above, many studies have shown the use of ISEM in analyzing the engineering problems. In this study, ISEM is applied to analyze the thermal stress distribution of multilayer structures. Not only that, a new polygon interface element is developed, called multiscale interface element method (MIEM), which could easily connect with many other polygon elements of much smaller scale. Finally, the developed model was applied to analyze thermal stress distribution of the chip scale package (CSP). In this case, MIEM display the accuracy of analyzing structural singularities.

2. Derivations of Interface-Element Method

2.1. The Governing Equation

The interface element method (IEM) is based upon the discipline of stress vector (tractions) continuity on the interfaces by assuming rigid element and satisfying interfacial constraint equations. In IEM model, a complicated domain can be subdivided into a series of smaller regions named block element. The block element is undeformable and it can be separated and produces rigid body motion when an external load is applied. The deformations are accumulated on the conjunct interface of block elements. The internal energy caused by deformation of the interface area can be related to the work done by the surface traction around the block element. For equilibrium solution, it satisfies the fact thatThe discretization form can be written as (2), where , , , and denote the element domain, stress boundary, interface region, and interfacial stress.

2.2. Displacement Distribution in Interfacial Element

A four-node 2D interface element for illustrated purpose is shown in Figure 1, in which and represent the local coordinate system that are normal direction and shear direction to the interface surface. The displacement field of the block element could be expressed by the rigid body motion, translational and rotational, such as where denotes the displacement vector at any point in the block element, is the degree of freedom at the mass of center, and is the shape function. With respect to the local system, the relative displacements at interface element are then derived as In the above equation, is normal displacement, is the shear displacement, and L is the coordinate transformation matrix which reads

2.3. Stress Distribution in Interfacial Element

A measure of the relative displacement of the block element is referred to as strain on interface. As shown in Figure 2, for any point M on interface surface, the normal strain contributed from block element 1 is expressed as the change in length per unit of the original length in the local system, where is the normal distance from the centroid of block element 1 to point M. The shear strain is evaluated from the ratio of slide distance from block element 1 to the original length . The strain in the region of block element 2 is calculated in a similar way. In summary, the components of strain are defined as follows:In plane stress situation, the components of stress can be evaluated by the following equation.where and are the normal stress and shear stress contributed by block element i. The interface stress should be balanced due to the system being in an equilibrium state. It satisfiedCombining (4), (6), and (8) gives constitutive equation of the interface elementRepresenting the above equations in matrix form, the interface stress where indicates the interface stress, is the elasticity matrix, and is displacement matrix in local coordinate system.Applying (11) into (2) gives the interface element domain integral, the element stiffness matrix for elasticity, and the traction term which can be evaluated. Finally, the governing equation of IEM is expressed as

2.4. Thermal Effects on Interfacial Stress

For an isotropic material, temperature change results in a body expansion or shrinkage but no distortion. In other words, temperature change affects the normal strains but not shear strain. Thus, by adding the thermal strain to (6), the total strain is expressed asin which is the coefficient of thermal expansion (CTE) and indicates the temperature change. Further substitution of (13) into (7) yieldsFollowing the derived procedure described in Section 2.3 will result in interface stress with thermal effects aswhere is the thermal elasticity matrix. The stiffness matrix and force matrix become

2.5. Boundary Condition

The fixed boundary condition restricted the interface deformation in the region of block element 2; those terms are neglected in the elasticity matrix, i.e.,The often used symmetry boundary condition corresponds to a zero shear stress () which givesAll of the derivations of IEM fundamental equations above are based on plane stress condition. For the plane strain condition, Young’s modulus, Poisson ratio, and CTE in matrices and are replaced by

2.6. Polygon Interface Element Formulation

In the interface element method, the block element need not necessarily be limited to the shape of the triangle or quadrangle. An arbitrary polygon block element is available. The position of centroid of block element can be calculated by partition of element to be n-2 triangle as shown in Figure 3. The position of centroid for each triangle is easy to compute by taking average of the coordinate of the 3 points, i.e.,and the areas for each triangle and the polygon are given asThe coordinates of the centroid of polygon element are

A novel method called multiscale interface element method (MIEM) is proposed in this paper based on the concept of polygon element. As shown in Figure 4, the nodes 3, 4, 5, and 6 in element 1 are rearranged in a straight line and the shape of the polygon element 1 looks to be a quadrangle which connects three elements with one edge. The advanced mesh technique allows the triangle or quadrangle element to connect more than one element on each edge with multiple interfaces. The method can rapidly transform the element density from fine to coarse to perform a multiscale modeling, as shown in Figure 5.

3. Numerical Examples

3.1. Five-Layered Stack Structure

In this section, two examples are used to verify the validity of the developed model. The first example is a five-layered stack structure that was originally considered by Pao [4], as shown in Figure 6. The values of the geometrical parameters and material properties are listed in Table 1. The uniform temperature decrement is -65°C in this example. Half of the model is constructed due to geometric symmetry about the vertical axis in this paper. Figure 7 depicts the half of model created using the MIEM, refining the local mesh which is very flexible at the interface and near the free end. The interfacial peel and shear stress at layer 1 and 2 conjunct interfaces (Al/Cu) are presented in Figures 8 and 9, respectively. Figures 10 and 11 present the interfacial peel and shear stress at layer 3 and 4 conjunct interfaces (Solder/BeO), respectively. The results indicate that the stress values calculated at free end by Pao are generally higher than those calculated by MIEM and FEM. Because Pao’s model assumes each layer as a Bernoulli beam, Timoshenko effect is ignored.

3.2. LCR Model

The second example is a leadless chip resistor (LCR), as shown in Figure 12. Ghorbani [9] provided an analytical solution to this problem, and his data is used for this work to directly compare the results. The values of the geometrical parameters and material properties are listed in Table 2. The shear modulus is determined from resulting from all materials being assumed to be isotropic and homogeneous. Figure 13 depicts the model created using the MIEM; the grid size is divided into 1/3 of the original at the interface of solder joints to improve the accuracy of the solution. Under a uniform temperature increase of 75°C, the interfacial peel and shear stress distributions at the lower conjunct interfaces in the sample are presented in Figures 14 and 15, respectively. It can be seen that Ghorbani’s analytical model predicts shear stresses that reach zero at the free end; this condition is not satisfied by the MIE and FE models. The results by MIEM are consistent with the trend of FEM, but these values are closer to the analytical solution. Overall, those are reasonably in agreement with the solutions of Ghorbani [9].

4. Applying Multiscale Interface Element Method to CSP

Figure 16 shows the chip scale package (CSP) model schema. The half-length of the structure is considered, and the values of the geometrical parameters and material properties are listed in Table 3. The multiscale interface element method is employed to analyze thermomechanical behaviors of CSP. The MIEM model for the CSP is depicted in Figure 17. The thermal stress distribution of conjunct interfaces of the last two solder joints and another layer is investigated in this study, because these conjunct interfaces are most likely to crack. Figure 18 depicts the FEM model with 26600 elements for the CSP. A 4-node biquadratic plane stress quadrilateral element with reduced integration (CPS4R in ABAQUS) is used. The numerical results of MIEM model are verified by comparing with FEM model, as shown in Figures 1922. The solid and dashed lines denote the numerical solutions by using MIEM and FEM, respectively. Figures 19 and 20 present the interfacial peel stress and shear stress on FR4 PCB and solder joint conjunct interfaces. Figures 21 and 22 present the interfacial peel stress and shear stress on FR4 substrate and solder joint conjunct interfaces. These analyses are conducted with a temperature range of 165°C (-40°C to 125°C). As shown in these results, the solutions of these two models are consistent.

From the shear stress distribution, the final point of calculated stress curve by MIEM is closer to the singularity than the FEM, and the value is higher, because MIEM has the characteristics of the smaller interface mesh. From the peel stress distribution, a good agreement has been obtained between the two sets of results. MIEM can directly obtain the stress value at the interface, which is also one of the advantages of this method in dealing with multilayer structures.

5. Conclusions

This research developed a multiscale interface element model to evaluate the interfacial stresses in multilayer packaging structures. The peel and shear stress functions are calculated in terms of interface element fractional distribution. Two typical examples are used to verify the validity of the developed model. The present analyses provide better accuracy and rationality in interfacial peel and shear stresses.

The developed model was then applied to analyze thermal stress distributions of the solder joint conjunct interfaces in a CSP model. Using this efficient numerical technique, a very fine mesh pattern can be established around each conjunct interface without increasing the degree of freedom of the global FEM solution. In general, the results presented in this study have shown that the proposed MIEM algorithm provides a fast, direct, and accurate tool for simulating the stress analysis of multilayer structures.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was partially supported by the Advanced Institute of Manufacturing with High-Tech Innovations (AIM-HI) from The Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan. This research was also supported by R.O.C. MOST Foundation Contract No. MOST107-3017-F-194 -001.