Probabilistic characterization of the effect of transient stochastic loads on the fatigue-crack nucleation time

by   Stephen Guth, et al.

The rainflow counting algorithm for material fatigue is both simple to implement and extraordinarily successful for predicting material failure times. However, it neglects memory effects and time-ordering dependence, and therefore runs into difficulties dealing with highly intermittent or transient stochastic loads with heavy tailed distributions. Such loads appear frequently in a wide range of applications in ocean and mechanical engineering, such as wind turbines and offshore structures. In this work we employ the Serebrinsky-Ortiz cohesive envelope model for material fatigue to characterize the effects of load intermittency on the fatigue-crack nucleation time. We first formulate efficient numerical integration schemes, which allow for the direct characterization of the fatigue life in terms of any given load time-series. Subsequently, we consider the case of stochastic intermittent loads with given statistical characteristics. To overcome the need for expensive Monte-Carlo simulations, we formulate the fatigue life as an up-crossing problem of the coherent envelope. Assuming statistical independence for the large intermittent spikes and using probabilistic arguments we derive closed expressions for the up-crossing properties of the coherent envelope and obtain analytical approximations for the probability mass function of the failure time. The analytical expressions are derived directly in terms of the probability density function of the load, as well as the coherent envelope. We examine the accuracy of the analytical approximations and compare the predicted failure time with the standard rainflow algorithm for various loads. Finally, we use the analytical expressions to examine the robustness of the derived probability distribution for the failure time with respect to the coherent envelope geometrical properties.


Probabilistic failure mechanisms via Monte Carlo simulations of complex microstructures

A probabilistic approach to phase-field brittle and ductile fracture wit...

On Physical Layer Security over Fox's H-Function Wiretap Fading Channels

Most of the well-known fading distributions, if not all of them, could b...

A Comprehensive Performance Evaluation of a DF-Based Multi-Hop System Over α-κ-μ and α-κ-μ-Extreme Fading Channels

In this work, an integrated performance evaluation of a decode-and-forwa...

Fishnet Model with Order Statistics for Tail Probability of Failure of Nacreous Biomimetic Materials with Softening Interlaminar Links

The staggered (or imbricated) lamellar "brick-and-mortar" nanostructure ...

Exchangeable Bernoulli distributions: high dimensional simulation, estimate and testing

We explore the class of exchangeable Bernoulli distributions building on...

Correlation structure in the elasticity tensor for short fiber-reinforced composites

The present work provides a profound analytical and numerical analysis o...

1 Introduction

Modern engineering applications increasingly rely on enormously capital intensive structures, which are placed in extreme conditions and subject to extreme loads that vary significantly throughout the structure expected lifetime. Failure costs are astronomical, including forgone profits, legal penalties, tort payouts, and reputation damage [8]. Minimizing lifetime costs require safe-life engineering and a conservative assessment of failure probabilities. Unfortunately, while material fatigue is a major contributor to failure, non destructively measuring fatigue is both difficult and expensive [30].

For many classes of structures, fatigue loads have a stochastic character with transient features that cannot be captured through a statistically stationary consideration. Examples include loads in wind turbines due to control and wind gusts, transitions between chaotic and regular responses in oil risers, and slamming loads in ship motions [13, 28, 29]

. For such applications, traditional frequency domain approaches have difficulty predicting the fatigue effects of intermittent loading, as those are inherently connected to time-ordering and therefore cannot be captured by the spectral content of the load.

An alternative approach for the prediction of fatigue-crack nucleation relies on hysteretic cohesive-law models [23, 11]. The appeal of this type of model of fatigue-crack nucleation is its applicability to arbitrary geometries, loading conditions, and loading histories, including random load cycling with transient spikes. However, quantifying the failure time under stochastic loading typically requires the generation of a large number of loading time histories and the computational solution of the hysteretic cohesive-law models for each of these loading scenarios. This Monte-Carlo scheme is associated with an enormous computational cost. Alternative approaches include statistical linearization [20, 16, 3], hierarchical modeling [31, 17], structured sampling methods [19], and optimal experimental design [1, 7, 27, 10, 12]. While there is a plethora of uncertainty quantification methods that can be potentially applied for the statistical characterization of the failure time there have been very few research endeavors on this direction. Moreover, one has to pay particular attention as very few uncertainty quantification methods can take into account the effect of transient intermittent events.

