Spectral Fatigue Analysis of Multiaxially

Loaded Components
W. Hinterberger1, O. Ertl1, C. Gaier1, H. Fleischer2
1

Engineering Center Steyr GmbH & Co KG, Magna Powertrain, Steyrer Strae 32, A-4300 St. Valentin;
{walter.hinterberger|otmar.ertl|christian.gaier}@ecs.steyr.com;
2
BMW AG, Knorrstrae 147, D-80788 Munich;
[email protected];

Abstract
Fatigue life analysis of deterministically loaded components based on Finite
Element results is already state of the art in automotive engineering and other
disciplines of mechanical engineering. However, realistic loads are more randomlike or purely probabilistic. Therefore the necessity of fatigue life assessment of
stochastically loaded components becomes more and more important, which
requires reformulation of fatigue hypotheses using probabilistic theory. In this
contribution the process and methods will be described, how to obtain local
damages from multiaxial loads described by power spectral densities. A practical
example will demonstrate the applicability of these new methods, which have been
implemented into the commercial fatigue solver FEMFAT.

1.

Introduction

In automotive engineering and other disciplines of mechanical engineering (e.g.

railway industry, wind energy plants, etc.) fatigue life analysis based on Finite
Element (FE) results is already state of the art. Especially for metal parts the
accuracy of fatigue life analysis has reached a level which is satisfactory for most
applications from a technical point of view.
The quality of fatigue life analysis depends strongly on the quality of input
parameters as e.g. material, loads and FE results. Loads are usually obtained from
direct measurements, or from multi-body simulations of dynamically loaded
components. Such kinds of loads used for fatigue life analysis are always of
deterministic nature. However, realistic loads are more random-like or purely
probabilistic as e.g. varying road and driver profiles for automotive applications,
varying wind and wave profiles for offshore plants, etc. Therefore the necessity of
fatigue life assessment of stochastically loaded components becomes more and
more important, which requires reformulation of fatigue hypotheses using
probabilistic theory. The way of thinking in terms of probability-density, frequency

domain, power spectral density (PSD), etc. differs strongly from the classical way
in terms of deterministic loads in the time domain.
Stochastic loads can be described in a practical and mathematically sound way by
so-called Power Spectral Density (PSD) functions. In this contribution the process
will be described, how to obtain local damages from multiaxial load-PSDs.
Different equivalent stress hypotheses, which have been reformulated for the use in
the frequency domain, will be presented. Furthermore methods have been
developed, to consider influences on fatigue life like stress gradient, constant
stresses, size, temperature, surface roughness, surface treatments, etc.
All these methods have been implemented into the fatigue solver FEMFAT [1]
forming the new module SPECTRAL. A practical example will demonstrate the
applicability of this new module. Comparisons to results of classical fatigue
analyses in time domain show similar accuracy, however, obtained in a strongly
reduced analysis time, which is a major advantage of spectral fatigue analysis
method over the classical approach.

2.

Fundamentals of multiaxial spectral fatigue analysis

For conventional damage analysis, loads applied to the structure are specified as a
function of time [2, 3]. If the loads are of stochastic nature and assumed to be
ergodic and normally distributed, the fatigue analysis can be performed entirely in
the frequency domain. Accordingly, the loads are specified by their power density
spectra instead of their time histories.
2.1.

The equation of motion

The general position-discrete, time-continuous equation of motion can be written

as


  + 
 +   =   .


(1)

The vector  describes the displacements relative to the stationary state while
,

and  represent the mass, damping and the stiffness matrix. The sum of all the

loads applied to the mechanical system is indicated on the right-hand side. The
vectors  are the so-called unit load cases and they describe the manner in which
the scalar load   influences the individual degrees of freedom. Based on the
assumption of small displacements and linear elastic deformation, there is a linear
correlation between the mechanical stress components  and the displacements 
 = C ,

(2)

whereby C is a constant matrix.

The method of modal reduction is commonly used in order to reduce the size of the
system of equations. Using the approach


    = 

(3)



the displacements are represented as a linear combination of constant modal forms

 with time-dependent coefficients, the modal coordinates  . Because generally
