Fatigue Module

Application Library Manual

Fatigue Module Application Library Manual
© 1998–2015 COMSOL
Protected by U.S. Patents listed on [Link]/patents, and U.S. Patents 7,519,518; 7,596,474;
7,623,991; 8,457,932; 8,954,302; 9,098,106; and 9,146,652. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, COMSOL Multiphysics, Capture the Concept, COMSOL Desktop, LiveLink, and COMSOL
Server are either registered trademarks or trademarks of COMSOL AB. All other trademarks are the property
of their respective owners, and COMSOL AB and its subsidiaries and products are not affiliated with,
endorsed by, sponsored by, or supported by those trademark owners. For a list of such trademark owners,
see [Link]/trademarks.
Version: COMSOL 5.2

Contact Information
Visit the Contact COMSOL page at [Link]/contact to submit general
inquiries, contact Technical Support, or search for an address and phone number. You can
also visit the Worldwide Sales Offices page at [Link]/contact/offices for
address and contact information.

If you need to contact Support, an online request form is located at the COMSOL Access
page at [Link]/support/case. Other useful links include:

• Support Center: [Link]/support

• Product Download: [Link]/product-download
• Product Updates: [Link]/support/updates
• COMSOL Blog: [Link]/blogs
• Discussion Forum: [Link]/community
• Events: [Link]/events
• COMSOL Video Gallery: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part number: CM023202

 Solved with COMSOL Multiphysics 5.2

Accelerated Life Testing

Introduction
Fatigue testing of nonlinear materials with a creep mechanism is a time consuming
process. In accelerated life testing the experiment time is greatly reduced by subjecting
the material to testing conditions in excess of the operating one. In the model an
aggressive thermal load cycle is simulated and its effect on the fatigue life of a solder
joint is examined.

This example demonstrates how to evaluate fatigue driven by a specific strain or energy
quantity and therefore a simple schematic representation of an electronic component
is used. Moreover, only one cycle is simulated as opposed to several that are required
to obtain a steady state cycle when working with nonlinear materials. This
simplification is motivated by the fact that the purpose of this model is to show how
to evaluate fatigue based on a user defined variable. The required strain and energy
variables are defined explicitly via ordinary differential equations and calculated during
the simulation of the thermal load cycle.

Model Definition
Microelectronic components consist of several materials. When subjected to thermal
changes the differences in coefficient of thermal expansion introduces a stress
concentration caused by the uneven material deformation. With the consecutive
heating and cooling, the stress state changes back and forth until the component finally
fails in fatigue.

The example uses a schematic geometry that captures the structural function of the
main components rather than the accurate geometry of the microelectronic
component. A symmetric model is used and half of its geometry is shown in Figure 1.
The total size of the resistor and of the printed circuit board (PCB) is 4 × 0.5 mm2.
The size of the solder is 0.5 × 0.25 mm2. It is also assumed that the model is in plane
strain conditions.

The resistor is made of alumina, the printed circuit board is a fiberglass based laminate
and the solder joint is an SnAgCu alloy that responds both linearly and nonlinearly to
an applied load. The elastic and the thermal properties of the materials are given in
Table 1.

1 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Resistor

Solder

PCB

Figure 1: Geometry of the electronic component. The dashed line denotes symmetry.

TABLE 1: ELASTIC AND THERMAL PROPERTIES OF MATERIALS

MATERIAL YOUNG’S MODULUS (GPA) POISSON’S RATIO COEFFICIENT OF THERMAL

EXPANSION (PPM/°C)

Fiberglass 22 0.4 21
SnAgCu 50 0.4 21
Alumina 300 0.22 8

The solder material nonlinearity is creep, which is a material behavior where

deformation changes over the time although the stress remain constant. Generally
materials creep with different rates in three distinct phases: primary creep, secondary
creep and tertiary creep; see Figure 2.

Figure 2: Creep strain development at constant stress. Primary creep denoted with 1,
secondary with 2 and tertiary with 3.

The secondary creep is also called the steady state creep since the strain rate is constant.
In the representation of the solder material the primary and tertiary creeps are

2 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

disregarded while the steady state creep is represented with a double Norton law
according to

σe nI σ e n II
ε c = A I  ------ + A II  ------
·
σn σn

·
where ε c is the creep rate, σe is the effective stress and AI, nI, σn, AII and nII are the
material constants given by Table 2. The first term represents a creep mechanism
observed at low stresses while the second term describes the dominating creep
behavior at high stresses, see Figure 3.
TABLE 2: CREEP MATERIAL CONSTANTS

MATERIAL CONSTANT VALUE

AI 8.03e-12 [1/s]
nI 3
σn 1 [MPa]
AII 1.96e-23 [1/s]
nII 12

Figure 3: Creep model for the solder SnAgCu solder material.

Two fatigue models are evaluated. In the first, the lifetime is based on a Coffin-Manson
type model. The development of the shear strain controls the fatigue according to

– 0.61
Δγ II = 0.587 ⋅ N

where ΔγII is the range of the creep shear strain of the second creep mechanism
experienced in a load cycle and N is the fatigue lifetime.

3 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

The second fatigue evaluation is dependent on the energy dissipation and follows a
Morrow type criterion according to

6 – 0.79
ΔW II = 74 ⋅ 10 ⋅ N

where ΔWII is the dissipated creep energy of the second creep mechanism during a load
cycle.

Results and Discussion

First the user defined creep strains and energies are verified. The creep strain calculated
in the structural analysis is a summation of creep contributions from different
mechanisms and thus

I II
εc = εc + εc
I II
where ε c is the creep strain of the first mechanism and ε c is the creep strain
contribution from the second mechanism. The individual contribution of each creep
mechanism is calculated in two separate ODE interfaces. A composition of creep
results is shown in Figure 4.

Figure 4: Creep strain history. The effective creep strain is shown.

4 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

From the figure it is clear that the creep strain calculated using the Nonlinear
Structural Materials Module equals the sum of the creep strains calculated with the
user defined method.

The creep energy dissipation can also be calculated using the additive decomposition
since

I II
δW c = δε c σ e = δε c σ e + δε c σ e

where δ denotes an increment and Wc denotes the energy dissipation density. The
results of the energy comparison in Figure 5 shows a perfect agreement between values
calculated by the built-in materials and the user defined material.

Figure 5: The energy dissipation history.

The fatigue life predicted by the Coffin-Manson type model and the Morrow-type
model is shown in Figure 6 and Figure 7 respectively. Both models predict
approximately the same life and indicate that the upper left corner is the critical point
in the assembly.

5 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Figure 6: Fatigue life based on the shear strain.

6 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Figure 7: Fatigue life based on the dissipated energy.

Notes About the COMSOL Implementation

Fatigue can be evaluated with many different models. Several models share the same
mathematical from. This is the case for the Coffin-Manson based and Morrow based
models where the only difference is the controlling strain or energy measure. In rocks
and rubbers it is commonly the elastic strain energy or the total strain energy that
controls fatigue while in materials subjected to cyclic creep the dissipated energy
controls the fatigue instead. The Coffin-Manson model that initially related plastic
strain to fatigue life has been modified by several researchers that propose a variety of
different strain measures as the fatigue controlling mechanism.

In COMSOL the Morrow and the Coffin-Manson relations

´ m
ΔW = W f ⋅ N
Δε ´ c
------ = ε f ⋅ N
2

can be modified so that the fatigue evaluation is based on a user defined energy density,
W, or strain, ε. If the fatigue controlling variable is already defined in a physics interface

7 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

it can be directly specified in the fatigue model node. The variable must be specified
together with the corresponding physics identifier for example, [Link] when using
the elastic strain energy density. When the fatigue controlling variable is not defined in
a physics interface, it can be calculated using an ODE Interface. This procedure is
demonstrated in the model where the fatigue is dependent on a specific creep strain
and a specific dissipated energy density.

The basic form for an ODE is

2
∂ u ∂u
ea + da = f (1)
∂t
2 ∂t

where u is the field variable, t is the time, ea denotes the mass coefficient matrix, da
denotes the damping coefficient matrix and f is the source term. The development of
creep following the Norton’s law is defined with

ij
dε c σ e n 3 s ij
= A  --------- ⋅ --- ⋅ ------
dt  σ ref 2 σ e
(2)
ij ij
dε c 2 dε c dε c
= --- ⋅
dt 3 dt dt

where ij denotes a specific component, σe is the effective von Mises stress, sij is the
deviatoric stress, and A, σref and n are material constants. Since the problem is a 2D
plane stress analysis only strains in x, y, z and xy directions can develop. This gives 4
dependent variables: ecx, ecy, ecz, ecxy. By comparing the first relation of
Equation 2 with Equation 1 it can be identified that there is no contribution from the
mass matrix, the damping matrix is an identity matrix and the source term equals the
right-hand side of the first relation in Equation 2. This relation contains both the
deviatoric stress tensor and the effective stress. Although both are already defined in
the Solid Mechanics interface, intuitively it might seem that they can be used in the
definition of the source term. COMSOL requires however that the derivatives of all
user defined variables can be evaluated at all calculation steps of the analysis. At zero
stress the numerical derivative of the effective stress

2 2 2 2
σe = σ x + σ y + σ z + σ x σ y + σ x σ z + σ y σ z + 3τ xy

tries to evaluate a negative power of zero. A remedy to this challenge is an own

definition of the effective stress with a small addition of stress.

8 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

user 2 2 2 2
σe = σ x + σ y + σ z + σ x σ y + σ x σ z + σ y σ z + 3τ xy + 1

The nonzero stress ensures a finite derivative at stress free conditions. In the example
a stress addition of 1Pa was chosen. This gives a negligible contribution to the effective
stress since the resulting stresses have the order of magnitude of MPa.

In addition to the dependent variables of the creep strain components also one
additional dependent variable, ece, for the effective creep strain defined by the second
relation in Equation 2 is needed. The source term for this relation can be defined with
the COMSOL built-in derivative operator. For example a partial derivative of u with
respect to x is simply defined d(u,x). By putting all together, the time derivatives of
the five dependent variables, ecx, ecy, ecz, ecxy, ece, becomes:

• alpha*[Link]
• alpha*[Link]
• alpha*[Link]
• alpha*[Link]
• (2/
3*(d(ecx,TIME)^2+d(ecy,TIME)^2+d(ecz,TIME)^2+2*(d(ecxy,TIME)^2))+1e-20)^0.5

where alpha=3/2*A/s_mises*(s_mises/s_ref)^n and

s_mises=sqrt([Link]^2+[Link]^2+[Link]^[Link]*[Link]-sol
[Link]*[Link]*[Link]+3*[Link]^2+(1e-6[MPa])^2).

The time derivative of the dissipated creep energy density is

dWc dε
----------- = --------c ⋅ σ e
dt dt

By using the derivative operator, the source term for the dependent variable, Wc, is
simply d(ece,TIME)*s_mises.

The calculation of the strain variables and the energy density must be made in separate
ODE Interfaces since the units of both variables differ.

Application Library path: Fatigue_Module/Strain_Life/

accelerated_life_testing

9 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Blank Model.

ROOT
On the Home toolbar, click Add Component and choose 2D.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 2.
4 In the Height text field, type 0.5.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 0.5.
4 In the Height text field, type 0.25.
5 Locate the Position section. In the y text field, type 0.5.

Rectangle 3 (r3)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 2.
4 In the Height text field, type 0.5.
5 Locate the Position section. In the y text field, type 0.75.

Point 1 (pt1)
1 On the Geometry toolbar, click Primitives and choose Point.
2 In the Settings window for Point, locate the Point section.

10 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

3 In the x text field, type 0.25.

4 In the y text field, type 0.625.
5 Right-click Point 1 (pt1) and choose Build Selected.
6 Click the Zoom Extents button on the Graphics toolbar.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Solid Mechanics (solid).
4 Click Add to Component in the window toolbar.
5 On the Home toolbar, click Add Physics to close the Add Physics window.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

A_I 8.03e-12 [1/s] 8.03E-12 1/s Low stress creep rate
coefficient
n_I 3 3 Low stress creep rate exponent
A_II 1.96e-23 [1/s] 1.96E-23 1/s High stress creep rate
coefficient
n_II 12 12 High stress creep rate exponent
s_ref 1[MPa] 1E6 Pa Reference stress

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 In the Function name text field, type thremLC.
4 In the table, enter the following settings:

t f(t)
0 25
15 100

11 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

t f(t)
30 100
45 25
60 25
5 Locate the Units section. In the Arguments text field, type min.
6 In the Function text field, type degC.

SOLID MECHANICS (SOLID)

Linear Elastic Material 1

1 In the Model Builder window’s toolbar, click the Show button and select Advanced
Physics Options in the menu.
2 In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid)
click Linear Elastic Material 1.
3 In the Settings window for Linear Elastic Material, click to expand the Energy
dissipation section.
4 Locate the Energy Dissipation section. Select the Calculate dissipated energy check
box.

Thermal Expansion 1
1 On the Physics toolbar, click Attributes and choose Thermal Expansion.
2 In the Settings window for Thermal Expansion, locate the Model Inputs section.
3 In the T text field, type thremLC(t).

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

Creep 1
1 On the Physics toolbar, click Attributes and choose Creep.
2 Select Domain 2 only.
3 In the Settings window for Creep, locate the Creep Data section.
4 In the A text field, type A_I.
5 In the σref text field, type s_ref.
6 In the n text field, type n_I.

12 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

Creep 2
1 On the Physics toolbar, click Attributes and choose Creep.
2 Select Domain 2 only.
3 In the Settings window for Creep, locate the Creep Data section.
4 In the A text field, type A_II.
5 In the σref text field, type s_ref.
6 In the n text field, type n_II.

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundaries 11 and 12 only.

Fixed Constraint 1
1 On the Physics toolbar, click Points and choose Fixed Constraint.
2 Select Point 8 only.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, type PCB in the Label text field.
3 Select Domain 1 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 22 [GPa] Pa Basic
Poisson's ratio nu 0.4 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 21e-6 1/K Basic

Material 2 (mat2)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Solder in the Label text field.

13 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

3 Select Domain 2 only.

4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 50 [GPa] Pa Basic
Poisson's ratio nu 0.4 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 21e-6 1/K Basic

Material 3 (mat3)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Alumina in the Label text field.
3 Select Domain 3 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 300[GPa] Pa Basic
Poisson's ratio nu 0.22 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 8e-6 1/K Basic

DEFINITIONS

Variables 1
1 In the Model Builder window, under Component 1 (comp1) right-click Definitions and
choose Variables.
Define help variables for the user defined strains.
2 In the Settings window for Variables, locate the Variables section.
3 In the table, enter the following settings:

Name Expression Unit

s_mises sqrt([Link]^2+[Link]^2+[Link]^[Link] N/m²
*[Link]*[Link]*[Link]+
3*[Link]^2+3*[Link]^2+3*[Link]^2+(1e-
6[MPa])^2)
alpha_I 3/2*A_I/s_mises*(s_mises/s_ref)^n_I m·s/kg
alpha_II 3/2*A_II/s_mises*(s_mises/s_ref)^n_II m·s/kg

14 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

ADD PHYSICS
1 On the Physics toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Mathematics>ODE and DAE Interfaces>Domain ODEs and
DAEs (dode).
4 Click Add to Component in the window toolbar.

DOMAIN ODES AND DAES (DODE)

1 In the Model Builder window, under Component 1 (comp1) click Domain ODEs and
DAEs (dode).
2 Select Domain 2 only.
3 In the Settings window for Domain ODEs and DAEs, locate the Units section.
4 Find the Source term quantity subsection. In the Unit text field, type 1/s.
5 Click to expand the Dependent variables section. Locate the Dependent Variables
section. In the Field name text field, type ec_I.
6 In the Number of dependent variables text field, type 5.
7 In the Dependent variables table, enter the following settings:

ecx_I
ecy_I
ecz_I
ecxy_I
ece_I

Set the same Gauss point integration order as in the built-in Norton material.
8 In the Model Builder window’s toolbar, click the Show button and select Discretization
in the menu.
9 Click to expand the Discretization section. From the Shape function type list, choose
Gauss point data.
10 From the Element order list, choose 4.

Distributed ODE 1
1 In the Model Builder window, under Component 1 (comp1)>Domain ODEs and DAEs
(dode) click Distributed ODE 1.
2 In the Settings window for Distributed ODE, locate the Source Term section.
3 In the f text-field array, type alpha_I*[Link] on the first row.

15 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

4 In the f text-field array, type alpha_I*[Link] on the 2nd row.

5 In the f text-field array, type alpha_I*[Link] on the 3rd row.
6 In the f text-field array, type alpha_I*[Link] on the 4th row.
7 In the f text-field array, type
2/3*(d(ecx_I,TIME)^2+d(ecy_I,TIME)^2+d(ecz_I,TIME)^2+
2*(d(ecxy_I,TIME)^2))+(1e-20))^0.5
on the 5th row.

SOLID MECHANICS (SOLID)

On the Physics toolbar, click Domain ODEs and DAEs (dode) and choose Solid Mechanics
(solid).

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Recently Used>Domain ODEs and DAEs (dode).
3 Click Add to Component in the window toolbar.

DOMAIN ODES AND DAES 2 (DODE2)

1 In the Model Builder window, under Component 1 (comp1) click Domain ODEs and
DAEs 2 (dode2).
2 Select Domain 2 only.
3 In the Settings window for Domain ODEs and DAEs, locate the Units section.
4 Find the Source term quantity subsection. In the Unit text field, type 1/s.
5 Click to expand the Dependent variables section. Locate the Dependent Variables
section. In the Field name text field, type ec_II.
6 In the Number of dependent variables text field, type 5.
7 In the Dependent variables table, enter the following settings:

ecx_II
ecy_II
ecz_II
ecxy_II
ece_II

8 Click to expand the Discretization section. From the Shape function type list, choose
Gauss point data.

16 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

9 From the Element order list, choose 4.

10 In the Model Builder window, click Distributed ODE 1.
11 In the Settings window for Distributed ODE, locate the Source Term section.
12 In the f text-field array, type alpha_II*[Link] on the first row.
13 In the f text-field array, type alpha_II*[Link] on the 2nd row.
14 In the f text-field array, type alpha_II*[Link] on the 3rd row.
15 In the f text-field array, type alpha_II*[Link] on the 4th row.
16 In the f text-field array, type
2/3*(d(ecx_II,TIME)^2+d(ecy_II,TIME)^2+d(ecz_II,TIME)^2+
2*(d(ecxy_II,TIME)^2))+(1e-20))^0.5
on the 5th row.

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Recently Used>Domain ODEs and DAEs (dode).
3 Click Add to Component in the window toolbar.
4 On the Physics toolbar, click Add Physics to close the Add Physics window.

DOMAIN ODES AND DAES 3 (DODE3)

1 In the Model Builder window, under Component 1 (comp1) click Domain ODEs and
DAEs 3 (dode3).
2 Select Domain 2 only.
3 In the Settings window for Domain ODEs and DAEs, locate the Units section.
4 Find the Dependent variable quantity subsection. From the list, choose Energy density
(J/m^3).
5 Find the Source term quantity subsection. In the Unit text field, type J/(s*m^3).
6 Click to expand the Dependent variables section. Locate the Dependent Variables
section. In the Field name text field, type Wc.
7 In the Number of dependent variables text field, type 2.
8 In the Dependent variables table, enter the following settings:

Wc_I
Wc_II

17 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

9 Click to expand the Discretization section. From the Shape function type list, choose
Gauss point data.
10 From the Element order list, choose 4.
11 In the Model Builder window, click Distributed ODE 1.
12 In the Settings window for Distributed ODE, locate the Source Term section.
13 In the f text-field array, type d(ece_I,TIME)*s_mises on the first row.
14 In the f text-field array, type d(ece_II,TIME)*s_mises on the 2nd row.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies.
4 In the Select study tree, select Preset Studies>Time Dependent.
5 Click Add Study in the window toolbar.
6 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 1

Step 1: Time Dependent

1 In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 From the Time unit list, choose min.
4 In the Times text field, type range(0,0.5,14.5) range(14.6,0.1,15.4)
range(15.5,0.5,29.5) range(29.6,0.1,30.4) range(30.5,0.5,44.5)
range(44.6,0.1,45.4) range(45.5,0.5,60).

Solution 1 (sol1)
1 On the Study toolbar, click Show Default Solver.
2 In the Model Builder window, expand the Solution 1 (sol1) node, then click
Time-Dependent Solver 1.
3 In the Settings window for Time-Dependent Solver, click to expand the Time
stepping section.
4 Locate the Time Stepping section. From the Steps taken by solver list, choose
Intermediate.
5 On the Study toolbar, click Compute.

18 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

RESULTS

Stress (solid)
Click the Zoom Extents button on the Graphics toolbar.

2D Plot Group 2
1 In the Model Builder window, under Results click 2D Plot Group 2.
2 In the Settings window for 2D Plot Group, type Creep strain I (dode) in the
Label text field.

Creep strain I (dode)

1 In the Model Builder window, expand the Results>Creep strain I (dode) node, then
click Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 In the Expression text field, type ece_I.

