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

Authors: D. S. Liu; Y. W. Chen; C. Y. Tsai
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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.


---

## Body

## 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 that(1)∭ΩσijδeijdΩ=-∬S0TiδuidS=∭ΩfiδuidΩ+∬SσpiδuidSThe discretization form can be written as (2), where Ωe, Sσe, S0e, and Ti denote the element domain, stress boundary, interface region, and interfacial stress.(2)-∑e=1n∬S0eTiδuidS=∑e=1n∭ΩefiδuidΩ+∬SσepiδuidS


### 2.2. Displacement Distribution in Interfacial Element

A four-node 2D interface element for illustrated purpose is shown in Figure1, in which n⃑ and s⃑ 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 (3)u=N×ugu=uvTug=ugvgθTN=10yg-y01x-xgwhere u denotes the displacement vector at any point in the block element, ug is the degree of freedom at the mass of center, and N is the shape function. With respect to the local system, the relative displacements at interface element are then derived as (4)δlocal1+δlocal2=δn1+δn2δs1+δs2T=-L1u1+L2u2=-L1N1ug1-N2ug2In the above equation, δn is normal displacement, δs is the shear displacement, andL is the coordinate transformation matrix which reads(5)L=cos⁡n⃑,x⃑cos⁡n⃑,y⃑cos⁡s⃑,x⃑cos⁡s⃑,y⃑L1=-L2=LFigure 1
Schematic diagram of IEM model.


### 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 Figure2, for any point M on interface surface, the normal strain contributed from block element 1 is expressed as the change in length δn1 per unit of the original length P1M¯=h1 in the local system, where P1M¯ 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 δs1 from block element 1 to the original length h1. 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:(6)εni=δnihiγsi=δsihiIn plane stress situation, the components of stress can be evaluated by the following equation.(7)σni=Eiεni=Eiδnihiτsi=Giγsi=E21+viδsihiwhere σni and τsi are the normal stress and shear stress contributed by block elementi. The interface stress should be balanced due to the system being in an equilibrium state. It satisfied(8)σn1=σn2=σnτs1=τs2=τsCombining (4), (6), and (8) gives constitutive equation of the interface element(9)δn1+δn2=σn1h1E1+σn2h2E2=σnE2h1+E1h2E1E2σn=E1E2E2h1+E1h2δn1+δn2δs1+δs2=21+v1h1τs1E1+21+v2h2τs2E2=τs21+v1h1E2+21+v2h2E1E1E2τs=E1E221+v1h1E2+21+v2h2E1δs1+δs2Representing the above equations in matrix form, the interface stress (10)T=Dδ=-DL1u1+L2u2=-DL1N1ug1+L2N2ug2=-DL1N1ug1-N2ug2=-DLN∗ug∗where T indicates the interface stress, D is the elasticity matrix, and δ is displacement matrix in local coordinate system.(11)T=σnτsT;δ=δniδsiT;D=dn00ds;N∗=N1-N2;ug∗=ug1ug2T;dn=E1E2E1h2+E2h1ds=E1E22E1h21+ν2+2E2h11+ν1Applying (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(12)KU=RK=∑e=1nKe,Ke=∫sjN∗TLTDLN∗dSU=ug1ug2ug3⋯TR=∑e=1nRe,Re=∬ΩeN∗TfdΩ+∫SσeN∗TPdSFigure 2
Schematic diagram of deformation.


### 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 as(13)εni=δnihi-αiΔTγsi=δsihiin which α is the coefficient of thermal expansion (CTE) and ΔT indicates the temperature change. Further substitution of (13) into (7) yields(14)σni=Eiεni=Eiδnihi-EiαiΔTτsi=Giγsi=E21+viδsihiFollowing the derived procedure described in Section 2.3 will result in interface stress with thermal effects as(15)T=-DLN∗ug∗-DTΔTDT=E1E2α1h1+α2h2E1h2+E2h10where DT is the thermal elasticity matrix. The stiffness matrix and force matrix become(16)KU=RK=∑e=1nKe,Ke=∫sjN∗TLTDLN∗dSU=ug1ug2ug3⋯TR=∑e=1nRe,Re=∬ΩeN∗TfdΩ+∫SσeN∗TPdS-∫S0eN∗TLTDTΔTdS


### 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.,(17)D=E1h100E12h11+v1The often used symmetry boundary condition corresponds to a zero shear stress (ds=0) which gives(18)D=E1h1000All 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 D and DT are replaced by(19)E=E1-v2v=v1-vα=a1+v


### 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 Figure3. The position of centroid for each triangle is easy to compute by taking average of the coordinate of the 3 points, i.e.,(20)xgi=x1i+x2i+x3i3ygi=y1i+y2i+y3i3and the areas for each triangle and the polygon are given as(21)ai=12det⁡1x1iy1i1x2iy2i1x3iy3iA=∑i=1n-2aiThe coordinates of the centroid of polygon element are(22)X-g=∑i=1n-2xgiaiAY-g=∑i=1n-2ygiaiAFigure 3
The polygon element.A novel method called multiscale interface element method (MIEM) is proposed in this paper based on the concept of polygon element. As shown in Figure4, 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.Figure 4
Multiscale interface element.Figure 5
The multiscale interface element model.