a small number of modes
in comparison to the number of degrees of freedom is
adequate for a sufficiently accurate description of the solution, the system of
equations Eq. 1 can be significantly simplified:


  +
 
 + 
  =   .




(4)

Here, the transpose of is represented by the notation and the matrices


,
 and 
  designate the modal mass, damping

and stiffness matrix respectively.

The modal approach in Eq. 2 yields

 = ,

(5)

for the mechanical stresses, whereby the modal stresses are incorporated using
C .

The modally reduced equation of motion Eq. 4 can be transformed into the
frequency domain by means of the Fourier transformation.


 + &2$
 +
 ' (  =  ) .
!2$

%



(6)

The hat symbol is used for the designation of Fourier transformed functions. Given
unique solvability, the transfer functions are defined as
 + &2$
 +
'
T+ f !2$%

and Eq. 6 can be written as

( = T+  ) 



-



(7)

(8)

Insertion in Eq. 5 yields a linear dependency of the Fourier transformed

mechanical stresses on the loads:

. = (  = T+  ) .



(9)

Common FEM program packages can output the necessary modal stresses as
well as the transfer functions T+  when a modally reduced frequency response
analysis is performed.
2.2.

From load to local power spectral densities

The power spectral density (PSD) specifies the power of a signal as a function of
frequency for an infinitesimally small frequency band. In other words, it specifies
the distribution of the total power of a signal over its constituent frequencies. The
PSD of a time signal  can be calculated using the Fourier transform .
according to

/0,0  1 . .  = 1 |.|% ,

(10)
whereby 1 designates a scale factor and . is the complex conjugate of .. In order
to arrive at a continuous spectrum starting from a discrete frequency spectrum of a
finite time signal, a local averaging in the frequency range is usually required (e.g.
averaging over the spectral lines of different time windows or the formation of a
moving average over the discrete spectral lines), which is shown here by the angle
brackets. For the estimation of the PSD of a given signal, see e.g. [4]. The scale
factor 1 is selected in accordance with accepted practice in such a way that the area
under the PSD is equal to the variance of the signal:
=

Var9: = ; /0,0  <.

>

(11)

The correlation between two different time signals   and %  can be
determined statistically with the help of the cross PSD, given by

/0? ,0@  1 .  .% 

(12)
Due to the scaling condition Eq. 11 the area under the real part of the cross PSD
equals the covariance of both signals
=

Cov , %  = ; !/0? ,0@ ' <.

>

(13)

For uncorrelated signals, the corresponding cross PSD vanishes. From Eq. 9
follows


..  = T+ ) DE  TF  ,

 E

(14)

whereby TF represents the complex conjugate transpose of TF . In further

consequence, with the help of Eq. 12, this yields a correlation between the PSDs of
the loads and the PSDs of the mechanical stress components:


/G,G  = T+   /HI ,HJ  TF 

 E

(15)

Using this equation, the 6x6-PSD matrix of the 6 stress components can be
determined for every node of the structure.
2.3.

Material characteristic equivalent stress PSDs

Fatigue assessment of local multiaxial stresses requires a damage hypothesis,

where the stress tensor is transferred to an equivalent scalar quantity, which is
directly comparable to the fatigue strength for uniaxial cyclic tension/compression
loading. This means for the spectral method, that the local 6x6 stress PSD matrix
needs to be mapped to a scalar equivalent stress PSD. In the following, several
equivalent stress definitions, which are already proved in the time domain, have
been investigated and extended for spectral fatigue:

For brittle material as e.g. cast iron, it is well known that the normal stress
hypothesis has to be applied. For multiaxial loadings it is useful to
combine it with a critical plane criterion:

KL MN
(16)
The equivalent stress is equal to the normal stress amplitude in material
planes. The plane, where reaches its maximum, is the critical one. Due to
the linear dependency on the stress tensor, the equivalent PSD of the
normal stress can be directly calculated.

For ductile materials as e.g. steel, it is common to use von Mises stress
according to distortion energy hypothesis:
KL O

%
%
%

+ %%
+ PP
 %% %% PP PP 
%
%
% 
+3R%
+ RP
+ R%%