2D Plot Group 3
1 In the Model Builder window, expand the Results>2D Plot Group 3 node, then click
2D Plot Group 3.
2 In the Settings window for 2D Plot Group, type Creep strain II (dode2 in the
Label text field.

Creep strain II (dode2

1 In the Model Builder window, expand the Results>Creep strain II (dode2 node, then
click Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 In the Expression text field, type ece_II.

2D Plot Group 4
1 In the Model Builder window, under Results click 2D Plot Group 4.
2 In the Settings window for 2D Plot Group, type Dissipated energy (dode3) in
the Label text field.
Create a plot that shows how strain develops during one cycle.

1D Plot Group 5
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Creep strain history in the
Label text field.

19 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Point Graph 1
On the Creep strain history toolbar, click Point Graph.

Creep strain history

1 Select Point 5 only.
2 In the Settings window for Point Graph, locate the y-Axis Data section.
3 In the Expression text field, type solid.ecGp11.
4 Click to expand the Legends section. Select the Show legends check box.
5 From the Legends list, choose Manual.
6 In the table, enter the following settings:

Legends
ec_x

7 Right-click Point Graph 1 and choose Duplicate.

8 In the Settings window for Point Graph, locate the y-Axis Data section.
9 In the Expression text field, type solid.ecGp22.
10 Locate the Legends section. In the table, enter the following settings:

Legends
ec_y

11 Right-click Results>Creep strain history>Point Graph 2 and choose Duplicate.

12 In the Settings window for Point Graph, locate the y-Axis Data section.
13 In the Expression text field, type solid.ecGp33.
14 Locate the Legends section. In the table, enter the following settings:

Legends
ec_z

15 Right-click Results>Creep strain history>Point Graph 3 and choose Duplicate.

16 In the Settings window for Point Graph, locate the y-Axis Data section.
17 In the Expression text field, type solid.ecGp12.
18 Locate the Legends section. In the table, enter the following settings:

Legends
ec_xy

20 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

19 In the Model Builder window, click Creep strain history.

20 In the Settings window for 1D Plot Group, click to expand the Legend section.
21 From the Position list, choose Lower right.
22 Click to expand the Title section. From the Title type list, choose None.
Verify that strains and energies are correctly calculated in the analysis.

1D Plot Group 6
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.

Point Graph 1
On the 1D Plot Group 6 toolbar, click Point Graph.

1D Plot Group 6
1 In the Settings window for Point Graph, locate the y-Axis Data section.
2 In the Expression text field, type [Link].
3 Select Point 5 only.
4 Locate the Legends section. Select the Show legends check box.
5 From the Legends list, choose Manual.
6 In the table, enter the following settings:

Legends
ec (solid)

7 Right-click Point Graph 1 and choose Duplicate.

8 In the Settings window for Point Graph, locate the y-Axis Data section.
9 In the Expression text field, type ece_I.
10 Locate the Legends section. In the table, enter the following settings:

Legends
ece_I (dode)

11 Right-click Results>1D Plot Group 6>Point Graph 2 and choose Duplicate.

12 In the Settings window for Point Graph, locate the y-Axis Data section.
13 In the Expression text field, type ece_II.
14 Locate the Legends section. In the table, enter the following settings:

Legends
ece_II (dode2)

21 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

15 Right-click Results>1D Plot Group 6>Point Graph 3 and choose Duplicate.

16 In the Settings window for Point Graph, locate the y-Axis Data section.
17 In the Expression text field, type ece_I+ece_II.
18 Locate the Legends section. In the table, enter the following settings:

Legends
ece_I+ece_II

19 Click to expand the Coloring and style section. Locate the Coloring and Style section.
Find the Line style subsection. From the Line list, choose Dotted.
20 In the Width text field, type 4.
21 In the Model Builder window, click 1D Plot Group 6.
22 In the Settings window for 1D Plot Group, type Effective creep history in the
Label text field.
23 Locate the Legend section. From the Position list, choose Upper left.
24 Locate the Title section. From the Title type list, choose None.

Effective creep history 1

1 Right-click Effective creep history and choose Duplicate.
2 In the Settings window for 1D Plot Group, type Creep dissipation history in
the Label text field.

Creep dissipation history

1 In the Model Builder window, expand the Results>Creep dissipation history node, then
click Point Graph 1.
2 In the Settings window for Point Graph, locate the y-Axis Data section.
3 In the Expression text field, type [Link].
4 Locate the Legends section. In the table, enter the following settings:

Legends
Wc (solid)

5 In the Model Builder window, under Results>Creep dissipation history click Point
Graph 2.
6 In the Settings window for Point Graph, locate the y-Axis Data section.
7 In the Expression text field, type Wc_I.

22 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

8 Locate the Legends section. In the table, enter the following settings:

Legends
Wc_I (dode3)

9 In the Model Builder window, under Results>Creep dissipation history click Point
Graph 3.
10 In the Settings window for Point Graph, locate the y-Axis Data section.
11 In the Expression text field, type Wc_II.
12 Locate the Legends section. In the table, enter the following settings:

Legends
Wc_II (dode3)

13 In the Model Builder window, under Results>Creep dissipation history click Point
Graph 4.
14 In the Settings window for Point Graph, locate the y-Axis Data section.
15 In the Expression text field, type Wc_I+Wc_II.
16 Locate the Legends section. In the table, enter the following settings:

Legends
Wc_I+Wc_II

17 On the Creep dissipation history toolbar, click Plot.

18 Click the Zoom Extents button on the Graphics toolbar.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

5 Click Add to Component in the window toolbar.

23 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

FATIGUE (FTG)

Strain-Life 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue (ftg)
and choose the domain evaluation Strain-Life.
2 Select Domain 2 only.
3 In the Settings window for Strain-Life, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Coffin-Manson.
5 From the Strain type list, choose User defined.
6 In the εi text field, type 2*ecxy_II.
7 Locate the Fatigue Model Parameters section. From the εf list, choose User defined.
In the associated text field, type 0.587.
8 From the c list, choose User defined. In the associated text field, type -0.61.

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Recently Used>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

4 Click Add to Component in the window toolbar.

5 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 2 (FTG2)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue 2 (ftg2).
2 Select Domain 2 only.

Energy-Based 1
1 Right-click Component 1 (comp1)>Fatigue 2 (ftg2) and choose the domain evaluation
Energy-Based.
2 Select Domain 2 only.
3 In the Settings window for Energy-Based, locate the Fatigue Model Selection section.
4 From the Energy type list, choose User defined.
5 In the Wd text field, type Wc_II.

24 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

6 Locate the Fatigue Model Parameters section. From the Wf’ list, choose User defined.
In the associated text field, type 74e6.
7 From the m list, choose User defined. In the associated text field, type -0.79.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)
Domain ODEs and DAEs (dode)
Domain ODEs and DAEs 2 (dode2)
Domain ODEs and DAEs 3 (dode3)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.
6 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 1, Time Dependent.
6 On the Home toolbar, click Compute.

25 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

26 | A C C E L E R A T E D L I F E TE S T I N G
 Solved with COMSOL Multiphysics 5.2

Bracket—Fatigue Evaluation
Introduction
The S-N curve, also called the Wöhler curve, is one of the most popular methods for
fatigue evaluation. The curve relates stress amplitude to the limiting fatigue life and
can be obtained directly from a set of standard fatigue test. Many times our
applications are subjected to conditions different from the experimental conditions of
the fatigue tests. The fatigue data must then be appropriately modified in order to take
the actual operating conditions into account.

In the example it is shown how to perform a fatigue evaluation when the material data
needs to account for harsh environmental conditions and poor material process. A
tutorial application of a bracket is examined.

Model Definition
The bracket geometry can be seen in Figure 1. The component is subjected to external
loads that lift one arm up and pull the other arm down. The load magnitude is cycled
between a zero load, a peak load and back to zero load. Additional information
regarding the model set up can be found in the documentation of the application
Bracket-Spring Foundation Analysis, found in the Structural Mechanics Module.

Figure 1: Bracket geometry.

The structural steel of the bracket has poor cleanliness and contains fairly large
inclusions. It is well known that with decreasing purity the endurance limit decreases

1 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

and thus large inclusions shorten the lifetime. The fatigue properties of the material
have been tested in a material laboratory and are summarized in Table 1. The material
has a distinct endurance limit at 110 MPa. Tests were preformed in nominal testing
conditions and on specimens with a very good surface finish.
TABLE 1: S-N CURVE DATA

FATIGUE LIFETIME (CYCLES) STRESS AMPLITUDE (MPA)

1.10e3 360
2.17e3 320
3.82e3 290
7.15e3 260
1.45e4 230
3.23e4 200
5.92e4 180
1.16e5 160
2.51e5 140
6.09e5 120
1.00e6 110

As opposite to the testing conditions, the bracket was machined and has a rough
surface as the result. Moreover the bracket operates in salt water that has a strong
influence on the fatigue resistance. Since there is lack of the test data for the operating
conditions the S-N curve data must be modified. One approach is to use

σ a = k ⋅ f SN ( N ) (1)

where σa is the stress amplitude, k is the modification factor, N is the fatigue lifetime,
and fSN is the S-N curve. Based on the experience, the modification factor for the
corrosion in salt water can be set to 0.4. The modification factor for the surface finish,
due to machining manufacturing technique, can be set to 0.7. The total modification
factor in Equation 1 is set to the product of the modification factors for all operating
conditions and thus k = 0.28.

2 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

Results and Discussion

The stress distribution in the bracket a the peak load is shown in Figure 2.

Figure 2: Stress in the bracket at the peak operating load.

A closeup of the part with stress concentrations reveals that also on the inner side of
the bracket there are significant stresses, see Figure 3.

3 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

Figure 3: A closeup of the part of the bracket that experiences highest stresses.

In most of the bracket, the stresses are not high enough to cause fatigue. Only in few
local points is there is a risk of fatigue, see Figure 4. The bracket most probably fails in

4 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

the connection between the two flat parts since the stresses on both the upper and the
lower sides are significantly high to initiate fatigue.

Figure 4: Fatigue lifetime.

In Figure 4 a finite life is predicted close to the holes that fasten the bracket. Those
results are highly dependent on the boundary conditions that were applied there. In
this model a spring foundations was used to fasten the bracket but the constraint could
have also been applied using a fixed connection of with a full contact analysis. This
would change the stress field in the vicinity of the holes somewhat and predict a
different fatigue lifetime.

Notes About the COMSOL Implementation

The S-N curve can be defined using different function types in COMSOL
Multiphysics. When using the interpolation function it is important to use many points
to specify the relation between fatigue life and stress amplitude since the range of the
fatigue life is very large and a small change in the stress amplitude results in a large
change in the fatigue life. The function response between two data points can be
evaluated in different ways but linear interpolation is the most common one. In
Figure 5 the S-N curve is specified with eleven and with four measurement points
respectively. Note that the fatigue life is displayed in a logarithmic scale and therefore

5 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

the interpolation is not a straight line although numerically it is linear. The two curves
4
display large differences. At 200 MPa one curve predicts 3.23 ⋅ 10 cycles while the
4
other one predicts 5.80 ⋅ 10 cycles. This is almost a factor two difference.

Figure 5: S-N curve based different number of measurement points.

Application Library path: Fatigue_Module/Stress_Life/bracket_fatigue

Modeling Instructions
In this example you will start from an existing model from the Structural Mechanics
Module.

From the File menu, choose Open.

Under the Application Library root, browse to the folder

Structural_Mechanics_Module/Tutorials and double-click the file
bracket_spring.mph.

ROOT
Fatigue is calculated based on a load cycle. Create a parameter to scale the magnitude
of the external load.

6 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

GLOBAL DEFINITIONS

Parameters
1 In the Model Builder window, expand the Global Definitions node, then click
Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

para 0 0 Load cycle control

Analytic 1 (load)
1 In the Model Builder window, under Global Definitions click Analytic 1 (load).
2 In the Settings window for Analytic, locate the Definition section.
3 In the Expression text field, type para*F*cos((x+0.3)/5e-2*pi).

COMPONENT 1 (COMP1)
Start by improving the existing mesh in order to get a better resolution of stress
gradients at the critical fillets.

MESH 1
In the Model Builder window, expand the Component 1 (comp1) node.

Free Tetrahedral 2
Right-click Mesh 1 and choose Free Tetrahedral.

Distribution 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Free
Tetrahedral 2 and choose Distribution.
2 Select Edges 28, 29, 97, and 98 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 10.

Size
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Finer.
4 Click the Build All button.

7 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

STUDY 1
Since it is an elastic case and the external load is proportional, only two load cases are
necessary to capture the stress amplitude evaluated in the S-N curve.

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list

para 0 1

6 On the Home toolbar, click Compute.

GLOBAL DEFINITIONS

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Find the Functions subsection. In the table, enter the following settings:

Function name Position in file

wohler 1

5 Click Browse.
6 Browse to the application’s Application Library folder and double-click the file
bracket_fatigue_sn_curve.txt.

7 Click Import.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.

8 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

4 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).

5 Click Add to Component in the window toolbar.
6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)

Stress-Life 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue (ftg)
and choose the boundary evaluation Stress-Life.
2 In the Settings window for Stress-Life, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.
4 Locate the Fatigue Model Selection section. From the σ list, choose Signed von Mises
(principal).
5 From the Modification list, choose Stress factor.
6 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
7 Locate the Fatigue Model Parameters section. From the fSN(N) list, choose wohler.
8 In the k text field, type 0.28.
9 Locate the Evaluation Settings section. In the Ncut text field, type 1e6.

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

9 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 1, Stationary.
6 On the Home toolbar, click Compute.

RESULTS

Cycles to Failure (ftg)

Rotate and zoom the bracket to get a better view of the critical points.

10 | BRACKET—FATIGUE EVALUATION
 Solved with COMSOL Multiphysics 5.2

Cycle Counting in Fatigue Analysis -

Benchmark
Introduction
In this benchmark model for the Rainflow counting algorithm, values computed by the
Fatigue Module in COMSOL are compared with the ASTM standard E1049-85 (Ref.
1). An extension to the benchmark compares the Palmgren-Miner cumulative damage
model to analytical expressions.

Model Definition
A flat test specimen is subjected to a repeated load cycle. The material has elastic
properties defined by Young’s modulus E = 69 GPa and Poisson’s ratio ν = 0.34. The
thickness of the specimen is 6.25 mm. The remaining dimensions are given in
Figure 1.

200

50 12.5
12.5

20

Figure 1: Flat test specimen. Dimensions are given in mm.

The Wöhler diagram (S-N curve) is given by the expression

2.35 – 0.091
σ a = 94e6 ( R ⁄ ( – 0.36 ) ) N

where σa is the stress amplitude, R is the R-value and N is the number of cycles to
failure for a constant stress cycle defined by σa and R. The relation is valid for
– 2 .5 ≤ R ≤ – 0 .2 , while the parameter N ≥ 1e8 is seen as infinite life and it should not
be taken into account in damage calculations.

1 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

An ASTM cycle (Ref. 1) is evaluated in the example. The load history of the cycle is
presented in Table 1.
TABLE 1: FATIGUE CYCLE

STEP LOAD UNITS

1 -2
2 1
3 -3
4 5
5 -1
6 3
7 -4
8 4
9 2

The unit load can be chosen arbitrarily and is here selected so that one unit load
corresponds to 10 MPa stress in the central cross section.

The ASTM example is further extended to examine the fatigue damage caused by
100000 blocks of the cycle.

In the test specimen, the stress state in the central part away from the fillets does not
vary with the position. Therefore, any point in the central thin cross section can be
chosen for evaluation.

Results and Discussion

Because of the symmetry, only a quarter of the test specimen is modeled. The 2D Solid
Mechanics interface is used with a plane stress assumption. Figure 2 shows the axial
stress caused by a unit load. Although a quarter of the specimen was modeled, the
results for the whole specimen can be examined using Data Sets of type Mirror 2D.

2 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Figure 2: Axial stress in the specimen caused by a unit load.

The Cumulative Damage feature generates results of the cycle counting algorithm as
well as a damage calculation. A point in the thin central cross section is evaluated. The
applied load cycle, see Table 1, is quantified with the Rainflow Counting algorithm
and shown in Figure 3. The ASTM results (Ref. 1) are shown in Table 2
TABLE 2: ASTM RAINFLOW COUNTING RESULTS

σa (MPA) σm (MPA) n
20 -10 0.5
15 -5 0.5
40 0 0.5
45 5 0.5
40 10 0.5
30 10 0.5
1 10 1.0

3 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

here, σm is the cycle mean stress amplitude and n is number of cycles.

Figure 3: Load cycle quantified with the Rainflow Counting option.

The Rainflow Counting results by ASTM and COMSOL are in perfect agreement.

The evaluation of the cumulative damage is now compared against analytical

expressions. In the Palmgren-Miner model, the damage is calculated for each stress
bin. The number of cycles to failure for a constant cycle is taken from the S-N curve
which is evaluated at the center of the bin.

With the chosen bin discretization the damage is evaluated at following bin stress
centers:

b
σ a = 17.1 MPa, 21.4 MPa, 25.7 MPa, 30.0 MPa, 34.3 MPa, 38.6 MPa,
42.9 MPa

and

b
σ m = -8.0 MPa, -4.0 MPa, 0.0 MPa, 4.0 MPa, 8.0 MPa

4 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

The key values for the damaging bins are presented in Table 3. The superscript ‘b’
denotes that the variable is evaluated at the bin center.
TABLE 3: DAMAGING BIN DATA

σa (MPA) σm (MPA) σab (MPA) σmb (MPA) Rb nb Nb

20.0 -10.0 21.4 -8.0 -2.19 0.500 Inf
15.0 -5.0 17.1 -4.0 -1.61 0.500 Inf
40.0 0.0 38.6 0.0 -1.00 0.500 3.44e7
45.0 5.0 42.9 4.0 -0.829 0.500 2.31e6
40.0 10.0 38.6 8.0 -0.656 0.500 5.84e5
30.0 10.0 30.0 8.0 -0.579 0.500 1.45e6
20.0 10.0 21.4 8.0 -0.456 1.00 2.47e6

The fatigue usage factor fus following the Palmgren-Miner linear damage rule is
calculated using the expression

 ni ⁄ Ni
b b
f us = m (1)
i=1

where m is the number or repeatable block and p is the number of bins.

ASTM
Following Equation 1 the fatigue usage factor of the ASTM cycles is f us = 0.184 .
The COMSOL result is fus= 0.182. The small discrepancy between the results is
attributed to the evaluation of the S-N curve.

5 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

The relative damage contribution from each bin is shown in Figure 4.

Figure 4: Contribution of each stress bin to the fatigue usage.

The Cumulative Damage model is based on the Palmgren-Miner linear damage rule.
It is a discrete model in the sense that the calculations are based on stress bins which
holds all stress cycles within a certain stress amplitude and mean stress range.

In this example, the stress amplitude range in a bin is 20 MPa/5=4 MPa and the mean
stress range in a bin is 30 MPa/7 MPa = 4.3 MPa. By changing the bin discretization,
the graph showing counted stress cycles in Figure 3, appears with a different
resolution.

As for the results of the relative damage in Figure 4, a change in discretization can
affect the results significantly. The damage is evaluated based on the bin stresses and
not true cycle stresses. Consider the cycle defined by σa = 45.0 MPa and
b b
σm = 5.0 MPa. It is evaluated in a bin defined by σ a = 42.9 MPa and σ m = 4.0 MPa.
Since bin stresses are lower than true stresses in that cycle, they predict less damage and
thus give a nonconservative contribution to the fatigue usage factor.

6 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Notes About the COMSOL Implementation

In the Cumulative Damaged feature, the arguments of the S-N curve must be specified
by first giving the R-value followed by the number of cycles, as it is done in the Analytic
function in the Arguments field.

Reference
1. ASTM International, Standard Practices for Cycle Counting in Fatigue Analysis,
Designation: E1049-85 (Reapproved 2011).

Application Library path: Fatigue_Module/Verification_Examples/

cycle_counting_benchmark

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.

7 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 10.
4 In the Height text field, type 100.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 10-6.25.
4 In the Height text field, type 50-sqrt(12.5^2-8.75^2).
5 Locate the Position section. In the x text field, type 6.25.

Circle 1 (c1)
1 On the Geometry toolbar, click Primitives and choose Circle.
2 In the Settings window for Circle, locate the Position section.
3 In the x text field, type 6.25+12.5.
4 In the y text field, type 50-sqrt(12.5^2-8.75^2).
5 Locate the Size and Shape section. In the Radius text field, type 12.5.

Difference 1 (dif1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 In the Settings window for Difference, locate the Difference section.
3 In the Relative repair tolerance text field, type 1e-6.
4 From the bigger rectangle subtract the smaller rectangle and the circle.
5 Click the Build All Objects button.

GLOBAL DEFINITIONS
Specify the load cycle.

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 In the table, enter the following settings:

t f(t)
1 -2
2 1

8 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

t f(t)
3 -3
4 5
5 -1
6 3
7 -4
8 4
9 -2

SOLID MECHANICS (SOLID)

1 In the Model Builder window, under Component 1 (comp1) click Solid Mechanics
(solid).
2 In the Settings window for Solid Mechanics, locate the 2D Approximation section.
3 From the list, choose Plane stress.
4 Locate the Thickness section. In the d text field, type 0.00625.

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundaries 1 and 2 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 3 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 From the Load type list, choose Total force.
5 Specify the Ftot vector as

0 x
F*int1(case) y

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.

9 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

3 In the table, enter the following settings:

Name Expression Value Description

F 10*6.25*12.5/2 390.6 Load unit
case 1 1 Load case

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 69e9 Pa Basic
Poisson's ratio nu 0.34 1 Basic
Density rho 2700 kg/m³ Basic

STUDY 1

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list

case range(1,1,9)

6 On the Home toolbar, click Compute.

RESULTS
Mirror solution of a quarter of a specimen and create results for a full specimen.

Mirror 2D 1
On the Results toolbar, click More Data Sets and choose Mirror 2D.

10 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Data Sets
1 In the Settings window for Mirror 2D, locate the Axis Data section.
2 In row Point 2, set Y to 100.

Mirror 2D 2
On the Results toolbar, click More Data Sets and choose Mirror 2D.

Data Sets
1 In the Settings window for Mirror 2D, locate the Data section.
2 From the Data set list, choose Mirror 2D 1.
3 Locate the Axis Data section. In row Point 1, set x to -6.25.
4 In row Point 2, set x to 6.25.
5 In row Point 2, set y to 0.
Display stress state in the whole specimen.

Stress (solid)
1 In the Model Builder window, under Results click Stress (solid).
2 In the Settings window for 2D Plot Group, locate the Data section.
3 From the Data set list, choose Mirror 2D 2.
4 From the Parameter value (case) list, choose 2.
5 Click to expand the Title section. From the Title type list, choose None.
6 In the Model Builder window, expand the Stress (solid) node, then click Surface 1.
7 In the Settings window for Surface, click Replace Expression in the upper-right corner
of the Expression section. From the menu, choose Component 1>Solid
Mechanics>Stress (Gauss points)>Second Piola-Kirchhoff stress, Gauss-point evaluation
(Material)>[Link] - Second Piola-Kirchhoff stress, Gauss-point evaluation, Y
component.
8 On the Stress (solid) toolbar, click Plot.

1D Plot Group 2
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Load Cycle Response in the Label
text field.

Point Graph 1
On the Load Cycle Response toolbar, click Point Graph.

11 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Load Cycle Response

1 Select Point 3 only.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid
Mechanics>Stress (Gauss points)>Second Piola-Kirchhoff stress, Gauss-point evaluation
(Material)>[Link] - Second Piola-Kirchhoff stress, Gauss-point evaluation, Y
component.
3 On the Load Cycle Response toolbar, click Plot.

GLOBAL DEFINITIONS
Specify the load cycle.

Analytic 1 (an1)
1 On the Home toolbar, click Functions and choose Global>Analytic.
2 In the Settings window for Analytic, locate the Definition section.
3 In the Arguments text field, type R, N.
4 In the Expression text field, type (94e6*(R/-0.36)^1.15)*N^-0.119.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)
In the Model Builder window, expand the Study 1 node.

12 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Cumulative Damage 1
1 Right-click Component 1 (comp1)>Fatigue (ftg) and choose Points>Cumulative
Damage.
Any point in the thin section away from the notch will give the same fatigue
response.
2 Select Point 3 only.
3 In the Settings window for Cumulative Damage, locate the Solution Field section.
4 From the Physics interface list, choose Solid Mechanics (solid).
5 Locate the Cycle Counting Parameters section. Find the Discretization subsection. In
the Nm text field, type 5.
6 In the Nr text field, type 7.
7 Locate the Damage Model Parameters section. From the σa(R,N) list, choose an1.
8 In the m text field, type 100000.
Set cutoff which can be seen as a limit for infinite life that does not contribute to the
damage.
9 Find the Evaluation settings subsection. In the Ncut text field, type 1e8.

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.

13 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 1, Stationary.
6 On the Home toolbar, click Compute.

14 | CYCLE COUNTING IN FATIGUE ANALYSIS - BENCHMARK

 Solved with COMSOL Multiphysics 5.2

Notch Approximation to Low-Cycle

Fatigue Analysis of Cylinder with a
Hole
Introduction
A load carrying component of a structure is subjected to multi-axial cyclic loading
during which a localized yielding of the material occurs. In this model you perform a
low cycle fatigue analysis of the part based on the Smith-Watson-Topper (SWT)
model. Due to localized yielding, you can use two methods to obtain the stress and
strain distributions for the fatigue evaluation. The first method is an elastoplastic
analysis with linear kinematic hardening, while the second is a linear elastic analysis
with Neuber correction for plasticity, based on the Ramberg-Osgood model. This
example explores the second method. In the model Elastoplastic Low-Cycle Fatigue
Analysis of Cylinder with a Hole, the same problem is solved using the full elastoplastic
approach.

Model Definition

