Thermal Stress Analysis of Chip Scale Packaging by Using Novel Multiscale Interface Element Method

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 developedmodel,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.


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 finiteelement 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 Mathematical Problems in Engineering  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.

Derivations of Interface-Element Method
. .e 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 that The discretization form can be written as (2), where Ω  ,    ,   0 , and   denote the element domain, stress boundary, interface region, and interfacial stress.
. .Displacement Distribution in Interfacial Element.A fournode 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 ] . .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  1   per unit of the original length  1  = ℎ 1 in the local system, where  1  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  1  from block element 1 to the original length ℎ 1 .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 satisfied Combining ( 4), (6), and (8) gives constitutive equation of the interface element Representing 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 . .ermal 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 as in which  is the coefficient of thermal expansion (CTE) and Δ indicates the temperature change.Further substitution of ( 13) into ( 7) yields Following the derived procedure described in Section 2.3 will result in interface stress with thermal effects as where   is the thermal elasticity matrix.The stiffness matrix and force matrix become . .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 ( = 0) which gives . .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 as The 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 (2) 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.

Numerical Examples
. .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.
. .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 Mathematical Problems in Engineering 7  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   =   /2(1 + ]  ) 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].

Applying Multiscale Interface
Element Method to CSP   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.

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.

Figure 14 :
Figure 14: Interfacial peel stresses at upper and lower conjunct interfaces.

Table 1 :
Geometrical parameters and material properties.

Table 2 :
Geometrical parameters and material properties.

Table 3 :
The geometrical parameters and material properties for CSP model.