Energy-Based Thermal Fatigue Prediction in a Ball Grid Array

In a cooling system, a microelectronic component has been identified as the critical link. Since the power is repeatedly switched on and off, the component is subjected to thermal cycling. As a result, a crack grows through a solder joint and disconnects the chip from the printed circuit board, causing a malfunction. In the simulation, the solder lifetime in two ball grid assemblies is predicted. The prediction is based on the Darveaux energy-based model. The fatigue model evaluates damage based on an averaged energy dissipation density in a thin layer, where a crack is expected to grow.

This example is based on a model from the Nonlinear Structural Materials Module — Viscoplastic Solder Joints. Since the model contains several solder joints that are modeled with a viscoplastic material, many degrees of freedom are required in order to simulate the correct creep behavior in all elements. From the fatigue point of view, only the critical part of the model is of interest. In order to capture it, the concept of submodeling is used. This technique consists of two steps. In the first one, the full model is analyzed with a coarse mesh in order to capture the general trends and to identify the critical part of the model. In the second step a fine submodel containing the critical part is made and the study is resolved. The global effects from the full model are transferred to the submodel via appropriate boundary conditions.

Model Definition

The microelectronic component consists of a flat printed circuit board that is covered with a thin copper layer. Two microprocessors are connected to the copper layer with a ball grid array of 60Sn-40Pb solder joints; see Figure 1.

Both microprocessors generate power when they are switched on and when they are in stand-by. This generates heat throughout its operational lifetime. The power cycle is at its maximum, 5·107 W/m3, during 4 h and at its minimum 1·107 W/m3, during 2 h. The switch between high and low power is not instantaneous, but takes few minutes.

Figure 1: Geometry of the microelectronic component.

The elastic and thermal properties of the material can be taken from the built-in material library in COMSOL Multiphysics. The heat transfer coefficient for all free surfaces can be approximated with 10 W/(m2·K). The nonlinear behavior of the solder material follows the Anand material model with parameters summarized in Table 1.

Table 1: Constants of the Anand model.
Property	value	description
A	1.49·107 1/s	Pre-exponential factor
Q	90046 J/mol	Activation energy/Boltzmann constant
ξ	11	Stress multiplier
m	0.303	Strain rate sensitivity of stress
ssat	80.42 MPa	Coefficient for deformation resistance saturation
sinit	56.33 MPa	Initial value of deformation resistance
h0	2640.75 MPa	Hardening coefficient
a	1.34	Strain rate sensitivity of hardening
n	0.0231	Sensitivity for deformation resistance

The Darveaux fatigue model is representative for life prediction of the solder material. The model combines life contributions from crack initiation and crack propagation using the expression

where N is the fatigue life given in number of cycles, ΔWave is the averaged dissipated energy density in a fatigue cycle, a is the distance the crack needs to propagate for the failure to occur, and K1, k2, K3, k4, and Wref are material constants. The numerical values of the material constants are given in Table 2.

Table 2: Fatigue material parameters.
Property	Value	Description
K1	13173	Crack initiation energy coefficient
k2	-1.45	Crack initiation energy exponent
K3	3.92·10-7 in	Crack propagation energy coefficient
k4	1.12	Crack propagation energy exponent
Wref	1 psi3	Reference energy density

The distance a is in this analysis be taken to be 2.6457·10−4 in. This is based on the assumption that the problem is not symmetric and a crack is expected to start on one side of the joint only and not all around the joint at the same time.

The concept of submodeling is utilized in this example. This technique requires that first an analysis of the full model is performed in order to capture general trends, followed by an analysis of a submodel that is studied in detail. The following steps are done:

A coupled structural and thermal analysis is performed during four load cycles on the full model. Since the model contains several solder joints, a coarse mesh is used. All joints are meshed in the same way in order to minimize numerical discrepancies that are mesh dependent.

A fatigue prediction is made on the fourth cycle. The energy dissipation volume average and corresponding life is evaluated for each individual solder joint. The critical joint is identified.

A submodel of the critical joint with a fine mesh is created. The displacments from step 1 are prescribed as on the boundaries where the submodel is cut out of the full model. Similarly, the temperature from step 1 is prescribed as a in the heat transfer analysis of the submodel. This is done for all time steps of the four simulated cycles.

A fatigue analysis is performed on the critical solder joint in the submodel. A life prediction is made based on the energy dissipation in a 50 μm thick layer. Two layers are evaluated. One that is in connection with the copper side and one that is in connection with the microchip side.

The submodel is shown in Figure 2. The purple color denotes the boundaries where the results from the structural analysis of the global model are prescribed. The solder joint is divided in three domains: a central one and two domains close to the interface to the other materials.

Figure 2: Submodel containing the critical joint.

Results and Discussion

The fatigue life prediction for all joints is shown in Figure 3. For both microchips, the critical joints are located in the corners of each ball grid array. This is expected, since those joints experience highest strains due to the differences in thermal properties. The solder joints of the microchip with the larger ball grid array shows a shorter life than the joints of the other microchip. The life prediction for all four corner joints is about the same, 104.8 cycles.

Figure 3: Fatigue life prediction for all joints based on the global model.

One of the four critical corner joints is reanalyzed in the submodel. In order to verify that the temperature is correctly prescribed in the submodel a comparison of the temperature history in a point on the upper sider of the joint is shown in Figure 4. The results of both models are in prefect match. Note that much fewer data points are stored for the submodel.