In this work, we first develop an efficient time-marching scheme for the fatigue model developed by Serebrinksy and Ortiz [23, 11]

with take into account the time-ordering effects. Based on this model we then derive analytical approximations for the probability mass function (pmf) of failure time in terms of the load probability density function (pdf) and the coherent envelope. We demonstrate the developed ideas in several examples involving intermittent loads. We also compare with standard fatigue life estimation methods and discuss their performance in various loading scenarios. Finally, we examine the robustness of the derived pmf with respect to geometrical properties of the coherent envelope.

2 The Serebrinsky-Ortiz (SO) model

We consider a single material element with one dimensional loading. The applied load is a random process given by , and the corresponding displacement, depends on the element constitutive relation. Further, this constitutive relation depends on the fatigue state of the material; after some number of loading/unloading cycles, , the material stiffness will degrade and eventually the material will fail. Our aim is in the relationship between the load and the failure time .

For the characterization of fatigue-crack initiation we consider the hysteretic cohesive-law model by Serebrinsky and Ortiz [23, 11]. In this model the fatigue-crack nucleation problem is formulated as a first-upcrossing problem of the opening displacement-load curve: with a critical boundary: the cohesive envelope.

We begin with a brief description of the cohesive model and explain in what sense it can be used for the quantification of fatigue-crack nucleation. A cohesive-law, i.e. the critical boundary is expressed as a curve in the plain. When the curve meets the descending branch of the monotonic cohesive envelope the material interface loses stability and we have the nucleation of fatigue-crack. The cohesive envelope is typically described by the relation


where are constants that characterize the material.

The next step is to characterize the evolution of the opening displacement in terms of a given loading history . We will employ a simple phenomenological model [23, 11]


where and are the loading and unloading incremental stiffnesses, respectively.

For simplicity we assume that unloading always takes place towards the origin and is determined by the unloading point. This condition fully defines . By contrast, the loading stiffness, , is assumed to evolve in accordance with the kinetic relation


where the fatigue endurance length, , is a characteristic opening displacement. The initial loading stiffness is given by , where is the first peak of the load function, and is the corresponding value of the ascending part of the cohesive envelope. Assuming that a typical initial load is much smaller than we can express in terms of the coherent envelope as


If a subsequent load peak results in intersection with the ascending part of the cohesive envelope, i.e. for then a new initial loading stiffness is defined and it evolves according to equation (3). This typically happens during extreme loading events and it is a mechanism that results in strong variability on the fatigue-lifetime of the material. By simple integration one can identify the evolution of and and estimate when the curve will eventually intersect the descending part of the cohesive envelope, when we have fatigue-nucleation.

The presented model allows to take into account not only the number and amplitude of the cycles but also their specific sequence. In this sense it can be employed to study the effect of noise but also high-frequency intermittent loads acting on the structure.

2.1 Time-discretization of the SO model

We assume that the loading function is a continuous random function with positive values. If the loading function becomes negative the opening displacement vanishes and the loading stiffness remains constant until the next positive load. Therefore, we assume a positive load and we approximate with a piece-wise linear function between local maxima/minima, i.e.

where is the number of cycle, () is the negative (positive) increment of the cycle, i.e. (), and () is the corresponding duration (Figure 1). The initial time of the cycle is denoted as , the initial load and opening displacement (both local maxima) are denoted as, and , respectively, while represents the loading stiffness. In addition, denotes the local minimum of the opening displacement during the cycle, i.e. when .

Figure 1: Piece-wise linear approximation of the loading function separates the loading and unloading part of the cycle and leads to an analytical approximation.

This approximation allows us to derive a solution in the form of an iterative map. We first express the evolution of the opening displacement during the descending part of the cycle, i.e. after the increment, . In this case,


Therefore, by direct integration we have the evolution of the opening displacement during unloading


