Theory reference

This page presents the background theory used in Fatlab. It is not meant as a detailed explanation, just an overview of the implementation details.

1. Introduction

Fig. 1.1 shows the stress related nomenclature used. We use ranges as far as possible, except when considering mean stress corrections.


Fig. 1.1: Stress related nomenclature.

The resistance related nomenclature is shown in Fig. 1.2. We use subscript R to denote resistance. Two SN curves are shown; the dashed one is supplied by the user (\Delta\sigma_{R1},N_{R1}) and the solid line is the one used in calculations (\Delta\sigma_{R2},N_{R2}), i.e. scaled down by the partial safety factor \gamma_{Mf}. The slopes of the two curves are identical, i.e. m_1 before the knee point and m_2 after. It is possible (but not necessary) to define also two cut-off levels for the SN curve; such that stress ranges below \Delta\sigma_{Rmin} will be ignored, and stress ranges above \Delta\sigma_{Rmax} will produce a warning and set D\geq1.0.


Fig. 1.2: Resistance related nomenclature.

2. Stress-time history generation

The stress-time history, i.e. the time varying stress tensor \boldsymbol{\sigma} (t), is generated from superposition of scaled FE stresses (unit load cases) as exemplified below in Fig. 2.1.

Fig. 2.1: Generation of stress time series.

Here, a number of unit LC stresses \boldsymbol{\sigma}_i (L_{FE,i}) (nodal stresses at unit load i) are scaled with the associated load-time history L_i(t) and combined using superposition.

\displaystyle\boldsymbol{\sigma} (t) = \sum \boldsymbol{\sigma}_i (L_{FE,i}) \cdot \frac{L_i(t)}{L_{FE,i}}

The time dependent stress tensor can then be established.

\boldsymbol{\sigma}(t) = \begin{bmatrix} \sigma_{xx}(t)& \tau_{xy}(t) & \tau_{xz}(t) \\ \tau_{xy}(t) & \sigma_{yy}(t)& \tau_{yz}(t) \\ \tau_{xz}(t) & \tau_{yz}(t) & \sigma_{zz}(t) \\ \end{bmatrix}

The process is repeated individually for each node in the model.

2.1 Non-linear models

Non-linear behavior in the FE model (e.g. contact, large displacement, varying load directions, etc.) is supported by allowing two or more separate unit LC stresses for each load component and using interpolation. Four load-stress relationship types are supported as shown in Fig. 2.2.


Fig. 2.2: Load-stress relationships.

2.1.2 Bilinear

The Bilinear setting is useful e.g. when modelling bearings or pinned connections, where the path of the internal forces in the structure is significantly different, depending on the direction of the load, but not the magnitude.