(17)

This hypothesis implies that fatigue limit for torsion is 3 lower than for
tension-compression load. This factor is unchangeable, which is a major
disadvantage of the von Mises stress. Due to the nonlinear relation the
PSD of the von Mises stress cannot be directly calculated like that for the
normal stress. Therefore the use of an equivalent PSD has been
proposed for spectral fatigue analyses [5].

As a measure of the materials ductility the ratio of the
tension/compression fatigue limit T and the shear fatigue limit RT ,

U = T /RT , can be used. According to the Germans FKM-Guideline [6],

possible values of k are:
Steel:
k = 1.73
Nodular cast iron:
k = 1.54
Cast aluminum alloys:
k = 1.33
Gray cast iron:
k = 1.18
The ratio U can be incorporated by a material characteristic equivalent
stress similar to the von Mises stress, however, related to material planes
with normal and shear stress amplitudes MN and RN :
% + U%R%
(18)
KL WMN
N
Again, the plane which maximizes KL is assumed to be critical for fatigue
failure. Since Eq. 18 is also nonlinear, an equivalent PSD is constructed
in a similar fashion as for the von Mises stress.

Unfortunately, Eq. 18 provides exact results for the equivalent stress in
the critical cutting plane only for U 2. A (generally small) error results
for U > 2 and tension/compression stresses. This is represented for
U = 3 at the left of Fig. 1, where the normalized Mohr's circle of stress
for tension can be seen, distorted to an ellipse due to the shear stress
factor. The critical cutting plane is now defined by the point on the ellipse
furthest removed from the origin. For tensional stresses this should be at
the intersection with the normal stress axis. However, because of the
heavy distortion the point wanders upwards (red arrow). This
subsequently leads to an oversized equivalent stress. For the case U = 3
the error is 6%. This can be remedied by a modification of the equivalent
stress for U > 2:
% + U%R%
KL WU % 1 U % /4MN
(19)
N
This also compresses the ellipse in the horizontal direction until, as shown
on the right of Fig. 1, it touches the unit circle and thus provides the
correct equivalent stress.

Figure 1: Distorted Mohrs circle according to equivalent stress Eq. 18 (left) and 19 (right).

2.4.

Damage calculation

For damage accumulation, a rainflow count is generally performed for the

equivalent stress history. Because the history is not available, a probability model
is required which makes an estimate of the amplitude distribution possible based
on the equivalent PSD of the equivalent stress. An analytical derivation is only
available for the narrow band spectra, whose rainflow amplitudes N obey a
Rayleigh distribution [7]:

N
N%
exp g
h
(20)
d>
2d>
Here d> is the 0-th spectral moment of the equivalent PSD /Geq ,Geq . The spectral
moments are defined as
\Rayleigh N  =
=

dM ;  M /Geq ,Geq  < .

>

(21)

For wide band spectra the Rayleigh distribution leads to very conservative damage
results. Therefore, various semi-empirical models have been proposed. Among
them, the Dirlik model has gained wide acceptance [7]:
\Dirlik N  = 

lWd>

exp g

N

h + %

lWd>
N
N%
+ P
exp g
h,
d>
2d>

N
N%
exp
g
h
m% d>
2m% d>

(22)

whereby  + % + P = 1 applies. The values of  , % , P , m and l are

determined from the moments d> , d , d% and dn [7].

Once the distribution of rainflow ranges is determined, the total damage opoN can
be straightforwardly determined by integration
=

opoN = ; \ N  N <N ,

>

(23)

where the damage function N  describes the contribution of a stress cycle N to
the total damage. The damage function can be derived from the S-N curve.
2.5.

Influence factors on fatigue life

For fatigue analysis synthetic S-N curves are constructed, which incorporate
influences such as those of the local stress gradient, constant stresses, rotating
principal stresses, temperature, or technological influences [6, 8-10]. These
influences have been adopted and reformulated to be convenient for the frequencybased approach.

3.

Example Brake Disc Cover

Figure 2: PSD of load (left) which acts via the 4 mounting holes (right) on the brake disc cover in
normal direction.