GEOMETRY
A cylinder contains a hole, drilled perpendicularly to its axis. The outer and inner
diameters of the cylinder are 200 and 180 mm respectively. Its height is 100 mm. The
diameter of the hole is 20 mm. The cylinder is loaded by an axial force which varies in
time.

As the structure and loading contains several symmetries, you may model only 1/8 of
the cylinder, a shown in Figure 1.

1 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Symmetry planes

Applied load

Figure 1: Model geometry with constrained and loaded faces.

MATERIAL PROPERTIES
• Elastic data: Isotropic with E = 210 GPa, ν = 0.3
• Cyclic Ramberg-Osgood plasticity data: K’ = 1550 MPa, n’ = 0.16
• Fatigue parameters for the SWT equation:

- σf' = 1323 MPa

- b = •0.097
- εf' = 0.375
- c = •0.60

CONSTRAINTS
Apply symmetry conditions on the three symmetry sections shown in Figure 1.

LOAD
The loaded boundary of the cylinder is subjected to a pressure varying between
+200 MPa and -200 MPa.

2 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Results and Discussion

The von Mises stress distribution at maximum load is shown in Figure 2. Notice that
the yield limit(380 = MPa) is extensively exceeded, so in reality significant plastic
strains can be expected.

Figure 2: von Mises stress level at the first maximum load

The computed number of cycles to fatigue is shown in Figure 3. It is slightly below

10000 cycles.

3 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Figure 3: The distribution of expected number of cycles at the hole.

Notes About the COMSOL Implementation

When using the elastic approach to low cycle fatigue, it is only the peak loads which
are important, not the load cycle as such. Include always both the maximum and the
minimum load, since some of the criteria are sensitive to the sign of the stresses. It is
also necessary to include the ‘zero’ load case, which contains no loads.

Application Library path: Fatigue_Module/Strain_Based/

cylinder_with_hole_elastic

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

4 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Cylinder 1 (cyl1)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 100.
4 In the Height text field, type 100.

Cylinder 2 (cyl2)
1 Right-click Cylinder 1 (cyl1) and choose Duplicate.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 90.