## 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 that(1)∭ΩσijδeijdΩ=-∬S0TiδuidS=∭ΩfiδuidΩ+∬SσpiδuidSThe discretization form can be written as (2), where Ωe, Sσe, S0e, and Ti denote the element domain, stress boundary, interface region, and interfacial stress.(2)-∑e=1n∬S0eTiδuidS=∑e=1n∭ΩefiδuidΩ+∬SσepiδuidS


## 2.2. Displacement Distribution in Interfacial Element

A four-node 2D interface element for illustrated purpose is shown in Figure1, in which n⃑ and s⃑ 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 (3)u=N×ugu=uvTug=ugvgθTN=10yg-y01x-xgwhere u denotes the displacement vector at any point in the block element, ug is the degree of freedom at the mass of center, and N is the shape function. With respect to the local system, the relative displacements at interface element are then derived as (4)δlocal1+δlocal2=δn1+δn2δs1+δs2T=-L1u1+L2u2=-L1N1ug1-N2ug2In the above equation, δn is normal displacement, δs is the shear displacement, andL is the coordinate transformation matrix which reads(5)L=cos⁡n⃑,x⃑cos⁡n⃑,y⃑cos⁡s⃑,x⃑cos⁡s⃑,y⃑L1=-L2=LFigure 1
Schematic diagram of IEM model.


## 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 Figure2, for any point M on interface surface, the normal strain contributed from block element 1 is expressed as the change in length δn1 per unit of the original length P1M¯=h1 in the local system, where P1M¯ 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 δs1 from block element 1 to the original length h1. 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:(6)εni=δnihiγsi=δsihiIn plane stress situation, the components of stress can be evaluated by the following equation.(7)σni=Eiεni=Eiδnihiτsi=Giγsi=E21+viδsihiwhere σni and τsi are the normal stress and shear stress contributed by block elementi. The interface stress should be balanced due to the system being in an equilibrium state. It satisfied(8)σn1=σn2=σnτs1=τs2=τsCombining (4), (6), and (8) gives constitutive equation of the interface element(9)δn1+δn2=σn1h1E1+σn2h2E2=σnE2h1+E1h2E1E2σn=E1E2E2h1+E1h2δn1+δn2δs1+δs2=21+v1h1τs1E1+21+v2h2τs2E2=τs21+v1h1E2+21+v2h2E1E1E2τs=E1E221+v1h1E2+21+v2h2E1δs1+δs2Representing the above equations in matrix form, the interface stress (10)T=Dδ=-DL1u1+L2u2=-DL1N1ug1+L2N2ug2=-DL1N1ug1-N2ug2=-DLN∗ug∗where T indicates the interface stress, D is the elasticity matrix, and δ is displacement matrix in local coordinate system.(11)T=σnτsT;δ=δniδsiT;D=dn00ds;N∗=N1-N2;ug∗=ug1ug2T;dn=E1E2E1h2+E2h1ds=E1E22E1h21+ν2+2E2h11+ν1Applying (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(12)KU=RK=∑e=1nKe,Ke=∫sjN∗TLTDLN∗dSU=ug1ug2ug3⋯TR=∑e=1nRe,Re=∬ΩeN∗TfdΩ+∫SσeN∗TPdSFigure 2
Schematic diagram of deformation.


## 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 as(13)εni=δnihi-αiΔTγsi=δsihiin which α is the coefficient of thermal expansion (CTE) and ΔT indicates the temperature change. Further substitution of (13) into (7) yields(14)σni=Eiεni=Eiδnihi-EiαiΔTτsi=Giγsi=E21+viδsihiFollowing the derived procedure described in Section 2.3 will result in interface stress with thermal effects as(15)T=-DLN∗ug∗-DTΔTDT=E1E2α1h1+α2h2E1h2+E2h10where DT is the thermal elasticity matrix. The stiffness matrix and force matrix become(16)KU=RK=∑e=1nKe,Ke=∫sjN∗TLTDLN∗dSU=ug1ug2ug3⋯TR=∑e=1nRe,Re=∬ΩeN∗TfdΩ+∫SσeN∗TPdS-∫S0eN∗TLTDTΔTdS


## 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.,(17)D=E1h100E12h11+v1The often used symmetry boundary condition corresponds to a zero shear stress (ds=0) which gives(18)D=E1h1000All 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 D and DT are replaced by(19)E=E1-v2v=v1-vα=a1+v


## 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 Figure3. The position of centroid for each triangle is easy to compute by taking average of the coordinate of the 3 points, i.e.,(20)xgi=x1i+x2i+x3i3ygi=y1i+y2i+y3i3and the areas for each triangle and the polygon are given as(21)ai=12det⁡1x1iy1i1x2iy2i1x3iy3iA=∑i=1n-2aiThe coordinates of the centroid of polygon element are(22)X-g=∑i=1n-2xgiaiAY-g=∑i=1n-2ygiaiAFigure 3
The polygon element.A novel method called multiscale interface element method (MIEM) is proposed in this paper based on the concept of polygon element. As shown in Figure4, 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.Figure 4
Multiscale interface element.Figure 5
The multiscale interface element model.