Moreover, during unloading we have, by combining equations (2) and (3) (note that ):

Integrating by parts we obtain,


or equivalently,


where is the loading stiffness right after the end of the negative (unloading) increment of the cycle. Next, we compute the opening displacement after the positive increment of the cycle. We first compute the evolution of the loading stiffness, during loading. We combine equations (2) and (3), to obtain


which is integrated to obtain the value of during the loading part of the cycle:


At the end of the cycle we will have


The initial loading stiffness is given by , where is the first peak of the load function, and is defined through the ascending part of the cohesive envelope, as the solution of the equation


If a subsequent load peak results in intersection with the ascending part of the cohesive envelope, i.e. for then a new initial loading stiffness is defined and it evolves according to equation (3).

We utilize expression (10) into eq. (2):


By integration we have,


which is the local maximum of the opening displacement history, , at the end of the cycle. Crossing of this quantity with the descending part of the cohesive envelope (1) is an indication of fatigue-crack nucleation. Equations (6), (7), (10) and (14) provide a piece-wise linear approximation of the opening displacement evolution, . A graphical summary of the described scheme is given in Figure 2.

Figure 2: Schematic of the hysteretic cohesive-law and corresponding definition of fatigue-crack nucleation. The numbers close to the cohesive envelope curve indicate the number of loading cycles. The solid lines show loading and the dashed lines represent unloading.

Based on the equation of the cohesive envelope (1), we have material failure for the minimum number, , for which we have an up-crossing of the descending part of the cohesive envelope:


A sample of the evolution of the stiffness, , is shown in Figure 3. We observe that the evolution is typically linear except of several discrete jumps. These jumps are associated with crossings of the ascending part of the coherent envelope due to intermittent loading events. In the next section we formulate an approximation scheme that takes into account the almost linear evolution of the stiffness away from extreme events and significantly accelerates the computation.

Figure 3: Sample time evolution of for the SO model. Discontinuities indicate up-crossings of the ascending branch of the coherent envelope.

2.2 Simulation of a loading time series using the probabilistic decomposition-synthesis method

Here we formulate an approximation scheme for the failure-time under arbitrary loading time series. We first decompose the loading signal into segments associated with extreme events and regular loading events. For segments of the loading time-series where there is no up-crossing of the coherent envelope, we will show that the evolution of can be linearly approximated:


where has a constant value that we will estimate later. The remaining segments of the loading time series are associated with the discontinuous jumps in Figure 3, corresponding to the intersections of the curve with the ascending part of the coherent envelope. This breakdown of an intermittent process into a linear region and an extreme region parallels the probabilistic decomposition-synthesis framework developed in [15, 14].

Let the sequence be the discretization of such that is the local maxima of , and let be a fixed threshold. This sequence may be broken into two sets:

  • – the linear (quiescent) region

  • – intermittent spikes

For points in the set , we will use the simplified update equation (16). For points in the set , we will use the full SO update step as described in the previous section, regardless of whether the curve actually crosses the coherent envelope. Finally, for technical reasons we will remove the first few peaks () from and add them to . This helps to initialize the algorithm for the case when there are otherwise few or no early spikes in the set .

2.2.1 Estimation of the slope

We can directly estimate the slope in (16) from the SO model and the input signal statistics. Combining equations (8) and (11) we have


The parameter , which has inverse dependence on the endurance length , is small unless a very large spike happens for very small , i.e., towards the end of the material lifetime. Therefore, expanding the exponential to first order produces


This first term on the right hand side of equation (18) depends on the descending increment of the input signal, but it also depends on the small quantity and therefore may be neglected in the small-increment regime. The second term is state independent, and is directly proportional to the size of the ascending increment of the input signal. When these approximations are made, we may relate the increments of to the increments of using the following expression:


Assuming known statistical characteristics for the loading process, while the condition holds, the pdf and expected mean for are given by


where is the expected value (ensemble mean) operator, and the conditional pdf is given by