Difference 1 (dif1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 Select the object cyl1 only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Active toggle button.
5 Select the object cyl2 only.
6 Click the Build All Objects button.

Cylinder 3 (cyl3)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 10.

5 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

4 In the Height text field, type 220.

5 Locate the Position section. In the y text field, type -110.
6 In the z text field, type 50.
7 Locate the Axis section. From the Axis type list, choose y-axis.
8 Click the Build All Objects button.

Difference 2 (dif2)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 Select the object dif1 only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Active toggle button.
5 Select the object cyl3 only.
6 Click the Build All Objects button.

Block 1 (blk1)
1 On the Geometry toolbar, click Block.
2 In the Settings window for Block, locate the Size and Shape section.
3 In the Width text field, type 100.
4 In the Depth text field, type 100.
5 In the Height text field, type 50.
6 Click the Build All Objects button.

Intersection 1 (int1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Intersection.
2 Select the objects dif2 and blk1 only.
3 Click the Build All Objects button.
4 Click the Zoom Extents button on the Graphics toolbar.

SOLID MECHANICS (SOLID)

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundaries 1, 6, and 7 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.

6 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

2 Select Boundary 3 only.

3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
0 y
-200[MPa] z

5 On the Physics toolbar, click Load Group and choose New Load Group.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 210[GPa] Pa Basic
Poisson's ratio nu 0.3 1 Basic
Density rho 0 kg/m³ Basic

MESH 1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Build All.

STUDY 1

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Define load cases check box.
In an analysis of this type you should always supply the maximum and the minimum
load, and also a zero solution.
4 Click Add three times.

7 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

5 In the table, enter the following settings:

6 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
1 In the Model Builder window, expand the Results>Stress (solid) node, then click
Surface 1.
2 In the Settings window for Surface, locate the Expression section.
3 From the Unit list, choose MPa.
4 On the Stress (solid) toolbar, click Plot.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)

Strain-Based 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue (ftg)
and choose the boundary evaluation Strain-Based.

8 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

2 Select Boundary 4 only.

As the elastic approach to strain based fatigue is based on a notch assumption, it only
makes sense to select boundaries at the hole.
3 In the Settings window for Strain-Based, locate the Fatigue Model Selection section.
4 From the Solution type list, choose Elastic solution with notch assumption.
5 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

MATERIALS

Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose All boundaries.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Fatigue ductility coefficient epsilonf_CM 0.375 1 Coffin-Manson
Fatigue ductility exponent c_CM -0.60 1 Coffin-Manson
Fatigue strength coefficient sigmaf_Basquin 1323[MPa] Pa Basquin
Fatigue strength exponent b_Basquin -0.097 1 Basquin
Young's modulus E solid.E Pa Basic
Poisson's ratio nu [Link] 1 Basic
Cyclic hardening coefficient K_ROcyclic 1550[MPa] Pa Ramberg-Osgood
Cyclic hardening exponent n_ROcyclic 0.16 1 Ramberg-Osgood

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.

9 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 1, Stationary.
6 On the Home toolbar, click Compute.

10 | NOTCH APPROXIMATION TO LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Elastoplastic Low-Cycle Fatigue

Analysis of Cylinder with a Hole
Introduction
A load carrying component of a structure is subjected to multi-axial cyclic loading
during which localized yielding of the material occurs. In this model you perform a
low cycle fatigue analysis of the part based on the Smith-Watson-Topper (SWT)
model. Due to localized yielding, you can use two methods to obtain the stress and
strain distributions for the fatigue evaluation. The first method is a full elastoplastic
analysis with kinematic hardening, while the second is a linear elastic analysis with
Neuber correction for plasticity. This example explores the first method. In the model
Notch Approximation to Low-Cycle Fatigue Analysis of Cylinder with a Hole, the
same problem is solved using the elastic approach.

Model Definition

GEOMETRY
A cylinder contains a hole, drilled perpendicularly to its axis. The outer and inner
diameters of the cylinder are 200 and 180 mm respectively. Its height is 100 mm. The
diameter of the hole is 20 mm. The cylinder is loaded by an axial force which varies in
time.

As the structure and loading contains several symmetries, you may model only 1/8 of
the cylinder, a shown in Figure 1.

1 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Symmetry planes

Applied load

Figure 1: Model geometry with constrained and loaded faces.

MATERIAL PROPERTIES
• Elastic data: Isotropic with E = 210 GPa, ν = 0.3
• Kinematic hardening plasticity data: Yield stress 380 MPa, Tangent modulus
75 GPa
• Fatigue parameters for the SWT equation:

- σf' = 1323 MPa

- b = •0.097
- εf' = 0.375
- c = •0.60

CONSTRAINTS
Apply symmetry conditions on the three symmetry sections shown in Figure 1.

LOAD
The loaded boundary of the cylinder is subjected to a pressure of 200 MPa having a
sinusoidal variation with time. Since the problem is quasi-static, time is not used
explicitly in the problem. Instead you model the load as a function of a parameter, and
use the parametric solver to trace the history.

2 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Results and Discussion

The von Mises stress when the maximum load is reached for the first time (at a
parametric value of 0.25) is shown in Figure 2. Notice that the yield limit 380 MPa is
exceeded by a large factor, and you can therefore expect significant plastic strains.

After completion of the first load cycle, there are no external loads. The internal
residual stresses are still high, as can be seen in Figure 3.

For a point located on the face of the hole and in the XY symmetry plane, the highest
stresses appear in the Z direction. In Figure 4, you can follow the development of the
normal stress and strain in this direction during the first load cycle. The accumulation
of plastic strain is shown in Figure 5. This repetitive plastic deformation can be viewed
as driving the generation of a fatigue crack.

After having run a second load cycle, the stress-strain loop repeats itself, and you can
consider the process to have reached a steady state condition. This is shown in
Figure 6. The continuation of the effective strain history is shown in Figure 7.

Figure 2: von Mises stress level at the first maximum load

3 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Figure 3: Residual effective stress after completion of first load cycle

Figure 4: Stress as function of strain during the first load cycle.

4 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Figure 5: Accumulated effective plastic strain during the first load cycle

Figure 6: Stress as function of strain after two load cycles

5 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Figure 7: Accumulated effective plastic strain during two load cycles

The result of the fatigue evaluation is shown in Figure 8 below. The most critical point
has a computed life of approximately 13000 load cycles.

6 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Figure 8: Lifetime plot

Notes About the COMSOL Implementation

You need to use two studies for the elastoplastic analysis. The second study contains
the strain cycle that is passed on to the fatigue evaluation. It takes the results from the
first study as initial conditions.

Application Library path: Fatigue_Module/Strain_Based/

cylinder_with_hole_plastic

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

7 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Cylinder 1 (cyl1)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 100.
4 In the Height text field, type 100.

Cylinder 2 (cyl2)
1 Right-click Cylinder 1 (cyl1) and choose Duplicate.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 90.

Difference 1 (dif1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 Select the object cyl1 only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Active toggle button.
5 Select the object cyl2 only.
6 Click the Build All Objects button.

Cylinder 3 (cyl3)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 10.

8 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

4 In the Height text field, type 220.

5 Locate the Position section. In the y text field, type -110.
6 In the z text field, type 50.
7 Locate the Axis section. From the Axis type list, choose y-axis.
8 Click the Build All Objects button.

Difference 2 (dif2)
1 On the Geometry toolbar, click Booleans and Partitions and choose Difference.
2 Select the object dif1 only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Active toggle button.
5 Select the object cyl3 only.
6 Click the Build All Objects button.

Block 1 (blk1)
1 On the Geometry toolbar, click Block.
2 In the Settings window for Block, locate the Size and Shape section.
3 In the Width text field, type 100.
4 In the Depth text field, type 100.
5 In the Height text field, type 50.
6 Click the Build All Objects button.

Intersection 1 (int1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Intersection.
2 Select the objects dif2 and blk1 only.
3 Click the Build All Objects button.
4 Click the Zoom Extents button on the Graphics toolbar.
Use an extra domain for controlling the mesh. The mesh should be fine in a region
close to the hole.

Cylinder 4 (cyl4)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Axis section.
3 From the Axis type list, choose y-axis.
4 Locate the Size and Shape section. In the Height text field, type 40.

9 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

5 In the Radius text field, type 15.

6 Locate the Position section. In the y text field, type 80.
7 In the z text field, type 50.
8 Click the Build All Objects button.

Intersection 2 (int2)
1 On the Geometry toolbar, click Booleans and Partitions and choose Intersection.
2 Select the objects cyl4 and int1 only.
3 In the Settings window for Intersection, locate the Intersection section.
4 Select the Keep input objects check box.
5 Click the Build All Objects button.

Delete Entities 1 (del1)

1 In the Model Builder window, right-click Geometry 1 and choose Delete Entities.
2 In the Settings window for Delete Entities, locate the Entities or Objects to Delete
section.
3 From the Geometric entity level list, choose Object.
4 Select the object cyl4 only.
5 Click the Build All Objects button.

Mesh Control Domains 1 (mcd1)

1 On the Geometry toolbar, click Virtual Operations and choose Mesh Control Domains.
2 On the object fin, select Domain 2 only.
3 On the Geometry toolbar, click Build All.

SOLID MECHANICS (SOLID)

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundaries 1, 6, and 7 only.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.

10 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

3 In the table, enter the following settings:

Name Expression Value Description

para 0 0 Load cycle control

SOLID MECHANICS (SOLID)

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 3 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 Specify the FA vector as

0 x
0 y
-200[MPa]*sin(2*pi*para) z

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

Plasticity 1
1 On the Physics toolbar, click Attributes and choose Plasticity.
2 In the Settings window for Plasticity, locate the Plasticity Model section.
3 From the Hardening model list, choose Kinematic.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 210[GPa] Pa Basic
Poisson's ratio nu 0.3 1 Basic
Density rho 0 kg/m³ Basic

11 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

Property Name Value Unit Property group

Initial yield stress sigmags 380[MPa] Pa Elastoplastic material
model
Kinematic tangent modulus Ek 75[GPa] Pa Elastoplastic material
model

MESH 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and
choose Build All.
The default mesh cannot resolve this problem, so you need to refine it.
Once a mesh is created the mesh control domains are removed. In order to get them
back, you need to clear the mesh.
2 On the Mesh toolbar, click Clear Mesh.
3 In the Model Builder window, click Mesh 1.
4 In the Settings window for Mesh, locate the Mesh Settings section.
5 From the Sequence type list, choose User-controlled mesh.

Size 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Free
Tetrahedral 1 and choose Size.
2 In the Settings window for Size, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 2 only.
5 Locate the Element Size section. Click the Custom button.
6 Locate the Element Size Parameters section. Select the Maximum element size check
box.
7 In the associated text field, type 1.5.
8 Click the Build All button.

Size
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Fine.
4 Click the Build All button.

12 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

STUDY 1

Step 1: Stationary
1 In the Model Builder window, under Study 1 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add to include para as Continuation parameter.
5 In the table, enter the following settings:

Parameter name Parameter value list

para range(0,0.025,1)

6 In the Model Builder window, click Study 1.

7 In the Settings window for Study, type Study 1 (First load cycle) in the Label
text field.
8 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
1 On the Stress (solid) toolbar, click Plot.
The default plot is from the last parameter value, which has no external loads. The
residual stress and corresponding displacements are shown. You may want to reduce
the deformation scale in order to get a better view.
2 In the Model Builder window, expand the Stress (solid) node.
3 In the Model Builder window, expand the Results>Stress (solid)>Surface 1 node, then
click Deformation.
4 In the Settings window for Deformation, locate the Scale section.
5 Select the Scale factor check box.
6 In the associated text field, type 2000.
7 On the Stress (solid) toolbar, click Plot.
Now, set back scaling to the default and examine the stress state at maximum
tension.
8 Clear the Scale factor check box.
9 In the Model Builder window, click Stress (solid).
10 In the Settings window for 3D Plot Group, locate the Data section.

13 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

11 From the Parameter value (para) list, choose 0.25.

12 On the Stress (solid) toolbar, click Plot.

1D Plot Group 2
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.

Point Graph 1
On the 1D Plot Group 2 toolbar, click Point Graph.

1D Plot Group 2
1 Select Point 6 only.
2 In the Settings window for Point Graph, locate the y-Axis Data section.
3 In the Expression text field, type [Link].
4 On the 1D Plot Group 2 toolbar, click Plot.

1D Plot Group 3
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.

Point Graph 1
On the 1D Plot Group 3 toolbar, click Point Graph.

1D Plot Group 3
1 Select Point 6 only.
2 In the Settings window for Point Graph, locate the y-Axis Data section.
3 In the Expression text field, type [Link].
4 Locate the x-Axis Data section. From the Parameter list, choose Expression.
5 In the Expression text field, type [Link].
6 On the 1D Plot Group 3 toolbar, click Plot.
Add another study for the next load cycle.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies>Stationary.
4 Click Add Study in the window toolbar.
5 On the Home toolbar, click Add Study to close the Add Study window.

14 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

STUDY 2

Step 1: Stationary
1 In the Model Builder window, under Study 2 click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Study Extensions section.
3 Select the Auxiliary sweep check box.
4 Click Add to include para as Continuation parameter.
5 In the table, enter the following settings:

Parameter name Parameter value list

para range(1,0.025,2)

6 Click to expand the Values of dependent variables section. Locate the Values of
Dependent Variables section. Find the Initial values of variables solved for subsection.
From the Settings list, choose User controlled.
7 From the Method list, choose Solution.
8 From the Study list, choose Study 1 (First load cycle), Stationary.
9 In the Model Builder window, click Study 2.
10 In the Settings window for Study, type Study 2 (Steady state load cycle) in
the Label text field.
11 On the Home toolbar, click Compute.

RESULTS

1D Plot Group 2
1 In the Model Builder window, expand the Results>1D Plot Group 2 node.
2 Right-click Point Graph 1 and choose Duplicate.
3 In the Settings window for Point Graph, locate the Data section.
4 From the Data set list, choose Study 2 (Steady state load cycle)/Solution 2 (sol2).
5 Click to expand the Coloring and style section. Locate the Coloring and Style section.
Find the Line style subsection. From the Line list, choose Dotted.
6 In the Width text field, type 4.
7 On the 1D Plot Group 2 toolbar, click Plot.

1D Plot Group 3
1 In the Model Builder window, expand the Results>1D Plot Group 3 node.
2 Right-click Point Graph 1 and choose Duplicate.

15 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

3 In the Settings window for Point Graph, locate the Data section.
4 From the Data set list, choose Study 2 (Steady state load cycle)/Solution 2 (sol2).
5 Locate the Coloring and Style section. Find the Line style subsection. From the Line
list, choose Dotted.
6 In the Width text field, type 4.
7 On the 1D Plot Group 3 toolbar, click Plot.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1 (First load cycle)
Study 2 (Steady state load cycle)

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)

Strain-Based 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue (ftg)
and choose the boundary evaluation Strain-Based.
2 In the Settings window for Strain-Based, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.
4 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
5 Locate the Evaluation Settings section. Find the Critical plane settings subsection. In
the Q text field, type 16.

16 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

MATERIALS

Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose All boundaries.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Fatigue ductility coefficient epsilonf_CM 0.375 1 Coffin-Manson
Fatigue ductility exponent c_CM -0.60 1 Coffin-Manson
Fatigue strength coefficient sigmaf_Basquin 1323[MPa] Pa Basquin
Fatigue strength exponent b_Basquin -0.097 1 Basquin
Young's modulus E 210[GPa] Pa Basic

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 3

Step 1: Fatigue
1 In the Model Builder window, under Study 3 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.

17 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 2 (Steady state load cycle), Stationary.
6 In the Model Builder window, click Study 3.
7 In the Settings window for Study, type Study 3 (Fatigue, SWT) in the Label text
field.
8 On the Home toolbar, click Compute.

RESULTS

Cycles to Failure (ftg)

1 In the Settings window for 3D Plot Group, type Fatigue life in the Label text
field.
2 Click to expand the Title section. From the Title type list, choose Manual.
3 In the Title text area, type Logarithm of lifetime in number of cycles.
4 On the Fatigue life toolbar, click Plot.

18 | ELASTOPLASTIC LOW-CYCLE FATIGUE ANALYSIS OF CYLINDER WITH A HOLE

 Solved with COMSOL Multiphysics 5.2

High-Cycle Fatigue Analysis of a

Cylindrical Test Specimen
Introduction
A benchmark model for the Fatigue Module. A cylindrical test is subjected to
non-proportional loading. Three stress based models: Findley, Matake, Normal stress,
are compared to analytical values and to each other. The non-smooth behavior of the
Matake model is captured and discussed.

This application presents verifications of the stress based fatigue models in COMSOL
using an example of a circular specimen subjected to normal and torsional loading.

Model Definition
A cylindrical test specimen with a geometry according to Figure 1 is subjected to
fatigue testing. Both normal load and twisting moment are applied in such a way that
the stress scenario shown in Figure 2 is affecting in the central thin part of the
specimen. The maximum shear stress is experienced on the outer radius of the bar.
50
20
r 50
r5 20

10
5

Figure 1: Geometry of the test specimen.

1 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

τ[MPA]

200

σ[MPA]
200

Figure 2: Stress history at the outer radius in the central part of the specimen.

This stress state can be simulated with cyclic normal force of 15708 N and cyclic
twisting moment of 22.672 Nm.

The specimen is made of mild steel with Young’s modulus E of 210 GPa and a
Poisson’s ratio ν of 0.3.

The examined fatigue criteria are Findley, Matake, and Normal stress. Material
parameters are calculated from uniaxial fatigue tests. Results of two tests are available.
In the first one material is tested in reversed tension-compression, which means that
the load amplitude is alternating around zero stress, and the endurance limit is
σ R = – 1 = 350 MPa. In the second one material is tested in pure tension, where load
is pulsating between zero to two times the load amplitude, and the endurance limit is
σ R = 0 = 288 MPa. The denominator of the endurance limit shows the R-value of the
test.

Findley parameters are related to the above fatigue data via

σR = –1 2
f F = ------------------ ( k F + 1 + k F )
2
σR = 0 2
f F = --------------- ( 2k F + 1 + 4k F )
2

and are fF = 213 MPa and kF = 0.20. Matake parameters are related to the above
fatigue data via

2 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

σR = –1
f M = ------------------ ( 1 + k M )
2
σR = 0
f M = --------------- ( 1 + 2k M )
2

and are fM = 223 MPa and kF = 0.27. In the normal stress model the stress limit equals
twice the stress amplitude of the fatigue test. The model does not take into account
the R-value dependence and therefor two different limits are obtained. In order to give
conservative results the lower limit, σ R = 0 ,is chosen and thus the model parameter is
fN = 576 MPa.

Results and Discussion

The load cycle is obtained in four load steps. The resulting stress at the outer radius of
the thin section follows the specification, see Figure 3.

Figure 3: Resulting stresses at the outer radius of the smallest section of the specimen.

The three fatigue criteria evaluated on the boundary of the specimen are presented in
Figure 4, Figure 5, and Figure 6.

3 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

In this example, it is the central part of the specimen which is in focus for the
evaluation, even though the maximum fatigue usage is slightly higher in the fillets.

The Findley model shows that the highest fatigue usage factor is found at the transition
between the thin cylinder and the fillet, with the value of 0.99 against 0.95 in the
center. A smooth transition in the results indicates that the specified search resolution
for the critical plane is sufficient to correctly capture the fatigue response.

Figure 4: Fatigue usage factor Findley criterion.

The Matake criterion gives less smooth results. The reason can be found in the
definition of the criterion, which selects the critical plane using the shear stresses only
and then adds the normal stress. This is discussed in more detail on page 9. The analysis
of the test specimen predicts a fatigue usage factor in the range from 0.78 to 0.98 in

4 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

the central part, as compared to 1.03 at some points at the transition from the thin
cylinder to the notch.

Figure 5: Fatigue usage factor Matake criterion.

The Normal Stress criterion predicts yet different fatigue usage factor than the two
previous models. Here again the highest fatigue usage factor is found on the transition

5 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

from the thin section to the notch, 0.93, as compared to the center of the specimen,
where it is 0.88.

Figure 6: Fatigue usage factor Normal stress criterion.

The stress state on the surface in the central part of the specimen can be evaluated
using analytical expressions. As there are no tractions on the boundary, it is in a state
of plane stress. Therefore, the search for the critical plane can be reduced to two
dimensions. In that case, the normal and shear stresses on a plane can be obtained as

2 2
σ ( ϕ ) = σ a cos ϕ + σ b sin ϕ + τ ab sin ϕ cos ϕ
σb – σa (1)
τ ( ϕ ) = ------------------ sin 2ϕ + τ ab cos 2ϕ
2

where σa and σb are orthogonal normal stresses, and τab is the shear stress. The
orientation angle ϕ is calculated from the axis of σa towards σb. In the current
example, the σa is the axial stress, and τab as the shear stress caused by the twisting
moment.

In the current example, a search resolution of Q = 11 is used. This gives a 9° angular

increment in the critical plane evaluation. Using Equation 1, the normal and shear
stresses evaluated on discrete planes are given in Table 1 and Table 2. The number in

6 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

the subscripts indicates the load case.

TABLE 1: NORMAL STRESS ON DISCRETE PLANES

ϕ (°) σn,1(MPa) σn,2(MPa) σn,3(MPa) σn,4(MPa) σn,max(MPa) Δσn(MPa)

0 200 200 -200 -200 200 400
9 231 159 -231 -159 231 462
18 249 113 -249 -133 249 498
27 252 65 -252 -65 252 504
36 241 21 -241 -21 241 481
45 215 -15 -215 15 215 431
54 179 -41 -17 41 179 358
63 135 -52 -135 52 135 269
72 87 -49 -87 49 87 174
81 41 -31 -41 31 41 81
90 0 0 0 0 0 0
99 -31 41 31 -41 41 81
108 -49 87 49 -87 87 174
117 -52 135 52 -135 135 269
126 -41 179 41 -179 179 358
135 -15 215 15 -215 215 431
144 21 241 -21 -241 230 481
153 65 252 -65 -252 252 504
162 113 249 -113 -249 249 498
171 159 231 -159 231 231 462

TABLE 2: SHEAR STRESS ON DISCRETE PLANES

ϕ (°) τ1 (MPa) τ2 (MPa) τ3 (MPa) τ4 (MPa) Δτ (MPa)

0 115 -115 -115 115 230
9 79 -141 -79 141 281
18 35 -152 -35 152 304
27 -13 -149 13 149 298
36 -59 -131 59 131 262
45 -100 -100 100 100 200
54 -131 -59 131 59 262

7 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

TABLE 2: SHEAR STRESS ON DISCRETE PLANES

ϕ (°) τ1 (MPa) τ2 (MPa) τ3 (MPa) τ4 (MPa) Δτ (MPa)

63 -149 -13 149 13 298
72 -152 35 152 -35 304
81 -141 79 141 -79 281
90 -115 115 115 -115 231
99 -79 141 79 -141 281
108 -35 152 35 -152 304
117 13 149 -13 -149 298
126 59 131 -59 -131 262
135 100 100 -100 -100 200
144 131 59 -131 -59 262
153 149 13 -149 -13 298
162 152 -35 -152 35 304
171 141 -79 -141 79 281
The Findley criterion searches for a critical plane where a combination between the
shear stress range and the normal stress is highest. This stress is shown in the second
column of Table 3, where it can be seen that 202 MPa is found at planes oriented at
18° and 162°. This results into a usage factor of 202/213 = 0.95, which is in good
agreement with the computed results.
TABLE 3: FATIGUE STRESS

ϕ (°) Δτ/2+0.20σn (MPa) Δτ/2+0.27σn (MPa) Δσn(MPa)

0 155 169 400
9 187 203 462
18 202 219 498
27 199 217 504
36 179 196 481
45 143 158 431
54 167 179 358
63 176 185 269
72 170 175 174
81 149 152 81
90 115 115 0
99 149 152 81

8 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

TABLE 3: FATIGUE STRESS

ϕ (°) Δτ/2+0.20σn (MPa) Δτ/2+0.27σn (MPa) Δσn(MPa)

108 170 175 174
117 176 185 269
126 167 179 358
135 143 158 431
144 179 196 481
153 199 217 504
162 202 219 498
171 187 203 462
The Matake criterion selects the critical plane as the one with the largest shear stress
range. This occurs at orientations 18°, 72°, 108°, and 162° where the range is
304 MPa, see Table 2. The maximum normal stress on those planes is either 87 MPa
or 249 MPa, see Table 1 and the Matake stress is either 219 MPa or 175 MPa, see
Table 3. Since the Matake criterion does not contain the normal stress in the selection
of the critical plane, the fatigue criteria is calculated with either one of them. In
Figure 5, this feature is demonstrated by the non smooth results. The analytical values
at the critical planes for the Matake usage factor are 219/223 = 0.98 and 175/
223 = 0.78, which is in good agreement with the computed results.

The Normal stress criterion considers the plane with the largest normal stress range.
The last column of Table 3 shows that this is found on planes with orientations 27°
and 153° where the range is 504 MPa. This results into a fatigue usage factor of 504/
576 = 0.88, which is in good agreement with the computed results.

Notes About the COMSOL Implementation

In the critical plane evaluation, a search resolution Q indicates the number of
evaluation point along 90° of a unit circle. The angle between two evaluation points is
then at most (90°)/(Q − 1).

Application Library path: Fatigue_Module/Verification_Examples/

cylindrical_test_specimen

9 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

Moment 22.672 [N*m] 22.67 N·m Twisting moment
Force 15708 [N] 1.571E4 N Normal force

GEOMETRY 1

Work Plane 1 (wp1)

On the Geometry toolbar, click Work Plane.

Bézier Polygon 1 (b1)

1 On the Geometry toolbar, click Primitives and choose Bézier Polygon.

10 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 1, set yw to -50.
5 In row 2, set yw to 50.
6 Find the Added segments subsection. Click Add Linear.
7 Find the Control points subsection. In row 2, set xw to 10.
8 Find the Added segments subsection. Click Add Linear.
9 Find the Control points subsection. In row 2, set yw to 30.
10 Find the Added segments subsection. Click Add Linear.
11 Find the Control points subsection. In row 2, set xw to 6.01.
12 Find the Added segments subsection. Click Add Quadratic.
13 Find the Control points subsection. In row 2, set xw to 5.25.
14 In row 2, set yw to 25.
15 In row 3, set yw to 20.
16 In row 3, set xw to 5.
17 Find the Added segments subsection. Click Add Linear.
18 Find the Control points subsection. In row 2, set yw to 0.
19 Find the Added segments subsection. Click Add Linear.
20 Find the Control points subsection. In row 2, set yw to -20.
21 Find the Added segments subsection. Click Add Quadratic.
22 Find the Control points subsection. In row 2, set xw to 5.25.
23 In row 2, set yw to -25.
24 In row 3, set yw to -30.
25 In row 3, set xw to 6.01.
26 Find the Added segments subsection. Click Add Linear.
27 Find the Control points subsection. In row 2, set xw to 10.
28 Find the Added segments subsection. Click Add Linear.
29 Find the Control points subsection. In row 2, set yw to -50.
30 Find the Added segments subsection. Click Add Linear.
31 Find the Control points subsection. In row 2, set xw to 0.
32 Right-click Bézier Polygon 1 (b1) and choose Build Selected.

11 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

Fillet 1 (fil1)
1 On the Geometry toolbar, click Fillet.
2 On the object b1, select Points 6 and 7 only.
3 In the Settings window for Fillet, locate the Radius section.
4 In the Radius text field, type 5.

Revolve 1 (rev1)
1 On the Geometry toolbar, click Revolve.
2 In the Model Builder window, right-click Revolve 1 (rev1) and choose Build All Objects.
3 Click the Zoom Extents button on the Graphics toolbar.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 210e9 Pa Basic
Poisson's ratio nu 0.3 1 Basic
Density rho 7800 kg/m³ Basic

SOLID MECHANICS (SOLID)

Fixed Constraint 1
1 On the Physics toolbar, click Boundaries and choose Fixed Constraint.
2 Select Boundaries 3, 4, 21, and 23 only.

Rigid Connector 1
1 On the Physics toolbar, click Boundaries and choose Rigid Connector.
2 Select Boundaries 11, 12, 40, and 41 only.

Applied Force 1
1 Right-click Rigid Connector 1 and choose Applied Force.
2 In the Settings window for Applied Force, locate the Applied Force section.

12 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

3 Specify the F vector as

0 x
Force y
0 z

4 On the Physics toolbar, click Load Group and choose New Load Group.

Applied Moment 1
1 In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid)
right-click Rigid Connector 1 and choose Applied Moment.
2 In the Settings window for Applied Moment, locate the Applied Moment section.
3 Specify the M vector as

0 x
Moment y
0 z

4 On the Physics toolbar, click Load Group and choose New Load Group.

GLOBAL DEFINITIONS
1 In the Model Builder window, expand the Global Definitions node, then click Load
Group 1 (lg1).
2 In the Settings window for Load Group, type lgf in the Parameter name text field.
3 In the Model Builder window, under Global Definitions click Load Group 2 (lg2).
4 In the Settings window for Load Group, type lgm in the Parameter name text field.
Create a load cycle consisting of four load cases.

STUDY 1

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Define load cases check box.
4 Click Add four times.

13 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

5 In the table, enter the following settings:

6 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
Verify stresses at the center of the test specimen.

1D Plot Group 2
On the Home toolbar, click Add Plot Group and choose 1D Plot Group.

Point Graph 1
On the 1D Plot Group 2 toolbar, click Point Graph.

1D Plot Group 2
1 Select Point 20 only.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid
Mechanics>Stress (Gauss points)>Second Piola-Kirchhoff stress, Gauss-point evaluation
(Material)>[Link] - Second Piola-Kirchhoff stress, Gauss-point evaluation, Y
component.
3 Click to expand the Legends section. Select the Show legends check box.
4 From the Legends list, choose Manual.
5 In the table, enter the following settings:

Legends
normal stress

6 Right-click Point Graph 1 and choose Duplicate.

7 In the Settings window for Point Graph, locate the y-Axis Data section.
8 In the Expression text field, type [Link].

14 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

9 Locate the Legends section. In the table, enter the following settings:

Legends
shear stress

10 Right-click Results>1D Plot Group 2>Point Graph 2 and choose Duplicate.

11 In the Settings window for Point Graph, locate the y-Axis Data section.
12 In the Expression text field, type [Link].
13 Locate the Legends section. In the table, enter the following settings:

Legends
von Mises

14 On the 1D Plot Group 2 toolbar, click Plot.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

5 Click Add to Component in the window toolbar.

FATIGUE (FTG)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue (ftg).
2 In the Settings window for Fatigue, type Fatigue Findley in the Label text field.
3 Right-click Component 1 (comp1)>Fatigue Findley and choose the boundary
evaluation Stress-Based.

FATIGUE FINDLEY (FTG)

On the Physics toolbar, click Fatigue (ftg) and choose Fatigue Findley (ftg).

Stress-Based 1
1 In the Settings window for Stress-Based, locate the Boundary Selection section.
2 From the Selection list, choose All boundaries.

15 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

3 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

4 Click Add to Component in the window toolbar.

FATIGUE 2 (FTG2)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue 2 (ftg2).
2 In the Settings window for Fatigue, type Fatigue Matake in the Label text field.
3 Right-click Component 1 (comp1)>Fatigue Matake and choose the boundary
evaluation Stress-Based.

FATIGUE MATAKE (FTG2)

On the Physics toolbar, click Fatigue 2 (ftg2) and choose Fatigue Matake (ftg2).

Stress-Based 1
1 In the Settings window for Stress-Based, locate the Boundary Selection section.
2 From the Selection list, choose All boundaries.
3 Locate the Fatigue Model Selection section. From the Criterion list, choose Matake.
4 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

16 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

4 Click Add to Component in the window toolbar.

5 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 3 (FTG3)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue 3 (ftg3).
2 In the Settings window for Fatigue, type Fatigue Normal Stress in the Label text
field.
3 Right-click Component 1 (comp1)>Fatigue Normal Stress and choose the boundary
evaluation Stress-Based.

FATIGUE NORMAL STRESS (FTG3)

On the Physics toolbar, click Fatigue 3 (ftg3) and choose Fatigue Normal Stress (ftg3).

Stress-Based 1
1 In the Settings window for Stress-Based, locate the Boundary Selection section.
2 From the Selection list, choose All boundaries.
3 Locate the Fatigue Model Selection section. From the Criterion list, choose Normal
stress.
4 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

MATERIALS

Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose All boundaries.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Normal stress sensitivity k_Findley 0.20 1 Findley
coefficient
Limit factor f_Findley 213 [MPa] Pa Findley
Normal stress sensitivity k_Matake 0.27 1 Matake
coefficient

17 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

Property Name Value Unit Property group

Limit factor f_Matake 223 [MPa] Pa Matake
Limit factor f_NormalS 576 [MPa] Pa Normal stress
tress

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, expand the Study 1 node, then click Study 2>Step 1:
Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 1, Stationary.
6 On the Home toolbar, click Compute.

RESULTS

Fatigue Usage Factor (ftg)

1 In the Settings window for 3D Plot Group, type Fatigue Usage Factor
(Findley) in the Label text field.

18 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

Fatigue Usage Factor (ftg2)

1 In the Model Builder window, under Results click Fatigue Usage Factor (ftg2).
2 In the Settings window for 3D Plot Group, type Fatigue Usage Factor (Matake)
in the Label text field.

Fatigue Usage Factor (ftg3)

1 In the Model Builder window, under Results click Fatigue Usage Factor (ftg3).
2 In the Settings window for 3D Plot Group, type Fatigue Usage Factor (Normal
Stress) in the Label text field.

19 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

20 | H I G H - C Y C L E F A T I G U E A N A L Y S I S O F A C Y L I N D R I C A L TE S T S P E C I M E N
 Solved with COMSOL Multiphysics 5.2

High-Cycle Fatigue of a
Reciprocating Piston Engine
Introduction
In a reciprocating piston engine the connecting rods transfer rotating motion into
reciprocating motion. The connecting rods are constantly under high stresses and the
load increases with the engine speed. A failure of one part in the engine usually results
in a replacement of the whole engine. It is therefore of crucial importance to design all
engine parts so that none of them fail during the operational lifetime of the engine.
The connecting rods are identified as the critical parts and are analyzed from the
fatigue perspective. The fatigue lifetime is predicted using the Basquin high-cycle
fatigue criteria.

This example is based on an application from the Multibody Dynamics Module,

Three-Cylinder Reciprocating Engine, where the critical part of the engine is modeled
as a flexible body while the remaining parts are modeled rigid. The connections
between different parts is obtained by using different type of joints. This technique
significantly reduces the model size while maintaining the force equilibrium in the
assembly.

Model Definition
The three-cylinder engine is presented in Figure 1 and it operates at 1000 RPM. Its
material data is taken from structural steel. Additional information regarding its set up
can be found in the documentation of the application Three-Cylinder Reciprocating
Engine, found in the Multibody Dynamics Module.

The material data from fatigue tests is summarized in Figure 2. The Basquin relation
with the material constants σf’ = 1043 MPa and b = -0.116 gives a good fit to the
experimental results. The material exhibits a fatigue limit at 210 MPa.

1 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

Figure 1: Geometry of the reciprocating piston engine.

Figure 2: Fatigue material curve.

2 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

Results and Discussion

At first, stresses in a fillet of the piston end of the connecting rod are examined, see
Figure 3. A fillet is chosen since a stress concentration due to geometrical change is
expected there. The engine needs few revolutions before it obtains a steady state
behavior. From cycle three, the stress history for each consecutive cycles seems to
repeat itself. The peak stress is about the same just as the rest of the stress cycle. The
stress history is dominated by the third principal stress since the connecting rod is in
compression. The two other principal stresses are small so that the stress state at the
fillet can be considered uniaxial. Therefore the principal stress is taken as the amplitude
stress in the Basquin relation as opposed to the von Mises stress which would be more
appropriate in a multiaxial loading. The fatigue life prediction is shown in Figure 4.

Figure 3: Stress history in the connecting rod.

3 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

Figure 4: Fatigue life prediction in the connecting rod.

The critical point is at the fillet close to the top end of the connecting rod where the
Basquin model predicts a fatigue life that is longer than twenty-five billion cycles. This
is an extremely long life. It can therefore be expected that the stress in the assembly is
below the endurance limit that for the used material is 210 MPa. By using the Basquin
relation

9 – 0.116
σ a = 1.043 ⋅ 10 ( N )

where σa is the stress amplitude and N is the number of cycles to failure, the fatigue
endurance limit life can be back-calculated to 245 million cycles. This is less than the
calculated value, see Figure 4, and therefore the connecting rod is designed for infinite
life. This could have been already observed in Figure 3 where the stress history is
shown. Since the principal stress range is about 110 MPa the stress amplitude is about
55 MPa and that is below the endurance limit of the material.

4 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

Notes About the COMSOL Implementation

In the Basquin evaluation the parameter Cycle cutoff can be used to incorporate the
effect of endurance limit. The number of cycles that gives the endurance limit must
then be back-calculated.

Application Library path: Fatigue_Module/Stress_Life/engine_fatigue

Modeling Instructions
In this example you will start from an existing model from the Multibody Dynamics
Module.

From the File menu, choose Open.

Under the Application Library root, browse to the folder

Multibody_Dynamics_Module/Automotive_and_Aerospace and double-click the
file reciprocating_engine.mph.

ROOT
If the model was stored without solutions, you will now have to run Study 1 and Study
2 before continuing.

Evaluate how stresses develop. Examine one point in a fillet close to the small end of
the connecting rod.

1D Plot Group 11
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Stress history: Connecting
rod in the Label text field.

3 Locate the Data section. From the Data set list, choose Study: Multibody analysis/
Solution 2 (3) (sol2).

Point Graph 1
On the Stress history: Connecting rod toolbar, click Point Graph.

Stress history: Connecting rod

1 Select Point 834 only.
2 In the Settings window for Point Graph, locate the y-Axis Data section.

5 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

3 In the Expression text field, type mbd.sp1.

4 Locate the x-Axis Data section. From the Parameter list, choose Expression.
5 In the Expression text field, type theta/(2*pi).
6 Click to expand the Legends section. Select the Show legends check box.
7 From the Legends list, choose Manual.
8 In the table, enter the following settings:

Legends
First principal stress

9 On the Stress history: Connecting rod toolbar, click Plot.

10 Right-click Point Graph 1 and choose Duplicate.
11 In the Settings window for Point Graph, locate the y-Axis Data section.
12 In the Expression text field, type mbd.sp2.
13 Locate the Legends section. Select the Show legends check box.
14 From the Legends list, choose Manual.
15 In the table, enter the following settings:

Legends
Second principal stress

16 On the Stress history: Connecting rod toolbar, click Plot.

17 Right-click Results>Stress history: Connecting rod>Point Graph 2 and choose
Duplicate.
18 In the Settings window for Point Graph, locate the y-Axis Data section.
19 In the Expression text field, type mbd.sp3.
20 Locate the Legends section. Select the Show legends check box.
21 From the Legends list, choose Manual.
22 In the table, enter the following settings:

Legends
Third principal stress

23 On the Stress history: Connecting rod toolbar, click Plot.

24 Right-click Results>Stress history: Connecting rod>Point Graph 3 and choose
Duplicate.

6 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

25 In the Settings window for Point Graph, locate the y-Axis Data section.
26 In the Expression text field, type [Link].
27 Locate the Legends section. Select the Show legends check box.
28 From the Legends list, choose Manual.
29 In the table, enter the following settings:

Legends
Effective von Mises stress

30 On the Stress history: Connecting rod toolbar, click Plot.

31 In the Model Builder window, click Stress history: Connecting rod.
32 In the Settings window for 1D Plot Group, locate the Plot Settings section.
33 Select the x-axis label check box.
34 In the associated text field, type Rotation of crankshaft (cycle).
35 Select the y-axis label check box.
36 In the associated text field, type Stress (Pa).
37 Click to expand the Title section. From the Title type list, choose Manual.
38 In the Title text area, type Stress history in a connecting rod fillet.
39 Click to expand the Legend section. From the Position list, choose Lower left.
40 On the Stress history: Connecting rod toolbar, click Plot.
The results indicate a mechanical response that repeats itself from the third cycle for
each consecutive load cycle. Recalculate the 6:th load cycle using a better time
discretization around the peak stress.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies.
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Heat Transfer in Fluids (ht)
Coefficient Form PDE (c)

7 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

5 Find the Studies subsection. In the Select study tree, select Preset Studies>Time
Dependent.
6 Click Add Study in the window toolbar.
7 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 3

Step 1: Time Dependent

1 In the Model Builder window, under Study 3 click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, type range(0.1368,4e-4,0.14)
range(0.1401,1e-4,0.1440) range(0.1444,4e-4,0.16).

4 Click to expand the Values of dependent variables section. Locate the Values of
Dependent Variables section. Find the Initial values of variables solved for subsection.
From the Settings list, choose User controlled.
5 From the Method list, choose Solution.
6 From the Study list, choose Study: Multibody analysis, Time Dependent.
7 From the Time (s) list, choose 0.1368.
8 In the Model Builder window, click Study 3.
9 In the Settings window for Study, type Study: Fatigue step in the Label text field.
10 On the Home toolbar, click Compute.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study: Thermodynamic analysis
Study: Multibody analysis
Study: Fatigue step

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

8 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

FATIGUE (FTG)

Stress-Life 1
1 In the Model Builder window, under Multibody Analysis (comp2) right-click Fatigue
(ftg) and choose the boundary evaluation Stress-Life.
2 Select boundaries of the connecting rod in the middle.
3 In the Settings window for Stress-Life, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Basquin.
5 Locate the Solution Field section. From the Physics interface list, choose Multibody
Dynamics (mbd).
6 Locate the Fatigue Model Parameters section. From the σf’ list, choose User defined.
In the associated text field, type 1.043e9.
7 From the b list, choose User defined. In the associated text field, type -0.116.
The cutoff value can be used to specify the endurance limit. In this example the
cutoff is set to a high value in order to examine how the Basquin model predicts
lifetime in case the material did not have an endurance limit.
8 Locate the Evaluation Settings section. In the Ncut text field, type 1e20.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies.
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Heat Transfer in Fluids (ht)
Coefficient Form PDE (c)
Multibody Dynamics (mbd)

5 Find the Studies subsection. In the Select study tree, select Preset Studies>Stationary.
6 Click Add Study in the window toolbar.

STUDY 4

Step 1: Stationary
1 In the Model Builder window, under Study 4 click Step 1: Stationary.

9 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

2 In the Settings window for Stationary, click to expand the Values of dependent
variables section.
3 Locate the Values of Dependent Variables section. Find the Values of variables not
solved for subsection. From the Settings list, choose User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study: Fatigue step, Time Dependent.
6 From the Time (s) list, choose All.
7 In the Model Builder window, click Study 4.
8 In the Settings window for Study, type Study: Fatigue analysis in the Label text
field.
9 On the Home toolbar, click Compute.

10 | HIGH-CYCLE FATIGUE OF A RECIPROCATING PISTON ENGINE

 Solved with COMSOL Multiphysics 5.2

Fa ti g ue Fa i lur e of an E yegl ass Fram e

Introduction
In the search for weight reduction, the cross-section of an eyeglass frame is
continuously reduced. The thin section over the nose transfers the entire load between
the two halves. This example predicts the fatigue life using the combined Basquin and
Coffin-Manson model when eyeglasses are subjected to bending.

Model Definition
The frame of the eyeglasses is made of the MONEL alloy 400. The eyeglass lenses are
made of a CR-39 material that is lighter in weight than glass. The risk of fatigue in the
lenses is avoided by using a coating that holds together any shards in case of fracture.
Therefore the fatigue life of the frame is only predicted.

The frame of the eyeglasses is very thin, 1 mm, and has a shape according to Figure 1.

Figure 1: Shape of the eyeglasses.

Young’s modulus of MONEL alloy 400 is taken from COMSOL Material Library and
the Poisson’s ratio 0.32. For the CR-39 material, the Young’s modulus is 2.1 GPa and
the Poisson’s ratio is 0.4.

Fatigue data for the MONEL alloy 400 has been obtained in rotating bending tests
and fitted to the combined Basquin and Coffin-Manson relation according to where
the strain amplitude, εa, is expressed as

′
σf –b ′ –c
ε a = ----- ⋅ N + ε f ⋅ N
E

where εa is the strain amplitude, E is the Young’s modulus, N represents the number
of cycles to failure, and σf’, εf’, b, and c are fatigue material parameters.

Bending of the eyeglasses is simulated by fixing one side of the eyeglasses and applying
an alternating vertical 4 N force on the other side of the eyeglasses.

1 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

Results and Discussion

The stress distribution in the eyeglasses at the peak load is shown in Figure 2.

Figure 2: Effective stress in eyeglasses.

Bending of eyeglasses causes high stresses on both sides of the thin central part. Since
the effective stress does not discriminate between tension and compression the same
stress levels are encountered on both sides. The resulting stress contours are similar if
the bending is reversed when the right part of the glasses is pulled down instead.

Since the fatigue data is obtained in a rotating bending test, the stresses and strains
alternate between tension and compression during one test cycle. This is exactly the
same structural behavior as in one load cycle of the eyeglasses and therefore the fatigue
curve is directly applicable to the resulting principal strains. The highest and smallest
principal strains at peak loads are shown in Figure 3.

2 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

Figure 3: Principal strains in the eyeglasses. Above first principal strain when pulling up
and below third principal strain when pulling down.

When pulling up, the lower side of the thin central section experiences peak tensile
strains. The peak compressive strains are experienced in the same point when pulling
the eyeglasses down. On the upper side of the thin central section over the nose, an
opposite situation is encountered, with peak compression when pulling up and peak
tension when pulling down. The difference in strain during both load events controls
the fatigue life that is predicted in Figure 4.

3 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

Figure 4: Fatigue life according to the combined Basquin and Coffin-Manson relation.

Application Library path: Fatigue_Module/Strain_Life/

eyeglass_frame_fatigue

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.

4 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

5 In the Select study tree, select Preset Studies>Stationary.

6 Click Done.

GEOMETRY 1
1 In the Model Builder window, click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Import 1 (imp1)
1 On the Home toolbar, click Import.
2 In the Settings window for Import, locate the Import section.
3 Click Browse.
4 Browse to the application’s Application Library folder and double-click the file
eyeglass_frame_fatigue.mphbin.

5 Click Import.
The thin central part of the eyeglasses transfers the entire load between the two
halves. Make a rectangle in the center to create a domain for fine structured mesh.

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size and Shape section.
3 In the Width text field, type 10.
4 In the Height text field, type 4.
5 Locate the Position section. In the x text field, type 55.
6 In the y text field, type -8.

Intersection 1 (int1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Intersection.
2 Select the objects imp1 and r1 only.
3 In the Settings window for Intersection, locate the Intersection section.
4 Select the Keep input objects check box.

Delete Entities 1 (del1)

1 Right-click Geometry 1 and choose Delete Entities.
2 In the Settings window for Delete Entities, locate the Entities or Objects to Delete
section.

5 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

3 From the Geometric entity level list, choose Domain.

4 On the object r1, select Domain 1 only.
5 Click the Build All Objects button.

MATERIALS
On the Home toolbar, click Windows and choose Add Material.

ADD MATERIAL
1 Go to the Add Material window.
2 In the Search text field, type monel 400.
3 Click Search.
4 In the tree, select Material Library>Nickel Alloys>Monel 400 (UNS N04400) (NW
4400)>Monel 400 (UNS N04400) (NW 4400) [solid]>Monel 400 (UNS N04400) (NW 4400)
[solid,annealed].
5 Select Domains 1–3 only.
6 Click Add to Selection in the window toolbar.

MATERIALS

Monel 400 (UNS N04400) (NW 4400) [solid,annealed] (mat1)

1 In the Model Builder window, under Component 1 (comp1)>Materials click Monel 400
(UNS N04400) (NW 4400) [solid,annealed] (mat1).
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Poisson's ratio nu 0.32 1 Young's modulus and Poisson's ratio

Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, type CR-39 in the Label text field.
3 Select Domains 4 and 5 only.

6 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 2.1[GPa] Pa Basic
Poisson's ratio nu 0.4 1 Basic
Density rho 1.3[g/cm^3] kg/m³ Basic

SOLID MECHANICS (SOLID)

1 In the Model Builder window, under Component 1 (comp1) click Solid Mechanics
(solid).
2 In the Settings window for Solid Mechanics, locate the Thickness section.
3 In the d text field, type 0.001.

Fixed Constraint 1
1 On the Physics toolbar, click Boundaries and choose Fixed Constraint.
2 Select Boundary 3 only.

Boundary Load 1
1 On the Physics toolbar, click Boundaries and choose Boundary Load.
2 Select Boundary 42 only.
3 In the Settings window for Boundary Load, locate the Force section.
4 From the Load type list, choose Total force.
5 Specify the Ftot vector as

0 x
1 y

6 On the Physics toolbar, click Load Group and choose New Load Group.

MESH 1

Mapped 1
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and
choose Mapped.
2 In the Settings window for Mapped, locate the Domain Selection section.
3 From the Geometric entity level list, choose Domain.
4 Select Domain 2 only.

7 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

Distribution 1
1 Right-click Component 1 (comp1)>Mesh 1>Mapped 1 and choose Distribution.
2 Select Boundary 1 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 15.

Distribution 2
1 Right-click Mapped 1 and choose Distribution.
2 Select Boundary 22 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 20.

Free Triangular 1
1 In the Model Builder window, right-click Mesh 1 and choose Free Triangular.
2 Right-click Free Triangular 1 and choose Build All.

STUDY 1

Step 1: Stationary
1 In the Model Builder window, under Study 1 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Define load cases check box.
4 Click Add four times.
5 In the table, enter the following settings:

Load case lg1 Weight

Load case 1 √ 0
Load case 2 √ 4
Load case 3 √ -4
Load case 4 √ 0

6 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
1 In the Settings window for 2D Plot Group, locate the Data section.
2 From the Load case list, choose Load case 2.

8 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

3 In the Model Builder window, expand the Stress (solid) node.

4 In the Model Builder window, expand the Results>Stress (solid)>Surface 1 node, then
click Deformation.
5 In the Settings window for Deformation, locate the Scale section.
6 Select the Scale factor check box.
7 In the associated text field, type 1.
8 On the Stress (solid) toolbar, click Plot.
The effective stress does not discriminate between tension and compression.
Evaluate how strains change at peak loads.

Stress (solid) 1
1 In the Model Builder window, right-click Stress (solid) and choose Duplicate.
2 In the Settings window for 2D Plot Group, type Principal strain (solid) in
the Label text field.
3 Locate the Plot Settings section. Clear the Plot data set edges check box.
4 Click to expand the Color legend section. Locate the Color Legend section. From the
Position list, choose Bottom.

Principal strain (solid)

1 In the Model Builder window, expand the Results>Principal strain (solid) node, then
click Surface 1.
2 In the Settings window for Surface, locate the Data section.
3 From the Data set list, choose Study 1/Solution 1 (sol1).
4 From the Load case list, choose Load case 2.
5 Locate the Expression section. In the Expression text field, type solid.ep1.
Display both the largest principal strain and the smallest principal strain in the same
figure.
6 In the Model Builder window, expand the Results>Principal strain (solid)>Surface 1
node, then click Deformation.
7 In the Settings window for Deformation, locate the Expression section.
8 In the Y component text field, type v+0.020.
9 In the Model Builder window, under Results>Principal strain (solid) right-click Surface
1 and choose Duplicate.
10 In the Settings window for Surface, locate the Data section.

9 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

11 From the Load case list, choose Load case 3.

12 Locate the Expression section. In the Expression text field, type solid.ep3.
13 In the Model Builder window, expand the Results>Principal strain (solid)>Surface 2
node, then click Deformation.
14 In the Settings window for Deformation, locate the Expression section.
15 In the Y component text field, type v-0.020.
16 On the Principal strain (solid) toolbar, click Plot.

COMPONENT 1 (COMP1)
On the Home toolbar, click Windows and choose Add Physics.

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1

4 Click Add to Component in the window toolbar.

FATIGUE (FTG)

Strain-Life 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue (ftg)
and choose the domain evaluation Strain-Life.
2 Select Domains 1–3 only.
3 In the Settings window for Strain-Life, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Combined Basquin and Coffin-Manson.
5 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
6 Locate the Evaluation Settings section. In the Ncut text field, type 5.75e6.

10 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

MATERIALS

Monel 400 (UNS N04400) (NW 4400) [solid,annealed] (mat1)

1 In the Model Builder window, under Component 1 (comp1)>Materials click Monel 400
(UNS N04400) (NW 4400) [solid,annealed] (mat1).
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Fatigue strength coefficient sigmaf_Basquin 970 [MPa] Pa Basquin
Fatigue strength exponent b_Basquin -0.077 1 Basquin
Fatigue ductility coefficient epsilonf_CM 0.738 1 Coffin-Manson
Fatigue ductility exponent c_CM -0.54 1 Coffin-Manson

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 In the Select study tree, select Preset Studies.
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

5 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
6 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.

11 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

5 From the Study list, choose Study 1, Stationary.

6 On the Home toolbar, click Compute.

12 | FATIGUE FAILURE OF AN EYEGLASS FRAME

 Solved with COMSOL Multiphysics 5.2

Fatigue Response of a Random

Non-Proportional Load
Introduction
A thin-walled frame member with a central cutout is subjected to a random load
scenario. Although the stresses are expected to be far below the yield level of the
material, a concern arises whether or not the component fails due to fatigue.

This example demonstrates an approach to damage quantification of a long load

history. The Rainflow counting algorithm is used to define the load scenario and
Palmgren-Miner linear damage model quantifies the damage.

Model Definition
The load carrying beam is shown in Figure 1. It has a length of 1.1 m, a thin-walled
square cross section having the dimension 160 mm x160 mm and a thickness of 6 mm.
The cutout is centrally placed on one of the faces and is 100 mm long, 80 mm wide
and has a fillet with radius 10 mm in each corner.

Figure 1: Geometry of the frame member.

The loading consists of bending moments in both directions and a twisting moment.
All three loads can vary independently. The information about the load is obtained
using three strain gauges which are glued on the bottom and the back sides of the
frame. The back side is the one opposite of the face with the hole. The location and
position of the strain gauges are shown in Figure 2. Strain gauges 1 and 2 are oriented
at a 45° angle to the global coordinate system while the 3rd gauge is aligned with the

1 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

length of the face (the x-axis). Both are located 850 mm from the edge along the
length and in the middle along the width of the faces.

Figure 2: Location and orientation of the strain gauges.

The frame is made out of steel with Young’s modulus, Poisson's ratio and density given
by E = 200 GPa, ν = 0.33 and ρ = 7800 kg/m3.

The fatigue behavior follows material library data for iron alloy 4340 for variant
defined by phase UTS 200 Ksi - 293K and variation unnotched.

The left end of the structure, x = 0, is clamped. One twisting moment, along x, and
two bending moments, along y and z, are applied on the right end at x = 1.1 m. The
lifetime for which the frame is designed is 10 000 longer than the what has been
recorded be the strain gauges. The response to one loading cycle block in all three
gauges is captured in Figure 3. From the history it is clear that the loading event is
non-proportional. However stress concentrations arise around the cutout where the
stress state is uniaxial. Therefore used models, Rainflow Counting and
Palmgren-Miner, can be seen as appropriate for fatigue evaluation.

2 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Figure 3: Strain history in tree gauges. The subscript indicates the gauge number.

In order to apply the strain history to the frame, the strain history needs to be
transferred to loads which can be applied on the boundary. The transformation can be
done in two ways. In the first one a relation between a unit moment and a response in
each strain gauge must be obtained using COMSOL. By repeating it for each unit
moment and inverting the relation a transformation matrix relating strains to moments
is obtained.

An alternative way is to use an approximative analytical relation based on beam theory

with a thin-walled assumption, Hooke’s law, and rotation of stresses in a plane. The
result is given below without a detailed derivation.

3 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

ε 1 = B 1 ( A 1 C 1 + A 2 D 1 )M Y + 2B 2 ( A 1 C 3 + A 2 D 3 )M X
ε 2 = B 1 ( A 1 D 1 + A 2 C 1 )M Y + 2B 2 ( A 1 D 3 + A 2 C 3 )M X (1)
ε 3 = A1 B1 MZ

3 2
where A 1 = 1 ⁄ E , A 2 = – ν ⁄ E , B 1 = 3 ( b + t ) ⁄ ( 4tb ) , B 2 = 1 ⁄ ( 2tb ) ,
C 1 = cos 45° cos 45° , C 2 = sin 45° sin 45° , C 3 = cos 45° sin 45° ,
D 1 = cos 135° cos 135° , D 2 = sin 135° sin 135° , and D 3 = cos 135° sin 135° .
The variables t = 6 mm and b = 154 mm define the thickness and the side (using
midsurfaces) of the cross sect ion.

Based on the geometrical and material constants the Equation 1 gives the following
relation

ε1 2.34e-8 -9.18e-9 0 MX
ε2 = -2.34e-8 -9.18e-9 0 MY (2)
ε3 0 0 -2.73e-8 M Z

Results and Discussion

According to FE analysis the applied moments and gauge strains are related by

ε1 2.62e-8 -9.45e-9 0 MX
ε2 = -2.62e-8 -9.55e-9 0 MY (3)
ε3 0 0 -2.86e-8 M Z

The coefficients differ somewhat when compared to the analytical results, Equation 2.
It is reasonable to assume that the final fatigue prediction also differs, depending on
used transformation matrix. Since Equation 2 is based on approximations further
results are based on the FE-relation, Equation 3.

4 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Stresses on the inner side of the shell are more critical from the fatigue point of view.
The variation of the fatigue usage factor along the fillets is shown in Figure 4. It
reaches about 0.11 and thus the examined frame should not fail in fatigue.

Figure 4: Fatigue usage factor along the cutout fillets. The angle is measured in the
xz-plane with the angle starting from the x-axis.

In the most loaded point in the structure the stress history seems to be fairly symmetric
around a zero mid stress. The mid stress is found in the range from -250 MPa to
250 MPa and the amplitude extend almost up to 600 MPa. The load distribution in
the critical point is captured and shown in the Rainflow histogram in Figure 5.

5 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Figure 5: Cycle counting, following the rainflow theory, in the point of the highest fatigue
usage factor.

In Figure 6 the relative contribution to the damage is shown for the same location.
The dark blue area indicates stress cycles which are non-damaging. This means that
they are below the endurance limit in the Wöhler curve, also called the S-N curve. The
damaging stress cycles are found for high stress amplitudes. Since the Palmgren-Miner
rule scales damage linearly with the number of cycles, when the number of cycles
increases with factor 1/0.11 the fatigue usage factor exceeds 1 and thus a failure occurs.
In practice, the linearity assumption of the Palmgren-Miner rule can be questioned, so
a proper safety factor should be applied.

An important information when evaluating Figure 5 and Figure 6 is that 37% of the
fatigue damage comes from one single event in the load history and that most damage
is caused by only few load cycles. This indicates that the load history recorded is to
short to make good predictions. Either new, longer, measurements should be made,
or a high safety factor should be used in combination with some statistical
considerations.

6 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Figure 6: Relative fatigue usage, following the Palmgren-Miner damage rule, in the
point of the highest usage factor.

Notes About the COMSOL Implementation

In COMSOL, several functions with the same argument list can be put in one file. This
is demonstrated in the example with three load functions define the moment
prescribed on one end of the frame. When an interpolation function is read from a file
it treats the first columns as arguments and the following one as function response. As
an example, in the interpolation function defining the moment around the y-axis the
number of arguments is assigned 1 and the function position is assigned 2. Thus the
first column is the argument and the third column is the function value (1+2).

COMSOL provides several options for specification of the S-N curve. Interpolation
function with options data type Grid, interpolation Linear and extrapolation Constant is
recommended. Those options are optimal for search of the life once R-value, and the
amplitude stress are known. The input for S-N curve defined in the example is shown
in Figure 7.

7 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Figure 7: Wöhler curve of the material.

The fatigue material data in the material library gives the maximum fatigue stress,
σ max , as the function of the number of cycles and R-value. In order to transfer it to
the stress amplitude, σ a ,which is required by the Cumulative Damage fatigue feature,
just apply following transformation

σ max ( 1 – R )
σ a = --------------------------------
2

Application Library path: Fatigue_Module/Damage/frame_with_cutout

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

8 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select physics tree, select Structural Mechanics>Shell (shell).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1

Block 1 (blk1)
1 On the Geometry toolbar, click Block.
2 In the Settings window for Block, locate the Size and Shape section.
3 In the Width text field, type 1.1.
4 In the Depth text field, type 0.154.
5 In the Height text field, type 0.154.
6 Locate the Position section. From the Base list, choose Center.

Work Plane 1 (wp1)

1 On the Geometry toolbar, click Work Plane.
2 In the Settings window for Work Plane, locate the Plane Definition section.
3 From the Plane type list, choose Face parallel.
4 Find the Planar face subsection. Select the Active toggle button.
5 On the object blk1, select Boundary 3 only.

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Position section.
3 From the Base list, choose Center.
4 Locate the Size and Shape section. In the Width text field, type 0.08.
5 In the Height text field, type 0.1.

Fillet 1 (fil1)
1 On the Geometry toolbar, click Fillet.
2 On the object r1, select Points 1–4 only.
3 In the Settings window for Fillet, locate the Radius section.

9 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

4 In the Radius text field, type 0.01.

Create points for strain evaluation where strain gauges are placed.

Point 1 (pt1)
1 On the Geometry toolbar, click More Primitives and choose Point.
2 In the Settings window for Point, locate the Point section.
3 In the x text field, type 0.3.
4 In the z text field, type -0.077.

Point 2 (pt2)
1 On the Geometry toolbar, click More Primitives and choose Point.
2 In the Settings window for Point, locate the Point section.
3 In the x text field, type 0.3.
4 In the y text field, type 0.077.
5 Click the Build All Objects button.

DEFINITIONS
Create a new coordinate system that is aligned with the strain gauge directions on the
bottom side of the frame.

Base Vector System 2 (sys2)

1 On the Definitions toolbar, click Coordinate Systems and choose Base Vector System.
2 In the Settings window for Base Vector System, locate the Settings section.
3 Find the Base vectors subsection. In the table, enter the following settings:

x y z
x1 cos(pi/4) sin(pi/4) 0
x2 -sin(pi/4) cos(pi/4) 0

4 Find the Simplifications subsection. Select the Assume orthonormal check box.
5 Click the Wireframe Rendering button on the Graphics toolbar.

SHELL (SHELL)
1 In the Model Builder window, under Component 1 (comp1) click Shell (shell).
2 Select Boundaries 2–5 only.
3 In the Settings window for Shell, locate the Thickness section.
4 In the d text field, type 0.006[m].

10 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Linear Elastic Material 2

1 On the Physics toolbar, click Boundaries and choose Linear Elastic Material.
2 Select Boundary 3 only.
3 In the Settings window for Linear Elastic Material, locate the Coordinate System
Selection section.
4 From the Coordinate system list, choose Base Vector System 2 (sys2).

Prescribed Displacement/Rotation 11
1 On the Physics toolbar, click Edges and choose Prescribed Displacement/Rotation.
2 Select Edges 1, 2, 4, and 6 only.
3 In the Settings window for Prescribed Displacement/Rotation, locate the Prescribed
Displacement section.
4 Select the Prescribed in x direction check box.
5 Select the Prescribed in y direction check box.
6 Select the Prescribed in z direction check box.
Apply one twisting and two bending unit moments and differentiate them using
load cases.

Rigid Connector 1
1 On the Physics toolbar, click Edges and choose Rigid Connector.
2 Select Edges 17–20 only.

Applied Moment 1
1 Right-click Rigid Connector 1 and choose Applied Moment.
2 In the Settings window for Applied Moment, locate the Applied Moment section.
3 Specify the M vector as

M x
0 y
0 z

4 In the Label text field, type Twisting Moment (x).

Applied Moment 2
1 In the Model Builder window, under Component 1 (comp1)>Shell (shell) right-click
Rigid Connector 1 and choose Applied Moment.
2 In the Settings window for Applied Moment, locate the Applied Moment section.

11 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

3 Specify the M vector as

0 x
M y
0 z

4 In the Label text field, type Bending Moment (y).

Applied Moment 3
1 Right-click Rigid Connector 1 and choose Applied Moment.
2 In the Settings window for Applied Moment, locate the Applied Moment section.
3 Specify the M vector as

0 x
0 y
M z

4 In the Label text field, type Bending Moment (z).

Twisting Moment (x)

1 In the Model Builder window, under Component 1 (comp1)>Shell (shell)>Rigid
Connector 1 click Twisting Moment (x).
2 On the Physics toolbar, click Load Group and choose New Load Group.

Bending Moment (y)

1 In the Model Builder window, under Component 1 (comp1)>Shell (shell)>Rigid
Connector 1 click Bending Moment (y).
2 Click Load Group and choose New Load Group.

Bending Moment (z)

1 In the Model Builder window, under Component 1 (comp1)>Shell (shell)>Rigid
Connector 1 click Bending Moment (z).
2 Click Load Group and choose New Load Group.

GLOBAL DEFINITIONS
1 In the Model Builder window, expand the Global Definitions node, then click Load
Group 1 (lg1).
2 In the Settings window for Load Group, type lgX in the Parameter name text field.
3 In the Label text field, type Load Group: Mx.

12 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

4 In the Model Builder window, under Global Definitions click Load Group 2 (lg2).
5 In the Settings window for Load Group, type lgY in the Parameter name text field.
6 In the Label text field, type Load Group: My.
7 In the Model Builder window, under Global Definitions click Load Group 3 (lg3).
8 In the Settings window for Load Group, type lgZ in the Parameter name text field.
9 In the Label text field, type Load Group: Mz.

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 In the table, enter the following settings:

Name Expression Value Description

M 1 [N*m] 1 N·m Unit moment

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 Select Boundaries 2–5 only.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 200e9 Pa Basic
Poisson's ratio nu 0.33 1 Basic
Density rho 7800 kg/m³ Basic

MESH 1
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Element size list, choose Extremely fine.
4 Click the Build All button.

13 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

STUDY 1

Step 1: Stationary
1 In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Define load cases check box.
4 Click Add three times.
5 In the table, enter the following settings:

6 In the Model Builder window, click Study 1.

7 In the Settings window for Study, type Study: Generalized Loads in the Label
text field.
8 On the Home toolbar, click Compute.

RESULTS

Derived Values
Create a table showing relation between gauge strain and each unit moment.

Point Evaluation 1
On the Results toolbar, click Point Evaluation.

Derived Values
1 Select Point 13 only.
2 In the Settings window for Point Evaluation, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Component
1>Shell>Strain>Strain tensor (local)>[Link] - Strain tensor (local), xx component.
3 Locate the Expression section. Select the Description check box.
4 In the associated text field, type e1.

14 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

5 Find the Parameters subsection. In the table, enter the following settings:

Name Value Unit Description

shell.z 1 Height of evaluation in shell, Height of evaluation in shell,
[-1,1] [-1,1]

6 Click the Evaluate button.

TA BL E
1 Go to the Table window.
2 Right-click Point Evaluation 1 and choose Duplicate.

RESULTS

Derived Values
1 In the Settings window for Point Evaluation, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Component
1>Shell>Strain>Strain tensor (local)>[Link] - Strain tensor (local), yy component.
2 Locate the Expression section. In the Description text field, type e2.
3 Click the Evaluate button.
4 In the Model Builder window, under Results>Derived Values right-click Point
Evaluation 1 and choose Duplicate.
5 Select Point 14 only.
6 In the Settings window for Point Evaluation, locate the Expression section.
7 In the Description text field, type e3.
8 Click the Evaluate button.

TA BL E
Go to the Table window.

GLOBAL DEFINITIONS
Load S-N curve.

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Click Browse.

15 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

5 Browse to the application’s Application Library folder and double-click the file
frame_with_cutout_SN_curve.txt.

6 Click Import.
Load twist and bending moments.

Interpolation 2 (int2)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Click Browse.
5 Browse to the application’s Application Library folder and double-click the file
frame_with_cutout_load.txt.

6 In the Number of arguments text field, type 1.

7 Click Import.
When an interpolation function is read from a file it treats the first columns as
arguments and the following one as function response. For example in the
interpolation function defining the moment around the y-axis the number of
arguments is assigned 1 and the function position is assigned 2. Thus the first
column is the argument and the third column is the function value (1+2).
8 Find the Functions subsection. In the table, enter the following settings:

Function name Position in file

fX 1
fY 2
fZ 3

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study: Generalized Loads

16 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue (ftg).
2 In the Settings window for Fatigue, type Fatigue Outside in the Label text field.
3 Right-click Component 1 (comp1)>Fatigue Outside and choose Edges>Cumulative
Damage.

FATIGUE OUTSIDE (FTG)

On the Physics toolbar, click Fatigue (ftg) and choose Fatigue Outside (ftg).

Cumulative Damage 1
1 Select all edges around the cutout.
2 In the Settings window for Cumulative Damage, locate the Solution Field section.
3 From the Physics interface list, choose Shell (shell).
4 Locate the Analysis section. From the Type list, choose Generalized loads.
5 Locate the Cycle Counting Parameters section. Find the Discretization subsection. In
the Nm text field, type 10.
6 In the Nr text field, type 20.
7 Locate the Damage Model Parameters section. In the m text field, type 10000.
8 From the σa(R,N) list, choose int1.
9 Locate the Generalized Load Definition section. In the sf text field, type 1000.
10 Click Add three times.
11 In the table, enter the following settings:

Generalized load history

fX
fY
fZ

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.

17 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Shell (shell)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2

Step 1: Fatigue
1 In the Model Builder window, under Study 2 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study: Generalized Loads, Stationary.
6 In the Model Builder window, click Study 2.
7 In the Settings window for Study, type Study: Fatigue Outside in the Label text
field.
8 On the Home toolbar, click Compute.

RESULTS

Fatigue Usage Factor (ftg)

Change the x-axis display to be an expression of the corner angles.

1 In the Model Builder window, expand the Fatigue Usage Factor (ftg) node, then click
Line Graph 1.
2 In the Settings window for Line Graph, locate the x-Axis Data section.
3 From the Parameter list, choose Expression.
4 In the Expression text field, type
atan2(z-0.03,x-0.04)*(dom==15)+atan2(z-0.03,x+0.04)*(dom==11)+ata
n2(z+0.03,x+0.04)*(dom==10)+atan2(z+0.03,x-0.04)*(dom==14).
5 From the Unit list, choose °.
6 Select the Description check box.

18 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

7 In the associated text field, type Cutout angle.

8 On the Fatigue Usage Factor (ftg) toolbar, click Plot.
9 In the Model Builder window, click Fatigue Usage Factor (ftg).
10 In the Settings window for 1D Plot Group, type Fatigue Usage Factor Outside
(ftg) in the Label text field.

Stress Cycle Distribution (ftg)

1 In the Model Builder window, under Results click Stress Cycle Distribution (ftg).
2 In the Settings window for 2D Plot Group, type Stress Cycle Distribution
Outside (ftg) in the Label text field.

Fatigue Usage Distribution (ftg)

1 In the Model Builder window, under Results click Fatigue Usage Distribution (ftg).
2 In the Settings window for 2D Plot Group, type Fatigue Usage Distribution
Outside (ftg) in the Label text field.

Examine fatigue at the bottom side of the shell elements.

SHELL (SHELL)
On the Physics toolbar, click Fatigue Outside (ftg) and choose Shell (shell).

1 In the Model Builder window, under Component 1 (comp1) click Shell (shell).
2 In the Settings window for Shell, click to expand the Height of evaluation in shell,
[-1,1] section.
3 Locate the Height of Evaluation in Shell, [-1,1] section. In the z text field, type -1.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study: Generalized Loads
Study: Fatigue Outside

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

19 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

FATIGUE 2 (FTG2)
1 In the Model Builder window, under Component 1 (comp1) click Fatigue 2 (ftg2).
2 In the Settings window for Fatigue, type Fatigue Inside in the Label text field.
3 Right-click Component 1 (comp1)>Fatigue Inside and choose Edges>Cumulative
Damage.

FATIGUE INSIDE (FTG2)

On the Physics toolbar, click Fatigue 2 (ftg2) and choose Fatigue Inside (ftg2).

Cumulative Damage 1
1 Select all edges around the cutout.
2 In the Settings window for Cumulative Damage, locate the Solution Field section.
3 From the Physics interface list, choose Shell (shell).
4 Locate the Analysis section. From the Type list, choose Generalized loads.
5 Locate the Cycle Counting Parameters section. Find the Discretization subsection. In
the Nm text field, type 10.
6 In the Nr text field, type 20.
7 Locate the Damage Model Parameters section. In the m text field, type 10000.
8 From the σa(R,N) list, choose int1.
9 Locate the Generalized Load Definition section. In the sf text field, type 1000.
10 Click Add three times.
11 In the table, enter the following settings:

Generalized load history

fX
fY
fZ

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.

20 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Shell (shell)
Fatigue Outside (ftg)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 3

Step 1: Fatigue
1 In the Model Builder window, under Study 3 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study: Generalized Loads, Stationary.
6 In the Model Builder window, click Study 3.
7 In the Settings window for Study, type Study: Fatigue Inside in the Label text
field.
8 On the Home toolbar, click Compute.

RESULTS

Fatigue Usage Factor (ftg2)

1 In the Model Builder window, expand the Fatigue Usage Factor (ftg2) node, then click
Line Graph 1.
2 In the Settings window for Line Graph, locate the x-Axis Data section.
3 From the Parameter list, choose Expression.
Define rounding position using an angle.
4 In the Expression text field, type
atan2(z-0.03,x-0.04)*(dom==15)+atan2(z-0.03,x+0.04)*(dom==11)+ata
n2(z+0.03,x+0.04)*(dom==10)+atan2(z+0.03,x-0.04)*(dom==14).
5 From the Unit list, choose °.
6 Select the Description check box.

21 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

7 In the associated text field, type Cutout angle.

8 On the Fatigue Usage Factor (ftg2) toolbar, click Plot.
9 In the Model Builder window, click Fatigue Usage Factor (ftg2).
10 In the Settings window for 1D Plot Group, type Fatigue Usage Factor Inside
(ftg2) in the Label text field.

Stress Cycle Distribution (ftg2)

1 In the Model Builder window, under Results click Stress Cycle Distribution (ftg2).
2 In the Settings window for 2D Plot Group, type Stress Cycle Distribution
Inside (ftg2) in the Label text field.

Fatigue Usage Distribution (ftg2)

1 In the Model Builder window, under Results click Fatigue Usage Distribution (ftg2).
2 In the Settings window for 2D Plot Group, type Fatigue Usage Distribution
Inside (ftg2) in the Label text field.

22 | FATIGUE RESPONSE OF A RANDOM NON-PROPORTIONAL LOAD

 Solved with COMSOL Multiphysics 5.2

Fatigue Analysis of a Wheel Rim

Introduction
During the development of safety critical components like a car wheel rim, making sure
fatigue cracks do not occur is one of the most important tasks. When a final prototype
is available this is ensured by testing, but prototype production and testing are time
consuming and expensive activities. Good predictions from simulations can keep down
the number the number of prototypes to a minimum.

In this example, you do a fatigue evaluation on a model of a wheel rim, subjected to

the load history from a simulated test.

Model Definition
For a definition of geometry, loads and boundary conditions, see the documentation
for the model Submodel in a Wheel Rim in the Structural Mechanics Module
Application Library.

The fatigue limit (in terms if the stress amplitude) is known for two cases with pure
axial loading. For pure tension it is 95 MPa, and for fully reversed loading it is
125 MPa. In this model, you use the Findley criterion, so the Findley parameters have
to be derived from these data.

In pure tension, the Findley criterion can be written as

 Δσ
2
------- + ( k ⋅ σ max ) + k ⋅ σ max = 2f
2
 2

This means that you have to solve the simultaneous equations

2 2
95 + ( k ⋅ 190 ) + k ⋅ 190 = 2f
2 2
125 + ( k ⋅ 125 ) + k ⋅ 125 = 2f

to get the Findley parameters f and k. The result is f = 84 MPa and k = 0.30.

Results and Discussion

The fatigue usage factor distribution is shown in Figure 1. The maximum value is
about 0.66, which should indicate that the design is good when taking into account

1 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

that the required safety factor has been included in the load. In Figure 2 the stress
histories in the critical point are displayed. The loading is slightly nonproportional, and
has a compressive mean stress, which is captured by the Findley criterion.

Figure 1: Fatigue usage factor using the Findley criterion.

2 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

Figure 2: Stress histories in the critical point.

Notes About the COMSOL Implementation

In this example, you perform the fatigue analysis as an additional study step in a model
which already contains the results from a stress analysis. Since the critical points for
fatigue crack initiation is on the free surface of the body, it is sufficient to do the fatigue
evaluation on the boundary, and not in the domain. This reduces CPU and memory
requirements significantly.

Application Library path: Fatigue_Module/Stress_Based/rim_fatigue

Modeling Instructions
In this example you will start from an existing model which is an example in the
Structural Mechanics Module.

From the File menu, choose Open.

3 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

Under the Application Library root, browse to the folder

Structural_Mechanics_Module/Tutorials and double-click the file
rim_submodel.mph.

RESULTS

Stress in Submodel
If the model was stored without solutions, you will now have to run Study 1 and Study
2 before continuing.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1
Study 2

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)
In the Model Builder window, expand the Results node.

Stress-Based 1
1 Right-click Component 2 (comp2)>Fatigue (ftg) and choose the boundary evaluation
Stress-Based.
2 Select Boundaries 2, 4, and 5 only.
3 In the Settings window for Stress-Based, locate the Solution Field section.
4 From the Physics interface list, choose Solid Mechanics 2 (solid2).

MATERIALS
In the Model Builder window, expand the Component 2 (comp2)>Materials node.

Material 3 (mat3)
1 Right-click Component 2 (comp2)>Materials and choose Blank Material.

4 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose All boundaries.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Normal stress sensitivity coefficient k_Findley 0.30 1 Findley
Limit factor f_Findley 84[MPa] Pa Findley

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)
Solid Mechanics 2 (solid2)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 3