## 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.Table 1
Geometrical parameters and material properties.
Layer Material Thickness(mm) CTE(ppm/°F) Young’s Modulus(GPa) Poissonratio 1 Al 1.626 23.1 70 0.345 2 Cu 0.635 16.5 126 0.343 3 Solder 0.152 28.4 22.17 0.400 4 BeO 0.635 6.3 345 0.300 5 Si 0.254 2.6 120 0.420Figure 6
The five-layered stack structure [4].Figure 7
The multiscale interface element model.Figure 8
The interfacial peel stresses at layer 1 and 2 conjunct interfaces (Al/Cu).Figure 9
The interfacial shear stresses at layer 1 and 2 conjunct interfaces (Al/Cu).Figure 10
The interfacial peel stresses at layer 3 and 4 conjunct interfaces (Solder/BeO).Figure 11
The interfacial shear stresses at layer 3 and 4 conjunct interfaces (Solder/BeO).


### 3.2. LCR Model

The second example is a leadless chip resistor (LCR), as shown in Figure12. 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 Gi=Ei/2(1+νi) 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].Table 2
Geometrical parameters and material properties.
Layer Width(mm) Thickness(mm) CTE(ppm/°C) Young’s Modulus(GPa) Poisson ratio Resistor 3.25 0.65 2.8 131.00 0.30 Solder Joints 0.76 0.12 22.0 64.54 0.40 Substrate 3.25 1.23 17.0 22.00 0.28Figure 12
Schematic of an LCR assembly [9].Figure 13
The interface element model.Figure 14
Interfacial peel stresses at upper and lower conjunct interfaces.Figure 15
Interfacial shear stresses at upper and lower conjunct interfaces.


## 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.Table 1
Geometrical parameters and material properties.
Layer Material Thickness(mm) CTE(ppm/°F) Young’s Modulus(GPa) Poissonratio 1 Al 1.626 23.1 70 0.345 2 Cu 0.635 16.5 126 0.343 3 Solder 0.152 28.4 22.17 0.400 4 BeO 0.635 6.3 345 0.300 5 Si 0.254 2.6 120 0.420Figure 6
The five-layered stack structure [4].Figure 7
The multiscale interface element model.Figure 8
The interfacial peel stresses at layer 1 and 2 conjunct interfaces (Al/Cu).Figure 9
The interfacial shear stresses at layer 1 and 2 conjunct interfaces (Al/Cu).Figure 10
The interfacial peel stresses at layer 3 and 4 conjunct interfaces (Solder/BeO).Figure 11
The interfacial shear stresses at layer 3 and 4 conjunct interfaces (Solder/BeO).


## 3.2. LCR Model

The second example is a leadless chip resistor (LCR), as shown in Figure12. 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 Gi=Ei/2(1+νi) 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].Table 2
Geometrical parameters and material properties.
Layer Width(mm) Thickness(mm) CTE(ppm/°C) Young’s Modulus(GPa) Poisson ratio Resistor 3.25 0.65 2.8 131.00 0.30 Solder Joints 0.76 0.12 22.0 64.54 0.40 Substrate 3.25 1.23 17.0 22.00 0.28Figure 12
Schematic of an LCR assembly [9].Figure 13
The interface element model.Figure 14
Interfacial peel stresses at upper and lower conjunct interfaces.Figure 15
Interfacial shear stresses at upper and lower conjunct interfaces.


## 4. Applying Multiscale Interface Element Method to CSP

Figure16 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 19–22. 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.Table 3
The geometrical parameters and material properties for CSP model.
Layer Thickness(mm) Elastic Modulus(GPa) Shear Modulus(GPa) CTE(ppm/°C) Poisson’s ratio 1 FR4 PCB 1.2 22 8.6 18 0.28 2 Solder joint 0.28 22.2 7.9 28.4 0.4 3 FR4 Substrate 0.25 22 8.6 18 0.28 4 Underfill 0.05 5 1.9 60 0.35 Solder bump 22.2 7.9 28.4 0.4 5 Si die 0.5 160 61.5 2.8 0.3Figure 16
The schematic of chip scale package (CSP) model.Figure 17
The MIEM model.Figure 18
The FEM model (total number of elements: 26600).Figure 19
Interfacial peel stresses on FR4 PCB and solder joint conjunct interfaces.Figure 20
Interfacial shear stresses on FR4 PCB and solder joint conjunct interfaces.Figure 21
Interfacial peel stresses on FR4 substrate and solder joint conjunct interfaces.Figure 22
Interfacial shear stresses on FR4 substrate and solder joint conjunct interfaces.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.



---

*Source: 4962498-2018-10-21.xml*