where is the indicator function. We note that one can also estimate the mean value empirically through a direct simulation of a short segment of the loading history. Under this approximation and assuming no up-crossing with the coherent envelope, the evolution of the stiffness can be approximated by:


This is a valid assumption as long as the spread of values of is not systematically large. In this case, the material life will extend until the material stiffness vanishes, i.e.


Below we consider some additional cases where analytical approximations can be obtained.

Narrow-banded Gaussian Process.

Suppose that

is a narrow-banded zero-mean Gaussian random process with standard deviation

. Since the negative values do not matter for the SO model, we need to characterize only the positive peaks. The pdf for the peaks of the process is given by (see, for instance [18])

Narrow-banded Non-Gaussian Process.

Suppose instead that is narrow-banded but not Gaussian. If the process is zero-mean then we can rely again on the amplitude of the positive peaks. The pdf of the peaks can be found as,


where is the -upcrossing rate, given by the Rice formula:


For detailed derivations of these relations we refer to [18].

Full Increment Calculation.

Suppose that we relax the requirement that the distribution of peaks be an adequate substitute for the distribution of increments. Following [25]

, we obtain an expression for the joint distribution of peak

, valley , and gap


an upper bound on the distribution of


and an expression for the distribution of increments :


Equation (28) resembles a Rice-type equation, and it in turn may be simplified by assuming a Gaussian closure for the second derivatives.

2.3 The rainflow counting algorithm

Figure 4: Sample plot for the coherent envelope in eq. (1) and the Serebrinksy-Ortiz fatigue model.

The simplest method to characterize material failure properties is to apply a harmonic load with known amplitude, , and count the number of cycles until material failure, . The relationship between and is a material characteristic; Figure 4 shows the curve for the SO constitutive model. Although realistic -curves are generally found to depend on the mean stress (compression versus tension) and loading (axial or torsional), simple cycle counting is found to work well for constant amplitude loading in the absence of intermittency [22]. This frequency domain approach is easily adapted to non-harmonic signals by the Palmgren–Miner rule, given by


which allows for calculation of equivalent fatigue by breaking the load signal into individual cycles, each of whose contributions is separately determined by reference to the -curve [4]. Based on the Palmgren–Miner rule the material fails when becomes greater than .

The well known rainflow counting algorithm [5] implements this rule by breaking a given signal into the corresponding set of increments. The use of equation (31) combined with a cycle-counting rule is standard in the literature on fatigue from random loads [26, 2, 6]. Note however that frequency domain methods such as rainflow counting completely ignores the time-ordering effects of the load history.

2.4 A stochastic load model with intermittency

To demonstrate the considered models we employ a load produced by a dynamical system exhibiting intermittent instabilities. This represents typical loads found in a wide range of engineering problems:



are white noise terms,

are constants, and are nonlinear coupling terms, and is the Young modulus. The model represents a structural mode, , which interacts nonlinearly with another mode, , creating intermittent energy transfers. It is a normal form arising in a wide range of systems related to fluid-structure interaction and mechanics, such as vortex induced vibrations, buckling, turbulent modes and waves [21]. The mode represents the ‘reservoir’ of energy (e.g. the axial mode in buckling) while the mode represents the mechanical mode which is excited by an additive noise as well as through intermittent energy transfers. The nonlinear coupling terms are associated with intermittent energy transfers from the mode to the mode. The displacement defines the evolution of the stress that the material exhibits.

Figure 5: Two loading scenarios drawn from the same stochastic system and for the same material parameters (left). The cohesive envelope (blue curves) together with the peaks of the opening displacement, , are shown (right). Although, in the first loading case the peaks are larger compared with the second case, yet the material fails with a smaller number of cycles. This early failure cannot be captured by the rainflow counting algorithm which underestimates the fatigue damage.