Step 1: Fatigue
1 In the Model Builder window, under Study 3 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 2, Stationary.
6 On the Home toolbar, click Compute.

5 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

RESULTS

Fatigue Usage Factor (ftg)

1 In the Settings window for 3D Plot Group, locate the Plot Settings section.
2 From the View list, choose View 4.

Max/Min Surface 1
On the Fatigue Usage Factor (ftg) toolbar, click More Plots and choose Max/Min Surface.

Fatigue Usage Factor (ftg)

1 In the Settings window for Max/Min Surface, click Replace Expression in the
upper-right corner of the Expression section. From the menu, choose Component
2>Fatigue>[Link] - Fatigue usage factor.
2 Click to expand the Advanced section. From the Display list, choose Max.
3 On the Fatigue Usage Factor (ftg) toolbar, click Plot.

TABLE
1 Go to the Table window.
2 Click the Zoom Extents button on the Graphics toolbar.
In order to get the location of the point with maximum fatigue usage, zoom in on
the maximum marker and click on it in the graphics window. You will then see the
value and the coordinates in the Table window. The location will be approximately
(0.016, 0.092, 0.088). Use these coordinates to create a Cut Point 3D data set for
detailed evaluation of the stress history in the critical point.

RESULTS

Cut Point 3D 1
On the Results toolbar, click Cut Point 3D.