Figure 4: Temperature history for both the full model and the submodel.

A comparison of the dissipated energy in a point on the upper side of the joint is shown in Figure 5. The results differs between the two models. This difference is caused by the difference in the mesh. The finer mesh in the submodel is better than the coarser mesh of the full model in capturing the strain gradient close to the upper side of the solder joint. Note that much fewer data points are stored in the submodel for the first three cycles, hence, the presented results are not a smooth. This does not affect the accuracy of the solution.

Figure 5: Comparison of the dissipated creep energy in both models.

The results of the fatigue analysis of the critical joint in the submodel are shown in Figure 6 and Figure 7. In the first figure, the fatigue life is based on the energy dissipation volume average evaluated over the whole joint, while in the second figure the volume average is performed over separate domains. The fatigue life predicted by different models is summarized in Table 3.

Table 3: Fatigue Life based on different modeling techniques.
Evaluation Method	Fatigue Life (cycles)	Source
Entire joint in the full model	5.8·104	Figure 3
Entire joint in the submodel	8.1·104	Figure 6
Thin layer in the submodel	9.1·103	Figure 7

The difference between the fatigue life prediction of the full model and of the submodel is small. It is however observed in real applications that a crack grows through a solder joint close to the interface and therefore a thin layer is required. If a thin layer were to be created in all joints of the full model, the simulation would require large computational resources. In the current example, the full model consists of about 170,000 DOFs while the submodel consists of 40,000 DOFs.

Figure 6: Fatigue prediction based on the volume average of the entire joint.

Figure 7: Fatigue prediction based on the volume average of different domains.

Notes About the COMSOL Implementation

In submodeling, the results are mapped from one component — a full model — to another one — a submodel. This is possible through the operator. In this example, the mapping of results is simple since it is 1 to 1. The operator allows also for mapping into other shapes.

Often, a load history in a fatigue study is repeatable. In order to write the cyclic function in a compact way, a modulo function can be used. This function calculates the reminder of a number, dividend, divided by an other number, divisor. In COMSOL Multiphysics it is defined as mod(f,p), where the first argument is the dividend and the second one is the divisor.

From the numerical point of view, an abrupt change in a load parameter can be challenging. In the current example, a sudden increase or decrease of the power provides such a challenge. In such a case it is favorable to use a smooth function that changes from one value to an other. This is possible in via flc2hs(t,p) function that is a Heaviside function with a smooth second derivative. The first argument defines a position when the step takes place and the second argument defines an interval on each side of the step position where the smooth transition takes place.

Fatigue_Module/Energy_Based/viscoplastic_solder_joints_fatigue

Modeling Instructions

In this example you will start from an existing model which is an example in the Nonlinear Structural Materials Module.

Application Libraries

From the menu, choose .

In the window, select in the tree.

Click

Full Model

In the window, click .

In the window for , type Full Model in the text field.

Define the crack size that will be used in crack models.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
lcrack	2.6457e-4[in]	6.7201E-6 m	Crack size

Define the power load cycle.

Analytic 1 (power)

In the window, click .

In the window for , locate the section.

In the text field, type (flc2hs(x-0.1,0.1)*50)*(x<6)-flc2hs(mod(x,6)-4.1,0.1)*40+(flc2hs(mod(x,6)-0.1,0.1)*40+10)*(x>=6).

Full Model (comp1)

Create a operator that will be used when the submodel is set up.

Definitions

General Extrusion 1 (genext1)

In the window, expand the node.

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Full Model: Load History

In the window, click .

In the window for , type Full Model: Load History in the text field.

Step 1: Time Dependent

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type 0 0.005 range(0.025,0.025,0.5) range(0.75,0.25,3.75) 3.975 4+{range(0,0.025,0.5) range(0.75,0.25,2)} 6+{range(0.025,0.025,0.5) range(0.75,0.25,3.75) 3.975 4+{range(0,0.025,0.5) range(0.75,0.25,2)}} 12+{range(0.025,0.025,0.5) range(0.75,0.25,3.75) 3.975 4+{range(0,0.025,0.5) range(0.75,0.25,2)}} 18+{range(0.025,0.025,0.5) range(0.75,0.25,3.75) 3.975 4+{range(0,0.025,0.5) range(0.75,0.25,2)}}.

Solver Configurations

In the window, expand the node.

Step 2: Time Dependent 2

In the text field, enter the same steps as .

In the window, expand the node.

Right-click and choose .

Results

Stress (solid)

Perform a fatigue study on the full model in order to find out which solder joint is the critical one.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Find the subsection. In the table, clear the check box for .

Click in the window toolbar.

In the toolbar, click

to close the window.

Fatigue (ftg)

Energy-Based 1

Right-click and choose the domain evaluation .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the a text field, type lcrack.

Root

In the window, expand the node.

Global Definitions

In the window, expand the node.

Solder, 60Sn-40Pb (mat4)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Crack initiation energy coefficient	K1_Darveaux	13173	1	Darveaux
Crack initiation energy exponent	k2_Darveaux	-1.45	1	Darveaux
Crack propagation energy coefficient	K3_Darveaux	3.92e-7[in]	m	Darveaux
Crack propagation energy exponent	k4_Darveaux	1.12	1	Darveaux
Reference energy density	Wref_Darveaux	1[psi]	J/m³	Darveaux