In Figure 5 we present two random samples (first column) of the load produced by the system (32). We can clearly observe the intermittent character of the time series. On the right column we present the cohesive envelope (1) together with the local maxima of the opening displacement (red circles). It is important to note that the two loading time series are qualitatively similar. However, while in the first case the load has a much larger maximum (grater than ) and several other intense peaks, the predicted number of cycles is , close to the one predicted by the standard rainflow counting prediction, since for this number of cycles we have , i.e. we have failure based on the Palmgren–Miner rule. On the other hand, for the second realization of the load, the peaks of the extreme events are much smaller in magnitude (the larger is close to ), yet we have fatigue-crack nucleation much earlier, . We emphasize that this material failure is inherently connected with the time-history of the load and cannot be captured by the rainflow-counting algorithm which predicts a significantly larger number of cycles until material failure. Specifically, for this number of cycles the Palmgren–Miner rule gives .

3 Analytical approximation of the failure time pmf

If not for intersections between the curve and the coherent envelope, the functional would be approximated the linear relationship (see eq. (24))


where is the initial stiffness and is the expected change in over a typical cycle corresponding to the load distribution. However, the coherent envelope changes things, by adding three effects:

  1. Extremely large loads that immediately cause material failure.

  2. Up-crossings with the ascending part of the envelope, which cause jumps (discontinuities) on the evolution of


  3. Up-crossings with the descending part of the envelope, which cause failure before has reached

Each of these effects will be considered in order to derive an analytical approximation for the pmf of the failure time. This is essential as a probabilistic quantification based on Monte-Carlo methods require, in order to accurately resolve a pmf that take takes values as small as , order of samples. While there are techniques to slightly improve on these numbers (e.g., importance sampling), this is a critical barrier of Monte-Carlo ideas to estimate the long tail for rare failure events.

3.1 Setup

We will write the coherent envelope, introduced in Section 2, in the generic form


which has both a monotonic ascending and monotonic descending branch:

and two corresponding inverses: and . Additionally, for what follows we define the functions:


which are assumed to be monotonic. These functions express the stiffness induced by an up-crossing with the ascending/descending part of the envelope and will be essential for our analysis. The monotonicity requirement is satisfied by typical coherent envelopes (e.g. eq. (1)).

3.2 Load statistics

We will assume that in the absence of envelope up-crossings the material stiffness, , evolves linearly with the number of cycles and with its gradient given by the mean value (see eq. (19)), which has been estimated through one of the methods detailed in Section 2.2.1. This is a valid assumption as long as the spread of values of is not systematically large. If this is not the case one can adopt a more complex model for the case of no envelope up-crossing using e.g. the Palmgren–Miner rule.

Further, we assume known pdf of local load maxima, , as well as a known cumulative distribution:

Finally, we will assume the independent spike hypothesis: that the amplitude of local maxima are uncorrelated. This is not true in general for narrow-banded processes, but it is approximately true for ‘large enough’ maxima, which is what is needed for our analysis. This will be a key assumption in determining the probability distribution for several intermediate quantities below.

3.3 Failure time distribution

3.3.1 Damage due to terminal loads larger than

The maximum value of is given by , and is the maximum load the material can sustain. When the load on the material exceeds , it will fail no matter the fatigue history. We call these loads terminal. The failure cycle for this terminal mechanism is given by


This is the expression for the first crossing time for the threshold . Following the independent spike hypothesis, the probability of seeing an extreme load above a given threshold may be modeled by sequential Bernoulli trials. In this case the pmf of the first cycle

when we have a spike of critical magnitude follows a geometric distribution




3.3.2 Damage due to up-crossings of the ascending part of the envelope

In general, the graph of may have multiple intersections with the ascending part, leading to multiple discontinuous jumps in (Figure 6). However, as we show below, the total fatigue lifetime effect depends only on the last such intersection with the ascending part of the envelope and at what cycle this occurs.

Invariance of the damage with respect to intermediate up-crossings. To prove this property we consider two scenarios where we have the same initial . The last point where we have intersection of the ascending part of the coherent envelope is assumed to be the same for both scenarios (having slope ), occurring at the same cycle, . In the first scenario we have an additional intersection of the envelope, at slope , occurring at some cycle, , while in the second case the only intersection occurs at (Figure 6). As we observe, for the first loading scenario, we have a jump at . Using simple geometry, the number of lost cycles due to this jump is given by (see Figure 6 right)