Data Sets
1 In the Settings window for Cut Point 3D, locate the Data section.
2 From the Data set list, choose Study 2/Solution 2 (3) (sol2).
3 Locate the Point Data section. In the X text field, type 0.0164.
4 In the Y text field, type 0.0924.
5 In the Z text field, type 0.0884.
6 Select the Snap to closest boundary check box.

6 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

1D Plot Group 4
On the Results toolbar, click 1D Plot Group.

Point Graph 1
On the 1D Plot Group 4 toolbar, click Point Graph.

1D Plot Group 4
1 In the Settings window for Point Graph, locate the Data section.
2 From the Data set list, choose Cut Point 3D 1.
3 Click Replace Expression in the upper-right corner of the y-axis data section. From
the menu, choose Component 2>Solid Mechanics 2>Stress>Stress tensor
(Spatial)>[Link] - Stress tensor, x component.
4 Locate the y-Axis Data section. From the Unit list, choose MPa.
5 Click to expand the Coloring and style section. Locate the Coloring and Style section.
Find the Line style subsection. From the Line list, choose Cycle.
6 Locate the x-Axis Data section. From the Axis source data list, choose All solutions.
7 Locate the Coloring and Style section. In the Width text field, type 3.
8 Click to expand the Legends section. Select the Show legends check box.
9 From the Legends list, choose Manual.
10 In the table, enter the following settings:

Legends
sx

11 Right-click Point Graph 1 and choose Duplicate.

12 In the Settings window for Point Graph, locate the y-Axis Data section.
13 In the Expression text field, type [Link].
14 Locate the Legends section. In the table, enter the following settings:

Legends
sy

15 Right-click Results>1D Plot Group 4>Point Graph 2 and choose Duplicate.

16 In the Settings window for Point Graph, locate the y-Axis Data section.
17 In the Expression text field, type [Link].

7 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

18 Locate the Legends section. In the table, enter the following settings:

Legends
sz

19 Right-click Results>1D Plot Group 4>Point Graph 3 and choose Duplicate.

20 In the Settings window for Point Graph, locate the y-Axis Data section.
21 In the Expression text field, type [Link].
22 Locate the Legends section. In the table, enter the following settings:

Legends
sxy

23 Right-click Results>1D Plot Group 4>Point Graph 4 and choose Duplicate.

24 In the Settings window for Point Graph, locate the y-Axis Data section.
25 In the Expression text field, type [Link].
26 Locate the Legends section. In the table, enter the following settings:

Legends
syz

27 Right-click Results>1D Plot Group 4>Point Graph 5 and choose Duplicate.

28 In the Settings window for Point Graph, locate the y-Axis Data section.
29 In the Expression text field, type [Link].
30 Locate the Legends section. In the table, enter the following settings:

Legends
sxz

31 In the Model Builder window, click 1D Plot Group 4.

32 In the Settings window for 1D Plot Group, click to expand the Title section.
33 From the Title type list, choose Manual.
34 In the Title text area, type Stress history in critical point.
35 On the 1D Plot Group 4 toolbar, click Plot.

8 | FATIGUE ANALYSIS OF A WHEEL RIM

 Solved with COMSOL Multiphysics 5.2

Fatigue Analysis of a
Non-Proportionally Loaded Shaft
with a Fillet
Introduction
This benchmark model is based on the example found in section 5.4.3 of Ref. 1. It
shows how to perform a high-cycle fatigue analysis for non-proportional loading using
critical plane methods.

Model Definition
The geometry is a circular shaft with two different diameters, 10 mm and 16 mm. At
the transition between the two diameters there is a fillet with a radius of 2 mm.

Figure 1: The notched shaft.

Two time-dependent loads are applied at the small end of the shaft: a transverse force
and a twisting moment. The force varies between 0 and 1.94 kN and the torque varies

1 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

between −28.7 and +28.7 Nm. Figure 2 shows the history of one loading cycle.

Figure 2: Load history.

The big end of the shaft is fixed. The material is Elastic with E = 100 GPa and ν = 0.

In Ref. 1 it is stated that the fatigue limit for completely reversed axial tension is
700 MPa, while the fatigue limit for pure tension is 560 MPa. These values are the
stress amplitudes. In uniaxial loading, the Findley criterion can be written as

2 2
( σ a ) + ( k ⋅ σ max ) + k ⋅ σ max = 2f

where σa is the stress amplitude and σmax is the maximum stress experienced in a
fatigue cycle. This means that you have to solve the simultaneous equations

2 2
700 + ( k ⋅ 700 ) + k ⋅ 700 = 2f
2 2
560 + ( k ⋅ 1120 ) + k ⋅ 1120 = 2f

to get the Findley parameters f and k. The result is f = 440 MPa and k = 0.23.

The Matake criterion is similar to Findley criterion, with the difference that the critical
plane is defined solely by the maximum shear stress. For a uniaxial case, the Matake
expression is

2 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

σa  kσ max
----- 1 + --------------
- = f
2  σa 

which gives the corresponding system of equations as

350 ⋅ ( 1 + k ) = f
280 ⋅ ( 1 + 2k ) = f

The solution is f = 467 MPa and k = 0.33 as parameters for the Matake case.

Results and Discussion

Figure 3 and Figure 4 show the stress distribution from the two basic load cases. The
location for the maximum effective stress is at the surface of the fillet, at a radius
slightly larger than the minimum radius of the shaft.

In Figure 5 the effective stress from the combined load case with transverse force and
positive torque is shown. It is symmetric with respect to the XY-plane, and is identical
also for the case when the torque is reversed.

3 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Figure 3: Axial stress from transverse force.

Figure 4: Effective stress from torque.

4 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Figure 5: Effective stress distribution for one of the combined load cases.

The results from the fatigue evaluation is shown in Figure 6 and Figure 7. With the
Findley criterion, the fatigue usage factor is computed to 0.98, in perfect agreement
with Ref. 1.

There is a large difference in the fatigue usage factor between the top and bottom side
of the bar, even though the effective stress is the same at both positions. This shows
how the criterion captures the difference between the predominantly tensile stress
states at the critical spot, and the compressive stress states on the other side.

Using the Matake criterion the fatigue usage factor decreases to 0.90, which shows that
there can be significant differences between results from seemingly similar models. The
critical plane computed in the Matake model differs from the one used in the Findley
model. As a consequence, the maximum normal stress on the critical plane can
significantly differ in the Matake case as compared to Findley case.

5 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Figure 6: Fatigue usage factor using the Findley criterion.

Figure 7: Fatigue usage factor using the Matake criterion.

6 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Notes About the COMSOL Implementation

In this model, you use the load case functionality in COMSOL to produce the load
cycle. In the first study the two basic load cases are analyzed. This study is not essential
for the analysis, but it allows you to inspect the results of the individual basic load cases.

Reference
1. D.F. Socie and G.B. Marquis, Multiaxial Fatigue, SAE, 1999.

Application Library path: Fatigue_Module/Stress_Based/shaft_with_fillet

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 3D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1
1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Work Plane 1 (wp1)

On the Geometry toolbar, click Work Plane.

7 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Bézier Polygon 1 (b1)

1 On the Geometry toolbar, click Primitives and choose Bézier Polygon.
2 In the Settings window for Bézier Polygon, locate the Polygon Segments section.
3 Find the Added segments subsection. Click Add Linear.
4 Find the Control points subsection. In row 2, set yw to 5.
5 Find the Added segments subsection. Click Add Linear.
6 Find the Control points subsection. In row 2, set xw to 30.
7 Find the Added segments subsection. Click Add Quadratic.
8 Find the Control points subsection. In row 2, set xw to 32.
9 In row 3, set xw to 32.
10 In row 3, set yw to 7.
11 Find the Added segments subsection. Click Add Linear.
12 Find the Control points subsection. In row 2, set yw to 8.
13 Find the Added segments subsection. Click Add Linear.
14 Find the Control points subsection. In row 2, set xw to 50.
15 Find the Added segments subsection. Click Add Linear.
16 Find the Control points subsection. In row 2, set yw to 0.
17 Right-click Bézier Polygon 1 (b1) and choose Build Selected.
18 Click the Zoom Extents button on the Graphics toolbar.

Revolve 1 (rev1)
1 On the Geometry toolbar, click Revolve.
2 In the Settings window for Revolve, locate the Revolution Angles section.
3 Clear the Keep original faces check box.
4 Locate the Revolution Axis section. Find the Direction of revolution axis subsection.
In the xw text field, type 1.
5 In the yw text field, type 0.
6 Right-click Revolve 1 (rev1) and choose Build Selected.
7 Click the Zoom Extents button on the Graphics toolbar.

SOLID MECHANICS (SOLID)

Fixed Constraint 1
1 On the Physics toolbar, click Boundaries and choose Fixed Constraint.

8 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

2 Select Boundaries 21–24 only.

Rigid Connector 1
1 On the Physics toolbar, click Boundaries and choose Rigid Connector.
2 Select Boundaries 1, 3, 5, and 7 only.

Applied Force 1
1 Right-click Rigid Connector 1 and choose Applied Force.
2 In the Settings window for Applied Force, locate the Applied Force section.
3 Specify the F vector as

0 x
0 y
-1.94[kN] z

4 On the Physics toolbar, click Load Group and choose New Load Group.

Applied Moment 1
1 In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid)
right-click Rigid Connector 1 and choose Applied Moment.
2 In the Settings window for Applied Moment, locate the Applied Moment section.
3 Specify the M vector as

28.7[N*m] x
0 y
0 z

4 On the Physics toolbar, click Load Group and choose New Load Group.

GLOBAL DEFINITIONS
1 In the Model Builder window, expand the Global Definitions node, then click Load
Group 1 (lg1).
2 In the Settings window for Load Group, type Transverse force in the Label text
field.
3 In the Parameter name text field, type lgF.
4 In the Model Builder window, under Global Definitions click Load Group 2 (lg2).
5 In the Settings window for Load Group, type Twisting moment in the Label text
field.

9 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

6 In the Parameter name text field, type lgM.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 100[GPa] Pa Basic
Poisson's ratio nu 0 1 Basic
Density rho 0 kg/m³ Basic

MESH 1
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Element size list, choose Fine.
4 Click the Build All button.
A finer mesh is needed in the fillet to resolve the stress concentration.
5 From the Sequence type list, choose User-controlled mesh.

Size 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Free
Tetrahedral 1 and choose Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Finer.

Size 2
1 Right-click Free Tetrahedral 1 and choose Size.
2 In the Settings window for Size, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Edge.
4 Select Edges 13, 14, 16, and 18 only.
5 Locate the Element Size section. Click the Custom button.

10 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

6 Locate the Element Size Parameters section. Select the Maximum element size check
box.
7 In the associated text field, type 0.5.
8 Select the Maximum element growth rate check box.
9 In the associated text field, type 1.2.
10 Click the Build All button.

STUDY 1

Step 1: Stationary
1 In the Model Builder window, under Study 1 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Define load cases check box.
4 Click Add two times.
5 In the table, enter the following settings:

6 In the Model Builder window, click Study 1.

7 In the Settings window for Study, type Study 1 (Basic load cases) in the Label
text field.
8 On the Home toolbar, click Compute.

RESULTS

Stress (solid)
Visualize the difference between the tension and compression on the opposite sides of
the shaft.

1 On the Stress (solid) toolbar, click Plot.

2 In the Settings window for 3D Plot Group, locate the Data section.
3 From the Load case list, choose Transverse force.
4 In the Model Builder window, expand the Stress (solid) node, then click Surface 1.

11 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

5 In the Settings window for Surface, click Replace Expression in the upper-right corner
of the Expression section. From the menu, choose Component 1>Solid
Mechanics>Stress>Stress tensor (Spatial)>[Link] - Stress tensor, x component.
6 On the Stress (solid) toolbar, click Plot.

ADD STUDY
1 On the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select study tree, select Preset Studies>Stationary.
4 Click Add Study in the window toolbar.
5 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 2