The results of a spectral fatigue analysis for a BMW brake disc cover (see Fig. 2
right) using FEMFAT SPECTRAL are compared with those of a time-based
analysis using FEMFAT ChannelMAX. A measured time series consisting of 1.2
million points has been used as input to describe the load over time. However, the
frequency-based calculation requires PSDs as input. Therefore an equivalent PSD
has been generated out of the time series data [4]. The PSD is sampled using 2000
equidistant spectral lines with 0.25 Hz distance up to a frequency of 500 Hz

(compare Fig. 2 left). For both calculations the influence of the stress gradient
has been incorporated.
The comparison of the damage results of the time-based (Fig. 3) and the frequencybased analysis (Fig. 4) shows good consistency since the most damaged areas are
congruent.

Figure 3: Results obtained by ChannelMAX. Most damaged regions are colored red.

Figure 4: Results obtained by FEMFAT SPECTRAL. Most damaged regions are colored red.

The absolute damage values from the FEMFAT SPECTRAL analysis are about
80% larger than the damage values of the FEMFAT MAX analysis. The
differences result from the different nature of loads (deterministic load in time
domain versus stochastic load in frequency domain) and the lost phase information
of PSDs. Beside the good result correlation a major advantage of FEMFAT
SPECTRAL is the computation time acceleration, in this case a factor of about 140
compared to ChannelMAX. Further the result PSD shows the mode shape which is
causal for the damage (at 175 Hz, see Fig. 5). The information which mode shape
contributes at most to the damage can be directly used for design improvement.

Figure 5: Equivalent stress PSD in critical plane of most damaged node.

4.

Conclusion and Outlook

The presented multiaxial frequency-based fatigue analysis approach which is

implemented in FEMFAT SPECTRAL is able to reproduce the results obtained by
a corresponding time-based analysis. The good accordance was achieved through
incorporation of influence factors such as the stress gradient.
For future work, it might be interesting to investigate, whether the spectral analysis
approach could be refined in order to incorporate also mean stress effects. A
promising technique for estimating the mean stress of damaging cycles was
proposed in [11].

5.

References

[1]

FEMFAT Homepage, www.femfat.com.

[2]

Gaier C., Dannbauer H., An Efficient Critical Plane Method for Ductile,
Semi-ductile and Brittle Materials, 9th Int. Fatigue Congress, Atlanta, 2006.

[3]

Gaier C., Dannbauer H., A Multiaxial Fatigue Analysis Method for Ductile,
Semi-Ductile, and Brittle Materials, Arabian Journal for Science &
Engineering, Volume 33, Number 1B, April 2008, pp. 223-235.

[4]

Welch, P.D. The Use of Fast Fourier Transform for the Estimation of Power
Spectra: A Method Based on Time Averaging Over Short, Modified
Periodograms. IEEE Transactions on Audio and Electroacoustics, AU-15,
pp. 70-73, 1967.

[5]

Pitoiset, X., Preumont, A. and Kernilis, A, Tools for a Multiaxial Fatigue

Analysis of Structures Submitted to Random Vibrations, Proceedings of the
European Conference on Spacecraft Structures, 1998.

[6]

FKM-Guideline Analytical Strength Assessment, Forschungskuratorium

Maschinenbau, VDMA Verlag, 5th Edition, Frankfurt am Main, Germany,
2003.

[7]

Nieslony, A. und Macha, E. Spectral Method in Multiaxial Random

Fatigue, Springer, Berlin, 2007.

[8]

Eichlseder W. and Unger B. (1994), Prediction of the Fatigue Life with the
Finite Element Method, SAE Paper 940245.

[9]

Unger B., Eichlseder W. and Raab G. (1996), Numerical Simulation of

Fatigue Life Is it more than a prelude to tests?, Fatigue96, Berlin.

[10]

Gaier C., Lukacs A., Hofwimmer K., Investigations on a statistical measure

of the non-proportionality of stresses, International Journal of Fatigue 26/4,
2004, pp.331-337.

[11]

Benasciutti D., Tovo R., Spectral methods for lifetime prediction under
wide-band stationary random processes, International Journal of Fatigue,
27, pp. 867-877, 2005.