Using the same argument we also obtain the number of lost cycles due to the second jump occurring at :


Finally, by direct computational of the total number of cycles, we observe that this is equal to the lost cycles in the second loading scenario:

Figure 6: Two loading scenarios shown in terms of the coherent envelope up-crossings and the stiffness in terms of the number of loading cycles.

The above argument can be generalized for an arbitrary number of intermediate jumps and proves that the damage due to intersections with the ascending part of the envelope depends only on the last such intersection.
Damage quantification due to up-crossings of the ascending part for random loading. We have concluded that the number of lost cycles due to up-crossing of the ascending part of the envelope does not depend on intermediate up-crossings but only on the last up-crossing of the ascending curve. To quantify this damage we define the damage quotient ,


which is meaningful only when it is positive, i.e. only when we have an up-crossing of the envelope. The damage quotient essentially measures the magnitude of each jump on the material stiffness, expressed in a number of cycles lost due to this jump, every time we have an up-crossing. For a generic loading sequence with peaks the number of lost cycles due to up-crossing with the ascending part of the envelope will be given by the maximum of this quantity:


with the condition that this maximum is a positive number, i.e. we have at least one up-crossing with the ascending part of the envelope. If no up-crossing occurs then .

The constant used in the normalization of is the ‘initial stiffness,’ and corresponds to an infinitesimal first load peak. Its value is given by (eq. (4)),


that is, the slope of the coherent envelope near the origin.

To quantify the pmf for , we will begin by using the pdf of the load peaks, , and the monotonicty of function to obtain the pdf for :


Based on this pdf we now obtain the pdf for the damage quotient :


and the corresponding cdf:


For any sequence of independent random variables

with corresponding pdf , and cdf

, we can show, using derived distributions, that the cumulative distribution function of the maximum is given by


while using a total probability argument we can prove that the pmf for the argument of the maximum is given by [9]


Applying this result for our problem leads to,


where is a function that is independent of (i.e. it has to be computed once):


Note that


To compute the function we consider its logarithm


We multiply both sides of the equation with and consider the limit of as for typical applications has a very small value which corresponds to a large number of cycles until failure. This allows us to express the right hand side in terms of an integral.


Therefore, for finite (and small) we have


Substituting the above into (50) results in the pmf of the cycles until the last up-crossing of the ascending part of the envelope:




3.3.3 Damage due to up-crossings of the descending part of the envelope

Any intersection with the descending part of the coherent envelope will immediately cause material failure. As such, there can only be one such intersection. In order to quantify the statistics of this event we define the anticipation function :


where , is the material stiffness coefficient before cycle , and . The material stiffness can be expressed in terms of the damage quotient (eq. (42)) and its maximum value as follows:


where the pdf for is given by eq. (52). The material failure time is given by the first zero up-crossing of the anticipation function:


We first compute the pdf for . This will be given by:


Therefore, conditioning on , which expresses the maximum lost stiffness of the material due to an up-crossing with the ascending part of the envelope, and , the cycle when this up-crossing occurs, we will have


To this end, the probability of having a material failure at cycles is


where, follows the cdf (eq. (52) and (55)), while follows the pmf in eq. (50). We consider the logarithm of this equation:

For the last term we set


We multiply and divide the right hand side with and express it as an integral


where the approximation is based on the assumption of small and large number of cycles until failure. Going back to the pmf for we have


The special case where no up-crossings with the ascending part of the envelope occur is when the magnitude of the maximum stiffness jump is zero, and . In this case we will have


Using eq. (50) and assuming that probable is much larger than the probable values of (so we do not have to formally condition on ), as well as a small , we have

Finally, we integrate over the variable after we multiply with the corresponding pdf :

This expression can be also written in a more compact form




Note that the integrands in equations (69) and (70) have compact support, subsets of the interval . Expression (69) together with the functions (eq. (66)), (eq. (55)), (eq. (57)) and consist of a full approximation of the cycles until material failure. These functions are given in terms of the coherent envelope shape and load peak statistics.

3.3.4 Combined failure time

From equations (37) and (69), we have expressions for the distribution of and