Step 1: Stationary
1 In the Model Builder window, under Study 2 click Step 1: Stationary.
2 In the Settings window for Stationary, locate the Study Extensions section.
3 Select the Define load cases check box.
4 Click Add three times.
5 In the table, enter the following settings:

6 In the Model Builder window, click Study 2.

7 In the Settings window for Study, type Study 2 (Combined load cases) in the
Label text field.
8 On the Home toolbar, click Compute.

RESULTS

Stress (solid) 1
1 On the Stress (solid) 1 toolbar, click Plot.
As a last step perform a fatigue analysis on the load cycle.

12 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1 (Basic load cases)
Study 2 (Combined load cases)

5 Click Add to Component in the window toolbar.

FATIGUE (FTG)
In the Model Builder window, expand the Stress (solid) 1 node.

Stress-Based 1
1 Right-click Component 1 (comp1)>Fatigue (ftg) and choose the boundary evaluation
Stress-Based.
2 In the Settings window for Stress-Based, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.
4 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
5 Locate the Evaluation Settings section. Find the Critical plane settings subsection. In
the Q text field, type 16.

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Recently Used>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Study 1 (Basic load cases)
Study 2 (Combined load cases)

4 Click Add to Component in the window toolbar.

13 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

5 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 2 (FTG2)
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue 2 (ftg2)
and choose the boundary evaluation Stress-Based.
2 In the Settings window for Stress-Based, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.
4 Locate the Fatigue Model Selection section. From the Criterion list, choose Matake.
5 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
6 Locate the Evaluation Settings section. Find the Critical plane settings subsection. In
the Q text field, type 16.

MATERIALS
Because the fatigue model is active only on the boundaries, you need to define a
material on the boundaries.

Material 2 (mat2)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 From the Geometric entity level list, choose Boundary.
4 From the Selection list, choose All boundaries.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Normal stress sensitivity k_Findley 0.23 1 Findley
coefficient
Limit factor f_Findley 440[MPa] Pa Findley
Normal stress sensitivity k_Matake 0.33 1 Matake
coefficient
Limit factor f_Matake 467[MPa] Pa Matake

ROOT
On the Home toolbar, click Windows and choose Add Study.

14 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

ADD STUDY
1 Go to the Add Study window.
2 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

3 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
4 Click Add Study in the window toolbar.

STUDY 3

Step 1: Fatigue
1 In the Model Builder window, under Study 3 click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Study 2 (Combined load cases), Stationary.
6 In the Model Builder window, click Study 3.
7 In the Settings window for Study, type Study 3 (Fatigue) in the Label text field.
8 On the Home toolbar, click Compute.

15 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

16 | FATIGUE ANALYSIS OF A NON-PROPORTIONALLY LOADED SHAFT WITH A FILLET

 Solved with COMSOL Multiphysics 5.2

Thermal Fatigue of a Surface Mount

Resistor
Introduction
A surface mount resistor is subjected to thermal cycling. The difference in the thermal
expansion between the materials introduces thermal stresses in the structure. The
solder, connecting the resistor to the printed circuit board, is seen as the weakest link
in the assembly. Because the operating temperature is high when compared to the
melting point of the solder, creep deformation occurs. In order to assure the structural
integrity of the component a fatigue analysis is performed where the life prediction
from two different fatigue models is compared.

Model Definition
A resistor is fastened on a printed circuit board (PCB) with SnAgCu solder. The solder
is connected to the printed circuit board through two copper pads and to the resistor
through a NiCr conductor. In reality there are additional thin films around the resistor
but they are disregarded in current analysis. A sketch of the surface mount assembly is
shown in Figure 1.

NiCr conductor Resistor

Solder
Copper pad

Printed circuit board

Figure 1: Schematic description of the surface mount resistor.

The resistor is made out of alumina and has dimensions 3.2 mm x 0.55 mm. It is
covered on both edges with a 0.025 mm layer of NiCr conductor. The thin layer
continues 0.325 mm along the lower and the upper side of the resistor. The printed
circuit board is large in comparison with the resistor and is here modeled 0.8 mm
thick. It has two copper pads on the top side that are 0.025 mm thick and 1.05 mm
wide. The thickness of the solder fillet between the copper pads and the NiCr

1 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

conductor is 0.05 mm. The remaining shape of the solder joint varies greatly between
each examined solder joint and is here modeled with two representative roundings.

Because the out-of-plane dimensions is 1.55 mm, which is significant in comparison

with the size of the resistor, the model is simulated in 2D with plane strain conditions.

The elastic properties of the materials are summarized in Table 1.

TABLE 1: ELASTIC AND THERMAL MATERIAL PROPERTIES

MATERIAL YOUNG’S MODULUS (GPA) POISSON’S RATIO COEFFICIENT OF THERMAL

EXPANSION (PPM/°C)

PCB laminate 22 0.4 21

Copper 141 0.35 17
SnAgCu 50 0.4 21
NiCr 170 0.31 13
Alumina 300 0.22 8

The SnAgCu solder material exhibits creep behavior, which can be modeled by a
Garofalo model where creep rate is described with

⋅ 10
4
c 6.19 –  5.32
-------------------------
dε ij 5  σe   RT  3 s ij
= 2.62 ⋅ 10 sinh  -------------------------6- e ⋅ --- ------ (1)
dT  39.1 ⋅ 10  2 σe

c
where ε ij is the creep strain tensor, T is the temperature, σe is the effective stress, R is
the universal gas constant, and sij is the deviatoric stress tensor.

The thermal load during an operating cycle is prescribed as a temperature which varies
between 20 °C and 70 °C. Each temperature change takes 2 minutes and is followed
by a 3 minutes dwell. This means that one fatigue cycle requires 10 minutes, see

2 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Figure 2.

Figure 2: Temperature load.

Since the solder material is nonlinear, see Equation 1, several cycles may need to
simulated before a stable cycle is obtained.

Two fatigue models are evaluated. The first is a strain-based Coffin-Manson type
model with the effective creep strain as the damage controlling mechanism and an
energy-based Morrow type model with the dissipated creep energy as the damage
controlling mechanism. The material constants for the Coffin-Manson model are
εf’=0.281 and c=-0.51. The material constants for the Morrow type model are
Wf’=55.0 J/m3 and m=-0.69.

Results and Discussion

The difference in the elastic and thermal properties introduces thermal stresses in the
device. Although they are not very high the solder experiences significant inelastic

3 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

strains. In Figure 3 accumulated effective creep strain after six cycles is shown.

Figure 3: Creep strains in the solder joint.

The highest strains occur in the thin solder layer just below the resistor. It is mainly the
shear strain component which contributes to the effective creep strain in that layer. The
slightly higher values around the edge are also affected by modeling of a sharp corner.
With a rounding instead the strains are somewhat lower. Nevertheless the location of
highest strain agrees well with the crack path in real applications. In order to evaluate
fatigue it is important to obtain a stable load cycle. In applications involving solder
joints, frequently either inelastic strain or dissipated energy is used to predict fatigue.
The change of creep strain during the first six cycles is therefore evaluated in a point
just below the resistor slightly shifted to the right form the sharp corner. The position
of this point can be debated. It is however located in the area where the largest strains
occurs and is therefore seen as the critical point. In Figure 4 the effective creep strain
and the shear creep strain component are shown.

4 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Figure 4: Creep strain development in a critical point below the resistor.

The dissipated energy represents a combined contribution of changes in stresses and

strains during a cycle and is in Figure 5 shown with a shear hysteresis. The shear
component has been chosen since it gives the dominating contribution to the effective
creep strain in the critical point.

5 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Figure 5: Shear hysteresis evaluated in a critical point just below the resistor.

It is clear from the two last figures that the first cycle is not representative for fatigue
analysis, since its response differs significantly from the one experienced in the
following cycles. Even after six cycles, the stress-strain loop has not stabilized. The
temperature cycling can by extended with additional cycles to evaluate whether the
state stabilizes further or not. In some structures it can be so that the hysteresis loop
is moving in stress-strain space. In this example additional cycles are not simulated
since the difference in the creep strain and the dissipated energy between cycle five and
six is small. Assuming that the consecutive cycles follow the trend and deform less as
well as dissipate less energy, the fatigue analysis based on the results of the sixth cycle
gives a conservative fatigue prediction.

The fatigue life based on the Coffin-Manson model is shown in Figure 6, and the
fatigue life based on the Morrow model is shown in Figure 7.

6 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Figure 6: Fatigue life based on the creep strain.

7 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Figure 7: Fatigue life based on the dissipated energy.

Fatigue based on strain gives lifetime of 1805 cycles while the energy prediction gives
2563 cycles.

Notes About the COMSOL Implementation

The shear strain in the thin section between the PCB and the resistor gives a
dominating contribution to creep. Therefor it might be so that it instead should be
evaluated in the fatigue expression of the Coffin-Manson model. This is possible using
the User defined option in the Strain type parameter. Type then the expression
solid.ec12 in the field for the Inelastic strain.

Application Library path: Fatigue_Module/Energy_Based/

surface_mount_resistor

8 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select physics tree, select Structural Mechanics>Solid Mechanics (solid).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Time Dependent.
6 Click Done.

GEOMETRY 1
Begin by changing the length unit to millimeters.

1 In the Model Builder window, under Component 1 (comp1) click Geometry 1.

2 In the Settings window for Geometry, locate the Units section.
3 From the Length unit list, choose mm.

Import 1 (imp1)
1 On the Home toolbar, click Import.
2 In the Settings window for Import, locate the Import section.
3 Click Browse.
4 Browse to the application’s Application Library folder and double-click the file
surface_mount_resistor.mphbin.
5 Click Import.

GLOBAL DEFINITIONS

Interpolation 1 (int1)
1 On the Home toolbar, click Functions and choose Global>Interpolation.
2 In the Settings window for Interpolation, locate the Definition section.
3 From the Data source list, choose File.
4 Click Browse.

9 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

5 Browse to the application’s Application Library folder and double-click the file
surface_mount_resistor_thermal_load_cycle.txt.

6 Click Import.
7 In the Function name text field, type thermLC.
8 Locate the Units section. In the Arguments text field, type min.
9 In the Function text field, type degC.

SOLID MECHANICS (SOLID)

1 In the Model Builder window, under Component 1 (comp1) click Solid Mechanics
(solid).
2 In the Settings window for Solid Mechanics, locate the Thickness section.
3 In the d text field, type 1.55 [mm].

Linear Elastic Material 1

1 In the Model Builder window’s toolbar, click the Show button and select Advanced
Physics Options in the menu.
2 In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid)
click Linear Elastic Material 1.
3 In the Settings window for Linear Elastic Material, click to expand the Energy
dissipation section.
4 Locate the Energy Dissipation section. Select the Calculate dissipated energy check
box.

Thermal Expansion 1
1 On the Physics toolbar, click Attributes and choose Thermal Expansion.
2 In the Settings window for Thermal Expansion, locate the Model Inputs section.
3 In the T text field, type thermLC(t).
4 Locate the Thermal Expansion Properties section. In the Tref text field, type 20
[degC].

Linear Elastic Material 1

In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click
Linear Elastic Material 1.

Creep 1
1 On the Physics toolbar, click Attributes and choose Creep.
2 In the Settings window for Creep, locate the Domain Selection section.

10 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

3 Click Clear Selection.

4 Select Domain 3 only.
5 Locate the Model Inputs section. In the T text field, type thermLC(t).
6 Locate the Creep Data section. From the Material model list, choose Garofalo
(hyperbolic sine).
7 In the A text field, type 262000.
8 In the σref text field, type 39.1[MPa].
9 In the n text field, type 6.19.
10 Select the Include temperature dependency check box.
11 In the Q text field, type 53200.
12 Click the Zoom Extents button on the Graphics toolbar.

Symmetry 1
1 On the Physics toolbar, click Boundaries and choose Symmetry.
2 Select Boundaries 19 and 20 only.

Roller 1
1 On the Physics toolbar, click Boundaries and choose Roller.
2 Select Boundary 2 only.

MATERIALS

Material 1 (mat1)
1 In the Model Builder window, under Component 1 (comp1) right-click Materials and
choose Blank Material.
2 In the Settings window for Material, type PCB in the Label text field.
3 Locate the Geometric Entity Selection section. Click Clear Selection.
4 Select Domain 1 only.
5 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 22 [GPa] Pa Basic
Poisson's ratio nu 0.4 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 21e-6 1/K Basic

11 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Material 2 (mat2)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Copper in the Label text field.
3 Select Domain 2 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 141 [GPa] Pa Basic
Poisson's ratio nu 0.35 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 16.6e-6 1/K Basic

Material 3 (mat3)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Solder in the Label text field.
3 Select Domain 3 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 50 [GPa] Pa Basic
Poisson's ratio nu 0.4 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 21e-6 1/K Basic

Material 4 (mat4)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type NiCr in the Label text field.
3 Select Domain 4 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 170 [GPa] Pa Basic
Poisson's ratio nu 0.31 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 13e-6 1/K Basic

12 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Material 5 (mat5)
1 Right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Alumina in the Label text field.
3 Select Domain 5 only.
4 Locate the Material Contents section. In the table, enter the following settings:

Property Name Value Unit Property group

Young's modulus E 300 [GPa] Pa Basic
Poisson's ratio nu 0.22 1 Basic
Density rho 1 kg/m³ Basic
Coefficient of thermal expansion alpha 8e-6 1/K Basic

MESH 1

Free Triangular 1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Free Triangular.

Distribution 1
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 right-click Free
Triangular 1 and choose Distribution.
2 Select Boundary 8 only.
3 In the Settings window for Distribution, locate the Distribution section.
4 In the Number of elements text field, type 30.

Size
1 In the Model Builder window, under Component 1 (comp1)>Mesh 1 click Size.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Finer.
4 Click the Build All button.

STUDY 1

Step 1: Time Dependent

Simulate a time history of 6 cycles.

1 In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.

13 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

3 In the Times text field, type range(0,10,60*60).

4 In the Model Builder window, click Study 1.
5 In the Settings window for Study, type Time history in the Label text field.

Solution 1 (sol1)
On the Study toolbar, click Show Default Solver.

TIME HISTORY

Solution 1 (sol1)
Force strict time stepping in order to improve the creep results.

1 In the Model Builder window, expand the Solution 1 (sol1) node, then click
Time-Dependent Solver 1.
2 In the Settings window for Time-Dependent Solver, click to expand the Time
stepping section.
3 Locate the Time Stepping section. From the Steps taken by solver list, choose Strict.
4 On the Study toolbar, click Compute.

RESULTS

Stress (solid)
1 In the Model Builder window, expand the Stress (solid) node.
2 In the Model Builder window, expand the Results>Stress (solid)>Surface 1 node.
3 Right-click Deformation and choose Disable.
4 Click the Zoom Extents button on the Graphics toolbar.
Display the accumulated creep strain in the solder joint.

Stress (solid) 1
1 Right-click Stress (solid) and choose Duplicate.
2 In the Settings window for 2D Plot Group, type Creep strain distribution in
the Label text field.

Creep strain distribution

1 In the Model Builder window, expand the Results>Creep strain distribution node, then
click Surface 1.
2 In the Settings window for Surface, click Replace Expression in the upper-right corner
of the Expression section. From the menu, choose Component 1>Solid
Mechanics>Strain (Gauss points)>[Link] - Effective creep strain.

14 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

3 Click the Zoom Extents button on the Graphics toolbar.

Display creep strain history. The shear component gives the largest contribution to
the effective creep strain.

1D Plot Group 3
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Creep strain in the Label text
field.

Point Graph 1
On the Creep strain toolbar, click Point Graph.

Creep strain
1 Select Point 7 only.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid
Mechanics>Strain (Gauss points)>[Link] - Effective creep strain.
3 Click to expand the Legends section. Select the Show legends check box.
4 From the Legends list, choose Manual.
5 In the table, enter the following settings:

Legends
Effective

6 Right-click Point Graph 1 and choose Duplicate.

7 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid
Mechanics>Strain (Gauss points)>Creep strain tensor, local coordinate
system>solid.ecGp12 - Creep strain tensor, local coordinate system, 12 component.
8 Locate the Legends section. In the table, enter the following settings:

Legends
Shear

9 In the Model Builder window, click Creep strain.

10 In the Settings window for 1D Plot Group, click to expand the Legend section.
11 From the Position list, choose Upper left.
12 On the Creep strain toolbar, click Plot.

15 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

13 Click the Zoom Extents button on the Graphics toolbar.

Display stress-strain hysteresis of the shear behavior.

1D Plot Group 4
1 On the Home toolbar, click Add Plot Group and choose 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Shear hysteresis in the Label
text field.

Point Graph 1
On the Shear hysteresis toolbar, click Point Graph.

Shear hysteresis
1 Select Point 7 only.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Component 1>Solid
Mechanics>Stress (Gauss points)>Stress tensor, Gauss-point evaluation
(Spatial)>[Link] - Stress tensor, Gauss-point evaluation, xy component.
3 Locate the x-Axis Data section. From the Parameter list, choose Expression.
4 Click Replace Expression in the upper-right corner of the x-axis data section. From
the menu, choose Component 1>Solid Mechanics>Strain (Gauss points)>Creep strain
tensor, local coordinate system>solid.ecGp12 - Creep strain tensor, local coordinate
system, 12 component.
5 On the Shear hysteresis toolbar, click Plot.
6 Click the Zoom Extents button on the Graphics toolbar.

COMPONENT 1 (COMP1)
On the Home toolbar, click Windows and choose Add Physics.

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Time history

4 Click Add to Component in the window toolbar.

16 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

FATIGUE (FTG)

Strain-Life 1
1 In the Model Builder window, right-click Fatigue (ftg) and choose the domain
evaluation Strain-Life.
2 Select Domain 3 only.
3 In the Settings window for Strain-Life, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Coffin-Manson.
5 From the Strain type list, choose Effective creep strain.
6 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

ADD PHYSICS
1 Go to the Add Physics window.
2 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Studies Solve
Time history

4 Click Add to Component in the window toolbar.

5 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 2 (FTG2)

Energy-Based 1
1 In the Model Builder window, under Component 1 (comp1) right-click Fatigue 2 (ftg2)
and choose the domain evaluation Energy-Based.
2 Select Domain 3 only.
3 In the Settings window for Energy-Based, locate the Fatigue Model Selection section.
4 From the Energy type list, choose Creep dissipation density.
5 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).

17 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

MATERIALS

Solder (mat3)
1 In the Model Builder window, expand the Component 1 (comp1)>Materials>Solder
(mat3) node, then click Solder (mat3).
2 In the Settings window for Material, click to expand the Material properties section.
3 Locate the Material Properties section. In the Material properties tree, select Solid
Mechanics>Fatigue Behavior>Energy-Based>Morrow.
4 Click Add to Material.
5 In the Material properties tree, select Solid Mechanics>Fatigue
Behavior>Strain-Based>Coffin-Manson.
6 Click Add to Material.
7 In the Model Builder window, under Component 1 (comp1)>Materials>Solder (mat3)
click Morrow (fatigueEnergyMorrow).
8 In the Settings window for Property Group, locate the Output Properties and Model
Inputs section.
9 Find the Output properties subsection. In the table, enter the following settings:

Property Variable Expression Unit Size

Fatigue energy coefficient Wf_Morrow 55e6 [J/m^3] J/m³ 1x1
Fatigue energy exponent m_Morrow -0.69 1 1x1

10 In the Model Builder window, under Component 1 (comp1)>Materials>Solder (mat3)

click Coffin-Manson (fatigueStrainCoffinManson).
11 In the Settings window for Property Group, locate the Output Properties and Model
Inputs section.
12 Find the Output properties subsection. In the table, enter the following settings:

Property Variable Expression Unit Size

Fatigue ductility coefficient epsilonf_CM 0.218 1 1x1
Fatigue ductility exponent c_CM -0.51 1 1x1

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.

18 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Physics Solve
Solid Mechanics (solid)

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.

STUDY 2
1 In the Model Builder window, click Study 2.
2 In the Settings window for Study, type Fatigue in the Label text field.

FATIGUE

Step 1: Fatigue
1 In the Model Builder window, under Fatigue click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Time history, Time Dependent.
6 From the Time (s) list, choose From list.
7 From the list select time steps from 3000s to 3600s.
8 On the Home toolbar, click Compute.

19 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

20 | THERMAL FATIGUE OF A SURFACE MOUNT RESISTOR

 Solved with COMSOL Multiphysics 5.2

Energy-Based Thermal Fatigue

Prediction in a Ball Grid Array
Introduction
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 results a crack grows through a solder joint and disconnects the
chip from the printed circuit board. The microelectronic component loses then its
operational functionality. 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 transfered 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 and takes few
minutes.

1 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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/m2K. 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 Multiplier of stress
m 0.303 Strain rate sensitivity of stress
s0 80.42 MPa Coefficient for deformation resistance saturation
sinit 56.33 MPa Initial value of deformation resistance
h0 2640.75 MPa Hardening constant
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

2 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

ΔW ave k 2
N = K 1  --------------- + --------------------------------
a -
 W ref  Δ W k

K 3 ---------------
ave 4
 W ref 

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.92e-7 m Crack propagation energy coefficient
k4 1.12 Crack propagation energy exponent
Wref 689 J/m3 Reference energy density

The variable a can in this analysis be taken as the diameter of the plane where a crack
propagates. 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:

1 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.
2 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.
3 A submodel of the critical joint with a fine mesh is created. The results from step 1
are prescribed via Prescribed Displacement on the boundaries where the submodel
is cut out of the full model. The result of the thermal analysis, the temperature, is
prescribed directly to the viscoplastic material node and via Thermal Expansion to
the whole submodel. This is done for all time steps of the four simulated cycles.

3 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

Since the results of the thermal analysis are prescribed to the model in each step,
only a structural study is performed. This reduces the initially multiphysics model to
a single physics model.
4 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 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 of
different parts of the component. 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, 103.4 cycles.

4 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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.

5 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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 at capturing the strain gradient close to the upper side of the
solder joint.

6 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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 103.4 Figure 3

Entire joint in the submodel 103.5 Figure 6
Thin layer in the submodel 102.6 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 is created in all joints of the full model, the simulation requires significant
computational resources. In the current example the full model consists of about
3.1 × 105 kDOFs while the submodel consists of 1.1 × 105 kDOFs.

7 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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

8 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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 General Extrusion operator. In this
example the mapping of results is simple since it is 1-1. The General Extrusion
operator allows also for mapping into other shapes through analytical functions.

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 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, a sharp change in a load parameter is 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 COMSOL via flc2hs(t,p) function that is a
Heaviside function with a smooth second derivative. The first argument defines a

9 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

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.

Application Library path: 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.

From the File menu, choose Open.

Under the Application Library root, browse to the folder

Nonlinear_Structural_Materials_Module/Viscoplasticity and double-click
the file viscoplastic_solder_joints.mph.

COMPONENT 1 (COMP1)
1 In the Model Builder window, click Component 1 (comp1).
2 In the Settings window for Component, type Full Model in the Label text field.
Define the power load cycle.

GLOBAL DEFINITIONS

Analytic 2 (an2)
On the Home toolbar, click Functions and choose Global>Analytic.

Analytic 1 (power)
1 In the Settings window for Analytic, type power in the Function name text field.
2 Locate the Definition section. In the Expression text field, type
(flc2hs(x-5*60,5*60)*5e7)*(x<6*60*60)-flc2hs(mod(x,6*60*60)-(4*60
*60+5*60),5*60)*4e7+(flc2hs(mod(x,6*60*60)-5*60,5*60)*4e7+1e7)*(x
>=6*60*60).
3 Locate the Units section. In the Arguments text field, type s.
4 In the Function text field, type W/m^3.