\boldsymbol{\sigma} = \left\{ \begin{matrix} \boldsymbol{\sigma}_i^+ & for & L_i \geq 0\\ \boldsymbol{\sigma}_i^- & for & L_i < 0 \\ \end{matrix} \right.


Fig. 2.3: Stress response in non-linear FE model.

1D interpolation

1D Interpolation is useful, e.g. when analyzing a model which experiences large deflections. In this case, the stresses will not depend linearly on the load, because the geometry of the structure changes during the analysis.

A an arbitrary number of FE unit load cases can be supplied and Fatlab will interpolate the stress components linearly between them on an individual basis, i.e. first \sigma_x, then \sigma_y and so forth.

2D interpolation

As the name suggests, 2D Interpolation is used when the stress-load relationship needs to be interpolated along two dimensions, e.g. in the case where the load changes direction and magnitude and both have a significant influence on the stress.

This setting is implemented using Matlab’s interp2, which does not support extrapolation beyond the supplied sample points. Hence, the loading must be constrained within the minimum and maximum of the supplied FE unit load cases.

The four di erent stress load relationships can be mixed and combined as required by the FE model at hand, i.e. one load component can be linear and another load can e.g. be bilinear or interpolated.

3. Fatigue stress calculation

For fatigue calculations, in particular the cycle counting process, only one stress component must be selected or a fatigue effective stress must be calculated. Selecting a stress component, e.g.  is straight forward, however this is only meaningful in certain cases, where the global coordinate system is appropriately oriented in relation to the fatigue critical location.

In other cases, a fatigue effective stress can be used, such as the (signed) von Mises, principal stress or maximum shear stress, which will be explained in the following.

In case of significant multiaxiality, Fatlab offers the possibility to perform critical plane analysis. However, this is a computationally expensive process and should thus only be applied if necessary.

3.1 (Signed) von Mises

Two versions of the von Mises equivalent stress are available; the standard strictly positive and the so-called signed von Mises, which can take negative values.

The first is generally a poor choice for fatigue assessment because the range calculated from the strictly positive history of the von Mises stress does not include any potentially negative part of the stress cycle.

\sigma_{vm} = \sqrt{1/2( (\sigma_x-\sigma_y)^2 + (\sigma_y-\sigma_z)^2 + (\sigma_z-\sigma_x)^2 +6(\tau_{xy}^2 + \tau_{yz}^2 + \tau_{xz}^2)) }

The signed von Mises \sigma_{svm} is thus used to mediate this shortcoming, i.e. by applying the sign of the first stress invariant.

\sigma_{svm} = sign(I_1) \cdot \sigma_{vm} where I_1 = \sigma_{x}+\sigma_{y}+\sigma_{z}.

3.2 Principal stresses

The calculation of principal stresses in 3D can be a relatively cumbersome process [1]; however they can also be determined from the eigenvalues of the stress tensor.

\boldsymbol{\sigma}_{p} = eig(\boldsymbol{\sigma}(t)).

The principal stresses are defined by their algebraic magnitude, i.e. \sigma_1 > \sigma_2 > \sigma_3.

Due to the way the principal stresses are defined, none of them are suitable for fatigue assessment alone. Therefore, the numerically largest principal stress is used

\sigma_{pnmax} = \left\{ \begin{matrix} \sigma_1 & for & |\sigma_1| \geq |\sigma_3|\\ \sigma_3 & for & |\sigma_3| > |\sigma_1| \\ \end{matrix} \right.

This is probably the most commonly used fatigue effective stress and is therefore used as default. It is not fool-proof for multiaxial loading though, as discussed here.

4. Critical plane analysis

The idea behind the critical plane approach is to find the plane in which the component will crack = the critical plane. In order to do so, a number of search planes intersecting the surface either orthogonally and/or at some inclination are searched for the maximum value of a damage parameter. The plane that maximizes the damage parameter is called the critical plane. Each search plane is defined by its unit normal vector \bold{n}_s, which is again defined by the angle to the local x-axis \theta and the inclination angle \phi,  as shown in Fig. 4.1 (center).

Fig. 4.1: Critical plane analysis.

In Fatlab, the user specifies the number of search planes to use in the analysis, but beware that the execution time for the analysis grows linearly with the number of search planes. The planes are distributed evenly between 0 and 180 around the surface normal, Fig. 4.1 (left).

The stress determination, cycle counting and damage accumulation is then carried out individually for each search plane, and the plane achieving the maximum damage is selected as the critical plane.

Establishing a rotation matrix describing the direction cosines for the local frame defining each search plane facilitates the calculations

\boldsymbol{A} = \begin{bmatrix} \boldsymbol{n} & \boldsymbol{u} & \boldsymbol{v} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix}

The normal stress \sigma_n on the search plane can then be determined from the stress tensor and the rotation matrix [11]

\sigma_n = \sigma_x a_{11}^2 + \sigma_y a_{12}^2 + \sigma_z a_{13}^2 + 2(\tau_{xy} a_{11}a_{12} + \tau_{xz} a_{11}a_{13} + \tau_{yz} a_{13}a_{12})

Resolving the amplitude or range of the normal stress is trivial using traditional cycle counting techniques both for constant and variable amplitude loading.

Similarly, the two shear stress components acting in the search plane is calculated as

\tau_u = \sigma_x a_{11}a_{21} + \sigma_y a_{12}a_{22} + \sigma_z a_{13}a_{23} + \tau_{xy}(a_{11}a_{22} + a_{12}a_{21}) + \tau_{yz}(a_{12}a_{23} + a_{13}a_{22}) + \tau_{zx}(a_{13}a_{21} + a_{11}a_{23})

\tau_v = \sigma_x a_{11}a_{31} + \sigma_y a_{12}a_{32} + \sigma_z a_{13}a_{33} + \tau_{xy}(a_{11}a_{32} + a_{12}a_{31}) + \tau_{yz}(a_{12}a_{33} + a_{13}a_{31}) + \tau_{zx}(a_{13}a_{31} + a_{11}a_{33})

In the most general case, i.e. non-proportional (out-of-phase) loading, the shear stress vector will describe some trajectory on the search plane, Fig. 4.1 (right). Resolving the shear stress range in this case is far from trivial and a multitude of techniques have been proposed for this purpose. Typically some geometric element is fitted to the trajectory, e.g. the longest chord, minimum circumscribed circle, maximum rectangular hull, etc., and the shear stress amplitude and mean is determined from here. Fatlab currently only supports the Longest Chord (LC) method.

A large number of damage parameters have been proposed in the literature, however only three is currently available in Fatlab;

  1. Normal stress criterion
  2. IIW equivalent
  3. (Modified shear stress criterion, under development)
  4. Findleys criterion

Normal stress criterion

The normal stress criterion is

\Delta\sigma_{CP-N} = \smash{\displaystyle\max_{0 \leq \theta \leq \pi}} \Delta\sigma_n

In case of uniaxial/proportional loading, the normal stress only criterion corresponds to the principal stress approach. In the case of non-proportional (out-of-phase) loading, however, the use of this criterion avoids the accumulation of damage from stresses in different directions.

Modified shear stress criterion

Under development.

IIW equivalent

The IIW recommends using the Gough-Pollard criterion for welded joints subjected to multiaxial loading:

\displaystyle\left(\frac{\Delta\sigma_x}{\Delta\sigma_R}\right)^2 +\left(\frac{\Delta\tau_{xy}}{\Delta\tau_R}\right)^2 \leq k_2

Where k_2 (denoted CV in [10]) is a comparison value that takes the value 1.0 for proportional loading and 0.5 for non-proportional. This criterion is rewritten to express a fatigue effective stress range

\displaystyle\Delta\sigma_{IIW} = \frac{1}{\sqrt{k_2}}\sqrt{\Delta\sigma_x^2+k_1\cdot\Delta\tau_{xy}^2}

where k_1=\Delta\sigma_R^2/\Delta\tau_R^2. The criterion is implemented in Fatlab using a critical plane approach, i.e. maximizing the above effective stress range around the surface normal. Currently, the criterion is only available for single cycle loads.

Findley criterion

The Findley criterion is a shear stress based critical plane criterion which predicts failure on the plane that maximizes the damage parameter

f = \tau + k\sigma_{max}

Here \tau is the shear stress amplitude, \sigma_{max} is the maximum normal stress occurring over a load cycle. k_1 is an experimentally determined material factor describing the sensitivity to normal stress.

In its original formulation, infinite life was predicted if the damage parameter is below some experimentally determined threshold value f \leq f_{crit}. The Findley criterion may be reformulated [2] in terms of an equivalent stress amplitude to be evaluated against the uniaxial fatigue resistance instead of f_{crit}, i.e.

\displaystyle\Delta\sigma_F = \frac{\Delta\tau + 2k_1\sigma_{max}}{\frac{1}{2}(k_1 + \sqrt{1+k_1^2})} \leq \Delta\sigma_R

This is the formulation applied in Fatlab.

5. Cycle counting

Fatlab provides the following options for cycle counting:

  1. Reservoir counting
  2. Rainflow counting (full and half cycles)
  3. Rainflow counting (full cycles only)
  4. Single cycle
  5. Single cycle LC (Longest Chord)
  6. Single cycle MCC (Minimum Circumscribed Circle, under development)
  7. Single cycle MRH (Maximum Rectangular Hull, under development)

The reservoir counting algorithm is set as default and is recommended for short stress-time histories, because it will only produce full cycles, i.e. no half-cycles as is the case for the rainflow counting algorithm. The rainflow counting algorithm is considered the best choice for long stress-time series. In such cases, most half-cycles will be paired up, and the residual amount of half-cycles relative to the total number of cycles is insignificant.
Before any cycle counting, the stress-time history of the fatigue stress is searched for extremes, i.e. turning points. This list of extreme values is then fed through a racetrack filter order to eliminate stress ranges smaller than some threshold value. The threshold is given as a fraction of the maximum range and defaults to 0.05.
The reduced stress-time history (list of extreme values/turning points) is then used for cycle counting.
No binning of the stress ranges is done, i.e. all stress ranges are evaluated individually.

5.1 Reservoir counting

In order to use this method, the stress-time history is rearranged, such that it begins and ends with the highest peak [5]. This is accomplished by moving the part of the stress-time history prior to the highest peak to the end of the stress history (i.e. moving A in Fig. 5.1).

The stress-time history is then imagined to be filled with water and subsequently drained, one valley at the time, starting from the lowest one. Each draining procedure results in one full stress cycle and the stress-time history given in Fig 5.1 will thus be counted as 7 cycles. The mean stress value is also recorded for each cycle.


Fig. 5.1: Reservoir counting technique, after [5].

5.2 Rainflow counting

Fatlab includes two rainflow counting algorithms, denoted Rainflow half and Rainflow full. The first is developed by Nieslony [3] and available from the Matlab File-exchange. The latter is implemented after Amzallag et al. [9].

5.2.1 Rainflow half

The rainflow counting algorithm by Nieslony [3] searches the stress-time history for half-cycles and pairs up opposing sets of these. In cases where half a stress cycle cannot be paired with an opposite half cycle, residual half cycles will occur as shown in Fig. 5.2. The handling of such residual half-cycles can be problematic; however in Fatlab they are treated as any other cycles just scaled by 0.5.


Fig. 5.2: Rainflow counting technique from [4].

For long stress-time histories, only a small amount of residual half cycles will remain. However, typically the largest peak and deepest valley will remain as residual half cycles, as is also shown in the above figure.

This algorithm is recommended for timeseries, where the residual half cycles will not be closed by a subsequent repeat of the timeseries, i.e. where the timeseries actually describes the full lifetime loading of the component (or at least will not be repeated).

5.2.2 Rainflow full

If the time series at hand is expected to repeat during the lifetime of the component under analysis, all half cycles will be closed eventually. In this case the Rainflow full algorithm is recommended, because it does not return any residual half cycles, i.e. they will all be closed.

The algorithm works by sequentially extracting cycles that are smaller than or equal to its neighbors, as illustrated in Fig. 5.3. After processing the entire timeseries this way, the residual cycles is duplicated and the process is repeated. This yields an extraction of full cycles corresponding to the residual and a new residual equal to the previous.


Fig. 5.3: Extraction of cycle in Rainflow full, after Amzallag [9].

5.3 Single cycle

If the load-time history defined by the user only contains one cycle, the cycle counting procedure can be bypassed using this option. Note that load histories covering less than one cycles is not supported.

5.4 Comparison

The table in Fig. 5.4 shows a comparison of the cycle counting techniques in order to illustrate the differences. As seen, Rainflow counting produces poor results for very short stress histories, where reservoir counting performs better. For long random stress histories, however, the results are much closer and for this case Rainflow counting is commonly accepted as the best technique.

Fig. 5.4: Comparison of cycle counting algorithms for different timeseries.

5.5 Single cycle LC (longest chord)

For shear stress based multiaxial criteria, e.g. Findley, the shear stress range on the search plane must be determined. This is achieved using the longest chord method [12]. The objective is to find the longest distance between two points on the curve described by the tip of the shear stress vector on the search plane, Fig. 4.1 (right), defined as:

\Delta\tau = \max_{t_i} \{ \max_{t_j} \lvert \tau(t_i) - \tau(t_j) \rvert \}

Currently, shear stress based multiaxial criteria can only be used with single cycle loads in Fatlab.

5.6 Single cycle MCC

Under development.

5.7 Single cycle MRH

Under development.

7. Damage accumulation

The linear damage accumulation rule of Palmgren-Miner is implemented in Fatlab. Here the load spectrum is divided in to blocks of n_i load cycles of a given stress range \Delta\sigma_i. The endurable number of cycles at this stress range is denoted N_i. The partial damage resulting from block is then taken as .  The total damage is the calculated by summation of all partial damage contributions.

\displaystyle D = \sum_{i}^{q} \frac{n_i}{N_i} = \frac{n_1}{N_1} + \frac{n_2}{N_2} +...+\frac{n_q}{N_q} \leq D_{al}

Failure is usually assumed when reaching a damage sum of D=1.0, however it is well known that the real damage sum can vary greatly, e.g. in the range of 0.1-10. Some codes therefore recommend using only 0.5 or 0.2.


Fig. 7.1: Principle of Palmgren-Miner damage accumulation

In order to take into account the bi-linear nature of SN curves, the calculation of the endurable number of cycles is calculated as

N_i = \left\{ \begin{matrix} \displaystyle \ \displaystyle \frac{C_1}{\Delta\sigma_i^{m_1}} & for & \Delta\sigma_i > \Delta\sigma_{R2@\sigma_m} \\ \ \displaystyle \frac{C_2}{\Delta\sigma_i^{m_2}} & for & \Delta\sigma_i \leq \Delta\sigma_{R2@\sigma_m} \\ \ \infty & for & \Delta\sigma_i < \Delta\sigma_{min} \\ \end{matrix} \right.

Where \Delta\sigma_{R2@\sigma_m} is the stress range at the knee point corrected for mean stress. The fatigue capacity C is also calculated for both parts of the SN curve (before/after knee).

C_1 = N_{R2} \cdot \Delta\sigma_{R2@\sigma_m}^{m_1}

C_2 = N_{R2} \cdot \Delta\sigma_{R2@\sigma_m}^{m_2}

Several investigations have shown the Palmgren-Miner linear damage accumulation principle to be unreliable, particularly under fluctuating mean stresses or when pronounced sequence effects occur, e.g. repeated compressive overloads.

Another issue that is overlooked by the linear damage accumulation is that stress ranges occurring towards the end of the service life is more damaging then in the beginning. Many more advanced damage accumulation principles have been proposed to remedy these shortcomings, however they typically rely on experimental parameters that are not available in an engineering context. Palmgren-Miner is therefore still the most widely used and recognized approach.

8. Result calculation

In addition to the damage D which is considered the primary result, several other values are calculated as described in the following.

8.1 Fatigue life

The endurable fatigue life is calculated from the damage and the applied load cycles as

\displaystyle N_{endurable}=\frac{N_{applied}}{D}

The calculated endurable life assumes the same distribution of load cycles as in the applied loads.

8.2 Damage equivalent stress range

The damage equivalent stress range, considering a bi-linear SN curve, is calculated as follows [7]. The considered stress ranges are corrected to zero mean (R=-1) stress using the currently selected mean stress correction algorithm.

\displaystyle \Delta\sigma_{eq} = \left(\frac{1}{D_{al}}\frac{\sum{\Delta\sigma_{i@\sigma_m=0}^{m_1}n_i}+ \Delta\sigma_{R2}^{m_1-m_2} \cdot \sum{\Delta\sigma_{j@\sigma_m=0}^{m_2}n_j}}{\sum{n_i}+\sum{n_j}}\right)^{\frac{1}{m_1}}

Where subscripts i and j refer to stress cycles before and after the knee point, respectively. The allowable damage sum is D_{al}=1.0 by default, but can be changed under SN curve setup.

The equivalent stress range is calculated at the following numbers of cycles:

  • At NR1, i.e. the SN curve definition point.
  • At NR2, i.e. the SN curve knee point.
  • At total number of applied cycles.

8.3 Utilization ratio

The utilization ratio (on stresses, rather than life) is calculated from

\displaystyle UR = \frac{\Delta\sigma_{eq@N_{R2},\sigma_m=0}}{\Delta\sigma_{R2}}

I.e. using the equivalent stress at the same number of cycles as the knee point of the SN curve.

9. References

[1] Ugural A.C. & Fenster S.K., Advanced Strength and Applied Elasticity, Pearson Education, 2003.

[2] Bruun Ø.A. & Härkegård G., A comparative study of design code criteria for prediction of the fatigue limit under in-phase and out-of-phase tension–torsion cycles, Int. Journal of Fatigue (73), pp. 1-16, 2015.

[3] Nieslony A., Rainflow counting algorithm,, accessed 22-05-2015.

[4] Nieslony A., Determination of fragments of multiaxial service loading strongly influencing the fatigue of machine components, Mechanical Systems and Signal Processing 23 (2009) 2712–2721.

[5] Maddox S.J., Fatigue strength of welded structures, 2nd ed., Woodhead Publishing, Cambridge, UK, 1991.

[6] Schijve J., Fatigue of structures and materials, 2nd ed., Springer, 2009.

[7] McDonald K. (editor), Fracture and fatigue of welded joints and structures, Woodhead Publishing, 2011.

[8] GL, Guideline for the Certification of Wind Turbines, 2010.

[9] Amzallag C., Gerey J.P., Robert J.L. and Bahuaudl  J., Standardization of the rainflow counting method for fatigue analysis, Int. Journal of Fatigue (16), pp. 287-293, 1993.

[10] Hobbacher A., IIW recommendations for fatigue design of welded joints and
components, IIW-doc. XIII-2460-13, 2013.

[11] Socie D. and Marquis G., Multiaxial Fatigue, SAE, 2000.

[12] Lemaitre J. and Chaboche J.-L., Mechanics of solid materials, Cambridge
University Press, 1994.

2 thoughts on “Theory reference

  1. Russle

    Do you have a typo in the 2nd equation in section “3.2 Principal stresses”.

    It says
    Smax = S1 for |S1| >= |S3|, or S3 for |S3| |S1|.

    Great site, thanks.


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.