FULL MODEL (COMP1)

In the Model Builder window, expand the Full Model (comp1) node.

10 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

H E A T TR A N S F E R I N S O L I D S ( H T )

Heat Source 1
1 In the Model Builder window, expand the Full Model (comp1)>Heat Transfer in Solids
(ht) node, then click Heat Source 1.
2 In the Settings window for Heat Source, locate the Heat Source section.
3 In the Q0 text field, type power(t).

DEFINITIONS

General Extrusion 1 (genext1)

1 On the Definitions toolbar, click Component Couplings and choose General Extrusion.
2 In the Settings window for General Extrusion, locate the Source Selection section.
3 From the Selection list, choose All domains.

STUDY 1
1 In the Model Builder window, click Study 1.
2 In the Settings window for Study, type Full Model: Load History in the Label
text field.

FULL MODEL: LOAD HISTORY

Step 1: Time Dependent

1 In the Model Builder window, expand the Full Model: Load History node, then click
Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, type 0 20 60 range(2*60,60,12*60)
range(15*60,3*60,30*60) range(40*60,10*60,3.5*60*60) 3.7*3600
3.9*3600 4*3600 4*3600+20 4*3600+60
range(4*3600+2*60,60,4*3600+12*60)
range(4*3600+15*60,3*60,4*3600+30*60)
range(4*3600+40*60,10*60,4*3600+1.5*60*60) 4*3600+1.7*3600
4*3600+1.9*3600 6*3600 6*3600+20 6*3600+60
range(6*3600+2*60,60,6*3600+12*60)
range(6*3600+15*60,3*60,6*3600+30*60)
range(6*3600+40*60,10*60,6*3600+3.5*60*60) 6*3600+3.7*3600
6*3600+3.9*3600 10*3600 10*3600+20 10*3600+60
range(10*3600+2*60,60,10*3600+12*60)

11 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

range(10*3600+15*60,3*60,10*3600+30*60)
range(10*3600+40*60,10*60,10*3600+1.5*60*60) 10*3600+1.7*3600
10*3600+1.9*3600 12*3600 12*3600+20 12*3600+60
range(12*3600+2*60,60,12*3600+12*60)
range(12*3600+15*60,3*60,12*3600+30*60)
range(12*3600+40*60,10*60,12*3600+3.5*60*60) 12*3600+3.7*3600
12*3600+3.9*3600 16*3600 16*3600+20 16*3600+60
range(16*3600+2*60,60,16*3600+12*60)
range(16*3600+15*60,3*60,16*3600+30*60)
range(16*3600+40*60,10*60,16*3600+1.5*60*60) 16*3600+1.7*3600
16*3600+1.9*3600 18*3600 18*3600+20 18*3600+60
range(18*3600+2*60,60,18*3600+12*60)
range(18*3600+15*60,3*60,18*3600+30*60)
range(18*3600+40*60,10*60,18*3600+3.5*60*60) 18*3600+3.7*3600
18*3600+3.9*3600 22*3600 22*3600+20 22*3600+60
range(22*3600+2*60,60,22*3600+12*60)
range(22*3600+15*60,3*60,22*3600+30*60)
range(22*3600+40*60,10*60,22*3600+1.5*60*60) 22*3600+1.7*3600
22*3600+1.9*3600 24*3600.
4 In the Model Builder window, expand the Full Model: Load History>Solver
Configurations node.

Step 2: Time Dependent 2

1 In the Model Builder window, expand the Full Model: Load History>Solver
Configurations>Solution 1 (sol1) node, then click Full Model: Load History>Step 2: Time
Dependent 2.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, enter the same steps as Step 1.
4 On the Home toolbar, click Compute.

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

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).

12 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

Click Add to Component in the window toolbar.

5 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE (FTG)

Energy-Based 1
1 In the Model Builder window, under Full Model (comp1) right-click Fatigue (ftg) and
choose the domain evaluation Energy-Based.
2 In the Settings window for Energy-Based, locate the Domain Selection section.
3 From the Selection list, choose Solder.
4 Locate the Fatigue Model Selection section. From the Criterion list, choose Darveaux.
5 From the Energy type list, choose Creep dissipation density.
6 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics (solid).
7 Locate the Fatigue Model Parameters section. In the a text field, type 2.6457e-4.

MATERIALS

Solder, 60Sn-40Pb (mat4)

1 In the Model Builder window, expand the Full Model (comp1)>Materials node, then
click Solder, 60Sn-40Pb (mat4).
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Crack initiation energy K1_Darveaux 13173 1 Darveaux
coefficient
Crack initiation energy k2_Darveaux -1.45 1 Darveaux
exponent
Crack propagation energy K3_Darveaux 3.92e-7 m Darveaux
coefficient

13 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

Property Name Value Unit Property group

Crack propagation energy k4_Darveaux 1.12 1 Darveaux
exponent
Reference energy density Wref_Darveaux 689 J/m³ Darveaux

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
5 Click Add Study in the window toolbar.
6 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 2
1 In the Model Builder window, click Study 2.
2 In the Settings window for Study, type Full Model: Fatigue Evaluation in the
Label text field.

FULL MODEL: FATIGUE EVALUATION

Step 1: Fatigue
1 In the Model Builder window, under Full Model: Fatigue Evaluation click Step 1:
Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Full Model: Load History, Time Dependent 2.

14 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

6 From the Time (s) list, choose From list.

7 In the Time (s) list, elect time steps from 6.48e4s to 8.64e4s.
8 On the Home toolbar, click Compute.

RESULTS

Cycles to Failure (ftg)

The solder joint with the shortest life is the critical one. Create a submodel of that
joint.

ROOT
1 In the Model Builder window, click the root node.
2 Click Add Component and choose 3D.

COMPONENT 2 (COMP2)
1 In the Model Builder window, click Component 2 (comp2).
2 In the Settings window for Component, type Submodel in the Label text field.

GEOMETRY 2

Import 1 (imp1)
1 On the Home toolbar, click Import.
2 In the Settings window for Import, locate the Import section.
3 Click Browse.
4 Browse to the application’s Application Library folder and double-click the file
viscoplastic_solder_joints.mphbin.

5 Click Import.

Cylinder 1 (cyl1)
1 On the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type 2.4e-4.
4 In the Height text field, type 9e-4.
5 Locate the Position section. In the x text field, type 68.28e-4.
6 In the y text field, type 55.38e-4.
7 In the z text field, type 10e-4.

15 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

Union 1 (uni1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Union.
2 Click the Select All button on the Graphics toolbar and remove the cylinder from the
selection.

Intersection 1 (int1)
1 On the Geometry toolbar, click Booleans and Partitions and choose Intersection.
2 Select the objects cyl1 and uni1 only.

Work Plane 1 (wp1)

1 On the Geometry toolbar, click Work Plane.
2 Click the Zoom Extents button on the Graphics toolbar.
3 In the Settings window for Work Plane, locate the Plane Definition section.
4 In the z-coordinate text field, type 13.5e-4.

Work Plane 2 (wp2)

1 Right-click Work Plane 1 (wp1) and choose Duplicate.
2 In the Settings window for Work Plane, locate the Plane Definition section.
3 In the z-coordinate text field, type 15.5e-4.

Partition Objects 1 (par1)

1 On the Geometry toolbar, click Booleans and Partitions and choose Partition Objects.
2 In the Settings window for Partition Objects, locate the Partition Objects section.
3 From the Partition with list, choose Work plane.
4 From the Work plane list, choose Work Plane 1 (wp1).
5 In the Relative repair tolerance text field, type 3E-4.
6 Select the object int1 only.

GEOMETRY 2

Partition Objects 2 (par2)

1 On the Geometry toolbar, click Booleans and Partitions and choose Partition Objects.
2 Select the object par1 only.
3 In the Settings window for Partition Objects, locate the Partition Objects section.
4 From the Partition with list, choose Work plane.
5 Click the Build All Objects button.
6 Click the Zoom Extents button on the Graphics toolbar.

16 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

7 In the Model Builder window, collapse the Geometry 2 node.

GEOMETRY 2
In the Model Builder window, collapse the Submodel (comp2)>Geometry 2 node.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Solid Mechanics (solid).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

SOLID MECHANICS 2 (SOLID2)

1 In the Model Builder window, under Submodel (comp2) click Solid Mechanics 2
(solid2).
2 In the Settings window for Solid Mechanics, locate the Structural Transient Behavior
section.
3 From the list, choose Quasi-static.

Linear Elastic Material 1

1 In the Model Builder window’s toolbar, click the Show button and select Advanced
Physics Options in the menu.
2 In the Model Builder window, under Submodel (comp2)>Solid Mechanics 2 (solid2)
click Linear Elastic Material 1.
3 In the Settings window for Linear Elastic Material, click to expand the Energy
dissipation section.
4 Locate the Energy Dissipation section. Select the Calculate dissipated energy check
box.

Viscoplasticity 1
1 On the Physics toolbar, click Attributes and choose Viscoplasticity.
2 In the Settings window for Viscoplasticity, locate the Domain Selection section.

17 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

3 Click Clear Selection.

4 Select Domains 4–6 only.
5 Locate the Model Inputs section. In the T text field, type
comp1.genext1(comp1.T).

6 Locate the Viscoplasticity Model section. In the A text field, type 1.49e7.
7 In the Q text field, type 90046.
8 In the ξ text field, type 11.
9 In the m text field, type 0.303.
10 In the s0 text field, type 80.42[MPa].
11 In the sinit text field, type 56.33[MPa].
12 In the h0 text field, type 2640.75[MPa].
13 In the a text field, type 1.34.
14 In the n text field, type 0.0231.

Linear Elastic Material 1

In the Model Builder window, under Submodel (comp2)>Solid Mechanics 2 (solid2) click
Linear Elastic Material 1.

Thermal Expansion 1
1 On the Physics toolbar, click Attributes and choose Thermal Expansion.
Prescribe results of the thermal analysis via thermal expansion in the entire
submodel.
2 In the Settings window for Thermal Expansion, locate the Model Inputs section.
3 In the T text field, type comp1.genext1(comp1.T).
4 Locate the Thermal Expansion Properties section. In the Tref text field, type T0.
Prescribe results of the structural analysis via displacements on the shared
boundaries.

Prescribed Displacement 1
1 On the Physics toolbar, click Boundaries and choose Prescribed Displacement.
2 Select Boundaries 1–5, 8, 9, 11, and 22–27 only.
3 In the Settings window for Prescribed Displacement, locate the Prescribed
Displacement section.
4 Select the Prescribed in x direction check box.
5 In the u0,x text field, type comp1.genext1(comp1.u).

18 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

6 Select the Prescribed in y direction check box.

7 In the u0,y text field, type comp1.genext1(comp1.v).
8 Select the Prescribed in z direction check box.
9 In the u0,z text field, type comp1.genext1(comp1.w).
10 Click to expand the Constraint settings section. Locate the Constraint Settings
section. From the Apply reaction terms on list, choose Current physics (internally
symmetric).

ADD MATERIAL
1 On the Home toolbar, click Add Material to open the Add Material window.
2 Go to the Add Material window.
3 In the tree, select Built-In>FR4 (Circuit Board).
4 Click Add to Component in the window toolbar.

MATERIALS

FR4 (Circuit Board) (mat5)

1 In the Model Builder window, under Submodel (comp2)>Materials click FR4 (Circuit
Board) (mat5).
2 In the Settings window for Material, locate the Geometric Entity Selection section.
3 Click Clear Selection.
4 Select Domain 1 only.

ADD MATERIAL
1 Go to the Add Material window.
2 In the tree, select Built-In>Copper.
3 Click Add to Component in the window toolbar.

MATERIALS

Copper (mat6)
1 In the Model Builder window, under Submodel (comp2)>Materials click Copper (mat6).
2 Select Domain 2 only.

ADD MATERIAL
1 Go to the Add Material window.
2 In the tree, select Built-In>Silicon.

19 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

3 Click Add to Component in the window toolbar.

MATERIALS

Silicon (mat7)
1 In the Model Builder window, under Submodel (comp2)>Materials click Silicon (mat7).
2 Select Domain 3 only.

ADD MATERIAL
1 Go to the Add Material window.
2 In the tree, select Built-In>Solder, 60Sn-40Pb.
3 Click Add to Component in the window toolbar.
4 On the Home toolbar, click Add Material to close the Add Material window.

MATERIALS

Solder, 60Sn-40Pb (mat8)

1 In the Model Builder window, under Submodel (comp2)>Materials click Solder,
60Sn-40Pb (mat8).
2 Select Domains 4–6 only.
Simulate the initial cycles using the submodel.

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.
2 Find the Studies subsection. In the Select study tree, select Preset Studies.
3 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

4 Find the Studies subsection. In the Select study tree, select Preset Studies>Time
Dependent.
5 Click Add Study in the window toolbar.

20 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

6 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 3
1 In the Model Builder window, click Study 3.
2 In the Settings window for Study, type Submodel: Load History in the Label text
field.

SUBMODEL: LOAD HISTORY

Step 1: Time Dependent

1 In the Model Builder window, click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, type 0 20 60 range(2*60,60,12*60)
range(15*60,3*60,30*60) range(40*60,10*60,3.5*60*60) 3.7*3600
3.9*3600 4*3600 4*3600+20 4*3600+60
range(4*3600+2*60,60,4*3600+12*60)
range(4*3600+15*60,3*60,4*3600+30*60)
range(4*3600+40*60,10*60,4*3600+1.5*60*60) 4*3600+1.7*3600
4*3600+1.9*3600 6*3600 6*3600+20 6*3600+60
range(6*3600+2*60,60,6*3600+12*60)
range(6*3600+15*60,3*60,6*3600+30*60)
range(6*3600+40*60,10*60,6*3600+3.5*60*60) 6*3600+3.7*3600
6*3600+3.9*3600 10*3600 10*3600+20 10*3600+60
range(10*3600+2*60,60,10*3600+12*60)
range(10*3600+15*60,3*60,10*3600+30*60)
range(10*3600+40*60,10*60,10*3600+1.5*60*60) 10*3600+1.7*3600
10*3600+1.9*3600 12*3600 12*3600+20 12*3600+60
range(12*3600+2*60,60,12*3600+12*60)
range(12*3600+15*60,3*60,12*3600+30*60)
range(12*3600+40*60,10*60,12*3600+3.5*60*60) 12*3600+3.7*3600
12*3600+3.9*3600 16*3600 16*3600+20 16*3600+60
range(16*3600+2*60,60,16*3600+12*60)
range(16*3600+15*60,3*60,16*3600+30*60)
range(16*3600+40*60,10*60,16*3600+1.5*60*60) 16*3600+1.7*3600
16*3600+1.9*3600 18*3600 18*3600+20 18*3600+60
range(18*3600+2*60,60,18*3600+12*60)
range(18*3600+15*60,3*60,18*3600+30*60)
range(18*3600+40*60,10*60,18*3600+3.5*60*60) 18*3600+3.7*3600

21 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

18*3600+3.9*3600 22*3600 22*3600+20 22*3600+60

range(22*3600+2*60,60,22*3600+12*60)
range(22*3600+15*60,3*60,22*3600+30*60)
range(22*3600+40*60,10*60,22*3600+1.5*60*60) 22*3600+1.7*3600
22*3600+1.9*3600 24*3600.

4 Click to expand the Values of dependent variables section. Locate the Values of
Dependent Variables section. Find the Values of variables not solved for subsection.
From the Settings list, choose User controlled.
5 From the Method list, choose Solution.
6 From the Study list, choose Full Model: Load History, Time Dependent 2.
7 From the Time (s) list, choose All.

Solution 3 (sol3)
1 On the Study toolbar, click Show Default Solver.
Force solver to calculate solution in specified time steps.
2 In the Model Builder window, expand the Solution 3 (sol3) node, then click
Time-Dependent Solver 1.
3 In the Settings window for Time-Dependent Solver, click to expand the Time
stepping section.
4 Locate the Time Stepping section. From the Steps taken by solver list, choose Strict.
5 Click to expand the Output section. Clear the Store time derivatives check box.
6 On the Study toolbar, click Compute.

Simulate the fatigue load cycle on the submodel.

Evaluate fatigue response of the critical solder joint.

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Structural Mechanics>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

22 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 2 (FTG2)

Energy-Based 1
1 In the Model Builder window, under Submodel (comp2) right-click Fatigue 2 (ftg2) and
choose the domain evaluation Energy-Based.
2 Select Domains 4–6 only.
3 In the Settings window for Energy-Based, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Darveaux.
5 From the Energy type list, choose Creep dissipation density.
6 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics 2 (solid2).
7 Locate the Evaluation Settings section. From the Volume average method list, choose
Entire selection.
8 Locate the Fatigue Model Parameters section. In the a text field, type 2.6457e-4.

MATERIALS

Solder, 60Sn-40Pb (mat8)

1 In the Model Builder window, expand the Submodel (comp2)>Materials node, then
click Solder, 60Sn-40Pb (mat8).
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Crack initiation energy K1_Darveaux 13173 1 Darveaux
coefficient
Crack initiation energy k2_Darveaux -1.45 1 Darveaux
exponent
Crack propagation energy K3_Darveaux 3.92e-7 m Darveaux
coefficient
Crack propagation energy k4_Darveaux 1.12 1 Darveaux
exponent
Reference energy density Wref_Darveaux 689 J/m³ Darveaux

23 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

ADD PHYSICS
1 On the Home toolbar, click Add Physics to open the Add Physics window.
2 Go to the Add Physics window.
3 In the Add physics tree, select Recently Used>Fatigue (ftg).
4 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

5 Click Add to Component in the window toolbar.

6 On the Home toolbar, click Add Physics to close the Add Physics window.

FATIGUE 3 (FTG3)

Energy-Based 1
1 In the Model Builder window, under Submodel (comp2) right-click Fatigue 3 (ftg3) and
choose the domain evaluation Energy-Based.
2 Select Domains 4–6 only.
3 In the Settings window for Energy-Based, locate the Fatigue Model Selection section.
4 From the Criterion list, choose Darveaux.
5 From the Energy type list, choose Creep dissipation density.
6 Locate the Solution Field section. From the Physics interface list, choose Solid
Mechanics 2 (solid2).
7 Locate the Fatigue Model Parameters section. In the a text field, type 2.6457e-4.

ROOT
On the Home toolbar, click Windows and choose Add Study.

ADD STUDY
1 Go to the Add Study window.

24 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

2 Find the Physics interfaces in study subsection. In the table, enter the following
settings:

3 Find the Studies subsection. In the Select study tree, select Preset Studies>Fatigue.
4 Click Add Study in the window toolbar.
5 On the Home toolbar, click Add Study to close the Add Study window.

STUDY 4
1 In the Model Builder window, click Study 4.
2 In the Settings window for Study, type Submodel: Fatigue Evaluation in the
Label text field.

SUBMODEL: FATIGUE EVALUATION

Step 1: Fatigue
1 In the Model Builder window, under Submodel: Fatigue Evaluation click Step 1: Fatigue.
2 In the Settings window for Fatigue, locate the Values of Dependent Variables section.
3 Find the Values of variables not solved for subsection. From the Settings list, choose
User controlled.
4 From the Method list, choose Solution.
5 From the Study list, choose Submodel: Load History, Time Dependent.
6 From the Time (s) list, choose From list.
7 From the list select time steps from 6.48e4s to 8.64e4s.
8 On the Home toolbar, click Compute.

RESULTS
Compare the creep dissipation in the full model and in the submodel.

Dissipation History
1 In the Model Builder window, under Results click Dissipation History.
2 In the Settings window for 1D Plot Group, click to expand the Title section.

25 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

3 From the Title type list, choose None.

4 Locate the Plot Settings section. Select the y-axis label check box.
5 In the associated text field, type Creep dissipation density (J/m^3).
6 Click to expand the Legend section. From the Position list, choose Upper left.
7 In the Model Builder window, expand the Dissipation History node, then click Point
Graph 1.
8 In the Settings window for Point Graph, click to expand the Coloring and style
section.
9 Locate the Coloring and Style section. From the Color list, choose Red.
10 Click to expand the Legends section. Select the Show legends check box.
11 From the Legends list, choose Manual.
12 In the table, enter the following settings:

Legends
Full model

13 Right-click Results>Dissipation History>Point Graph 1 and choose Duplicate.

14 In the Settings window for Point Graph, locate the y-Axis Data section.
15 In the Expression text field, type [Link].
16 Locate the Data section. From the Data set list, choose Submodel: Load History/
Solution 3 (4) (sol3).
17 Locate the Selection section. Select the Active toggle button.
18 Select Point 9 only.
19 Locate the Coloring and Style section. Find the Line style subsection. From the Line
list, choose Dashed.
20 From the Color list, choose Green.
21 Locate the Legends section. In the table, enter the following settings:

Legends
Submodel

Dissipation History 1
1 In the Model Builder window, right-click Dissipation History and choose Duplicate.
2 In the Settings window for 1D Plot Group, type Shear Creep History in the Label
text field.

26 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

3 Locate the Plot Settings section. In the y-axis label text field, type Shear creep.
4 Click to expand the Title section. From the Title type list, choose None.

Shear Creep History

1 In the Model Builder window, expand the Results>Shear Creep History node, then
click Point Graph 1.
2 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Full Model>Solid
Mechanics>Strain (Gauss points)>Creep strain tensor, local coordinate
system>solid.ecGp13 - Creep strain tensor, local coordinate system, 13 component.
3 In the Model Builder window, under Results>Shear Creep History click Point Graph 2.
4 In the Settings window for Point Graph, click Replace Expression in the upper-right
corner of the y-axis data section. From the menu, choose Submodel>Solid Mechanics
2>Strain (Gauss points)>Creep strain tensor, local coordinate system>solid2.ecGp13 -
Creep strain tensor, local coordinate system, 13 component.
5 On the Shear Creep History toolbar, click Plot.
Verify that the temperature is the same in both studies.

Temperature History
1 In the Model Builder window, under Results click Temperature History.
2 In the Settings window for 1D Plot Group, click to expand the Title section.
3 From the Title type list, choose None.
4 Locate the Plot Settings section. Select the y-axis label check box.
5 In the associated text field, type Temperature (C).
6 In the Model Builder window, expand the Temperature History node, then click Point
Graph 1.
7 In the Settings window for Point Graph, click to expand the Coloring and style
section.
8 Locate the Coloring and Style section. From the Color list, choose Red.
9 Click to expand the Legends section. Select the Show legends check box.
10 From the Legends list, choose Manual.
11 In the table, enter the following settings:

Legends
Full model

27 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY

 Solved with COMSOL Multiphysics 5.2

12 Right-click Point Graph 1 and choose Duplicate.

13 In the Settings window for Point Graph, locate the y-Axis Data section.
14 In the Expression text field, type solid2.T.
15 Locate the Data section. From the Data set list, choose Submodel: Load History/
Solution 3 (4) (sol3).
16 Locate the Selection section. Select the Active toggle button.
17 Select Point 9 only.
18 Locate the Coloring and Style section. Find the Line style subsection. From the Line
list, choose Dashed.
19 From the Color list, choose Green.
20 Locate the y-Axis Data section. From the Unit list, choose degC.
21 Locate the Legends section. In the table, enter the following settings:

Legends
Submodel

28 | ENERGY-BASED THERMAL FATIGUE PREDICTION IN A BALL GRID ARRAY