Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion

Engineering Fracture Mechanics xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion Matteo Corbetta ⇑, Claudio Sbarufatti, Andrea Manes, Marco Giglio Politecnico di Milano, Dipartimento di Meccanica, via La Masa 1, 20156 Milan, Italy

a r t i c l e

i n f o

Article history: Received 3 November 2014 Received in revised form 25 June 2015 Accepted 6 July 2015 Available online xxxx Keywords: Fatigue crack growth Bayesian filtering Random loading Remaining useful life

a b s t r a c t Particle filters are effective tools for the monitoring of damage propagation phenomena. However, a common hypothesis of particle filters for damage prognosis is the constant-amplitude fatigue loading affecting the damage growth. This work constitutes a preliminary analysis of the performance of particle filtering in case of random load, relaxing the hypothesis of constant-amplitude fatigue. Two case studies referring to stationary random loads are introduced: the first concerning a narrow-band stress history, while the second focusing on a wide-band stress spectrum. A solution for each case study is provided and validated using numerical simulations of fatigue cracks in a metallic plate. Ó 2015 Published by Elsevier Ltd.

1. Introduction Bayesian filtering methods prove to be powerful tools for the monitoring and prediction of degradation processes. Especially in the field of structural degradation, their ability to explore and filter the state-space of the system allows a reliable and fast estimation of the remaining useful life (RUL) of the structure. The Bayesian updating of the RUL can be performed by means of traditional non-destructive techniques or real-time structural health monitoring systems that provide new information on the degradation state. Among several Bayesian filtering algorithms, particle filter (PF) (also known as sequential Monte Carlo sampling, or sequential importance sampling/resampling) is particularly suitable for degradation processes because of its capability to describe nonlinear systems affected by non-Gaussian probability density functions (PDFs) of the system variables or model parameters. It has been applied to fatigue-induced degradation processes frequently showing promising results. Among the most relevant applications of PF for structural integrity assessment, Orchard et al. applied PF algorithms to monitor and predict fatigue crack propagation of a turbine engine as well as a planetary carrier plate [1–3]. Other applications of PF for fatigue crack growth (FCG) can be found in [4–7]. All of them grounded on the well-known Paris’ law [8] for the estimation of the FCG rate defined as length per load cycle. PFs grounding on the NASGRO model [9] are available in [10,11]. Focusing on PF for the prognosis of composite materials, several works have been recently published by Chiachio et al. [12–14], to name just a few. In these papers, the authors made use of a modified Paris’ law based on the strain energy release rate to estimate the matrix crack density growing into cross-ply laminates. They predict the evolution of the matrix crack density as well as the stiffness reduction of the material, and these predictions were combined to estimate the RUL of carbon fiber-reinforced polymer coupons. In [14], they used a PF algorithm with the technique of subset simulation to enhance the prognostic capabilities. Instead, the work in [15] shows the application of PF for creep problems. All of the cited publications try to solve one or several specific problems of the real-time prognosis of ⇑ Corresponding author. Tel.: +39 02 2399 8213. E-mail address: [email protected] (M. Corbetta). http://dx.doi.org/10.1016/j.engfracmech.2015.07.008 0013-7944/Ó 2015 Published by Elsevier Ltd.

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

2

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

Nomenclature C empirical parameter of the fatigue crack growth rate model d parameter of the measurement system f r ; f r;1 ; f r;2 stress history frequencies F crack shape function k discrete time index K stress intensity factor m empirical parameter of the fatigue crack growth rate model N load cycle Ns number of particles Nf failure cycle t time u input of the system R stress ratio S stress S0 mean stress Sa stress amplitude w weight of the particles x state variable (semi-crack length) x0 initial semi-crack length xCR critical semi-crack length z observation of the system state b0 ; b1 parameters of the measurement system c empirical parameter of the fatigue crack growth rate model d; Kronecker delta g measurement noise h parameter vector m process noise r2DS stress signal variance

r2D

variance of the stress disturbance process

dx dN

fatigue crack growth rate

f ðÞ

evolution equation

gðÞ hðÞ i Nð; Þ Nð0; r2D Þ

observation equation general fatigue crack growth rate model particle index Gaussian probability density function stress disturbance process

MVNð; Þ multivariate Gaussian probability density function Uð½; Þ pðjÞ

uniform probability density function conditioned probability density function

CðÞ fS D

Gamma function equivalent stress range

fatigue-induced damage propagation. However, none of them assume a variable load condition, typical of aeronautical, civil, and some mechanical structures. The first application of PF for variable-amplitude load conditions has been recently proposed in [16]. The authors used the Huang’s model to mitigate the effect of the load ratio R and the retardation effect due to overloads.1 As indicated in the paper, the load condition is created by one load block with two different maximum load levels, and a well-defined number of cycles for each load level. The minimum load remains the same throughout the whole test, and the overstress (i.e. the value of the maximum load in the block) changes according to the experiments. Another recent tutorial on PF shows an identical application [17]. Even if the load condition is defined as a variable-amplitude fatigue, the load applied on the structure is deterministic, thereby the instantaneous value of the fatigue load in terms of mean and amplitude is known. Relevant works in the field of real-time prognosis and uncertainty in prognostics under variable-amplitude and random loading conditions has already been published by Sankararaman et al. [18,19]. These works underline the importance of a proper

1 The mitigation of the retardation effect in the Huang’s model is based on the model initially proposed by Wheeler in order to take the overload effects into account.

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

3

evaluation of uncertainties in prognostics, especially for the estimation of the remaining lifetime of degrading systems. However, an application of a PF algorithm for structures subject to non-deterministic loads is still missing in literature. Hence, this paper discusses the potential use of particle filter algorithms to monitor fatigue-induced degradation in the presence of random load conditions. Given the complexity of the topic, the present work analyzes two case studies under the hypothesis of stationary random loads. The main issues and difficulties potentially encountered during the implementation of the algorithm are analyzed. The first case study relates to a narrow-band stress history with null average stress. If the stress history satisfies the hypothesis of narrow-band and zero mean value, the stress ranges DS are distributed as a Rayleigh PDF. A sampling method grounding on the Rayleigh PDF is then inserted in the PF framework to monitor the crack growth process and to predict the RUL of the structure. The second case refers to a wide-band stress history with a mean value different from zero. The methodology based on the Rayleigh PDF cannot be applied here. Thus, the load range DS required to run the PF is estimated within the algorithm using the artificial dynamics method. The approach is grounded on a different idea, i.e., the load range becomes a parameter of the crack propagation model. Further details are provided in the related sections. The use of PF to predict the RUL of random loaded structures constitutes a novelty in the real-time prognostics panorama, and the main contribution of the work is to provide efficient strategies to model the stress range DS that randomly changes with time. Simulated FCG tests on a virtual metallic plate validate the proposed approaches (the experimental validation of the proposed methods will be a matter of future research). The performance of each method is discussed in terms of crack growth monitoring and remaining lifetime estimation. The paper is structured as follow. Section 2 briefly describes the problem of real-time damage prognosis with PF in the presence of random loading conditions, while Section 3 refers to the PF algorithm tailored to fatigue crack growth problems. Section 4 analyzes two cases of FCG under random loads. A method to effectively describe the stress range within the PF is presented and analyzed for each case study. A critical analysis of the potential use of PF in a real-time prognostic perspective is provided. Section 5 concludes the paper with general guidelines for future and further research. 2. Real-time damage prognosis using PF in the presence of random loading conditions As mentioned in Section 1, PF is a consolidated method for dealing with the prognosis of structures subjected to constant-amplitude fatigue. As a matter of fact, the prediction of the RUL requires the step-by-step simulation of each sample of the system state until a critical condition, using the same stress range and stress ratio during the entire fatigue life. When all the samples representing the degradation state reach the threshold, the algorithm stops and the PDF of the remaining lifetime is computed. The case of random load conditions is rather more complicated. If the structure is subjected to random loads, the instantaneous value of the load is unknown, and the future loads are unknown as well. Therefore, the stress range required to run the algorithm is not deterministically known, but it has to be estimated (and predicted). In this case, the real-time prognosis of FCG becomes a challenging problem. The list below reports some of the complications encountered during the implementation of the PF algorithm for random loaded structures. 2.1. Uncertainties on the current load condition In a structure subjected to random load, the load features are known in statistical terms. If a real-time usage monitoring system is available, then the well-established rainflow method (presented in a simplified form in the standards [20]) can be applied to extract the mean load and the load range. However, the rainflow load counting is an approximated extraction of the actual load affecting the damage, which remains unknown. As a result, the FCG rate in terms of length per cycle can only be approximated using the values derived from the rainflow load counting procedure. 2.2. Synchronization of the real-time PF with the real system Since the PF algorithm is supposed to operate in parallel with the real system (i.e. the structure subjected to random loads), a perfect synchronization between the acquisition system of the load and the prognostic system operation is mandatory. Thus, the algorithm is able to associate the load cycle occurrences to a specific time window, and to reproduce a virtual crack evolution with the same load cycles in the same time window. 2.3. Retardation and acceleration effects due to the load sequence This problem is well-known in the field of fatigue crack growth, and some authors report that the load sequence affects also the damage propagation in composite materials [21]. The most widespread model for the retardation effect due to overloads is the Wheeler’s model (which is applicable to metallic structures); it defines a modified crack growth rate because of the negative stress value in the plastic region around the crack tip after an overload. However, the effect of under-loads is not accounted for with this model. The load sequence effect is well-represented in case of a constant-amplitude load with one or several overloads altering the load history, but the implementation of retardation models remains difficult for random load Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

4

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

histories. In addition, the number of future overloads or under-loads is unknown; this number should be estimated from previous experimental data and consequently embedded in the prognostic unit using statistical updating approaches. 2.4. Uncertainty of future loads and operating regimes One of the main difficulties affecting the algorithms for real-time prognosis is the uncertainty affecting the future load conditions. If the load history has some stationary statistical features, the real-time monitoring of the load allows the reproduction of several load histories with the same features [22]. These simulated load histories can be used to simulate the damage propagation from the current state until the critical condition. Nonetheless, the operation regime can change in the future according to the usage of the structure or as a consequence of an unexpected change of the external loads. Even though studies on the load spectrum of the monitored system can be performed in advance, they are statistical representations of the actual load history. The matters reported above relate to the random load conditions only, and are an addition to the other critical complexities affecting the real-time prognosis of structures under constant-amplitude fatigue (which has not been solved yet). The interested reader can refer to [23] for an extensive discussion on the uncertainties affecting prognostics and the remaining life prediction. A simple example of uncertainties influencing the RUL prediction of structures is represented by the uncertainty of the FCG model parameters [24], which markedly affects the prediction capability of PF [25]. Consequently, the load condition becomes an additional source of uncertainty crucial for the efficiency and the effectiveness of real-time prognostic approaches. Given the complexity of the matter, a comprehensive and exhaustive dissertation of the topic goes beyond the scope of a single paper. Hence, this work imposes some restrictions on the possible scenarios: (i) Hypothesis of stationary random loads. Therefore, the random process describing the stress history has some properties including the mean and variance that do not change with time. The most general case of non-stationary random loads is not considered in this work. (ii) Stress histories have a finite number of relevant values in the frequency domain. The first case study adheres to the hypothesis of a narrow-band stress history with a zero-mean stress. Then, the hypothesis of a narrow-band stress history is relaxed in the second case study, in which the load history has a non-zero mean value, and two main frequencies. (iii) The stress ranges DS extracted from the simulated stress histories with the rainflow method are representative of the true stresses affecting the crack propagation. The target crack growths (simulated cracks in a thin plate) are built with the stress ranges extracted from the simulated stress history. This hypothesis is necessary to validate the results using virtual crack growth tests, as performed in this work. (iv) Proportionality between the number of load cycles and the considered time window. Considering a time window T, the number of load cycles affecting the crack growth is proportional to the length of T. Thereby, a direct proportionality between the number of load cycles (that are extracted with the rainflow method) and the time-to-failure of the system is available, and the lifetime estimation of the algorithm can be presented as the number of load cycle-to-failure, as in the case of constant-amplitude fatigue loads. Despite the several restrictive hypotheses reported above, the development of a Bayesian framework for random loaded structures remains a challenging problem. The nonlinearity describing the damage propagation prevents closed form solutions for the residual lifetime of the system. Approximations of the crack growth evolution using simple features of the load history including the mean value and the load range produce large errors in the lifetime prediction, even if the load sequence effect is not present. Hence, the introduction of a time-varying input in the PF algorithm is addressed using different techniques aimed at monitoring the evolution of the crack and correctly predicting the residual lifetime of the structure. 3. Particle filtering framework for fatigue crack propagation This section summarizes the PF algorithm providing the fundamental equations for the conditioned estimation of the state and the RUL. The model of the system evolution, the observation equation, and the theoretical solution of the conditioned posterior probability are shown. The approximation of the posterior PDF using the PF algorithm is briefly described, as well as the prediction equation used to predict the system state in the future. Then, the FCG rate model used to simulate the crack propagation, and the lognormal random process suitable for monotonic degradation phenomena are presented. Eventually, the equation to estimate the RUL is given. The explanation below abstains from deep mathematical analyses of the method. The interested reader can refer to the specific literature for a detailed treatment of the algorithm ([26,27], and the references therein). 3.1. Dynamic state-space model and posterior distribution estimation Let us suppose a time-varying system state x defined by the evolution equation f in the discrete time domain (1), and an observation equation g able to link the measured quantities z with the system state x (2). The first order Markovian assumption holds. Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

5

xk ¼ f ðxk1 ; hk1 ; uk1 ; m k1 Þ

ð1Þ

zk ¼ gðxk ; gk Þ

ð2Þ

The subscript k indicates the discrete time-step index, while the term uk1 is the input of the system. The couple of Eqs. (1) and (2) define the dynamic state-space (DSS) model. A direct measure of the actual system state is not feasible but it can be estimated assuming that the observation equation g is known. Since the observation equation is altered by one or several sources of uncertainty (a reasonable assumption concerning common measurement systems), and the real system evolution can be affected by unknown processes or events not accounted for in the model f, two artificial random processes m and g are inserted in the DSS formulation. These random processes are usually called noises altering the deterministic evolution of the system and the perfect measurement of the system state. Thus, the monitoring of the system state becomes a conditioned probability problem. It means that the probability distribution of the system state depends on the observations pðxk jz0:k Þ. Here, the subscript 0 : k indicates the series of observation from the beginning of the degradation process to the step k. Assuming a known posterior distribution at the previous time step, pðxk1 jz0:k1 Þ, the mathematical formulation of the problem is represented by the Chapman–Kolmogorov Eq. (3) (prediction) and by the Bayesian updating procedure (4) (updating).

Z

pðxk jz0:k1 Þ ¼

pðxk jxk1 ; z0:k1 Þpðxk1 jz0:k1 Þdxk1

ð3Þ

X

pðxk jz0:k Þ ¼

pðxk jz0:k1 Þpðzk jxk Þ pðzk jz0:k1 Þ

ð4Þ

pðxk jz0:k1 Þ is the prediction PDF of the system state given the observations, while X is the domain of the state variable xk1 . The updating of the prediction PDF is made according to Bayes’ rule depending on the likelihood of the observations pðzk jxk Þ, and the evidence pðzk jz0:k1 Þ ensuring that the integration of the posterior PDF is equal to one. The analytical solution of (3) and (4) is not available except for the linear-Gaussian cases, but the FCG process is nonlinear, and the PDFs involved in the process can be non-Gaussian, too. As a result, the Monte Carlo approach allows the approximation of the posterior distribution in (4) by means of a series of samples of the system state xðiÞ (5), where the superscript i defines the ith sample. These samples xðiÞ are usually called particles in the PF framework. Each particle is properly weighted according to the uncertainties affecting the DSS model, which are the uncertainties of the evolution equation and the observation equation. The term wðiÞ in (5) is the normalized weight associated to the ith particle xðiÞ . The subscript k indicates that ðiÞ

the estimation refers to the kth time-step. The normalized particle weights wk needed to estimate the posterior PDF in (5) can be calculated using (6) and (7).

pðxk jz0:k Þ 

X ðiÞ wk dx i ðiÞ

~ ðiÞ ~ ðiÞ w k ¼ wk1

ð5Þ

ðiÞ k ;xk

ðiÞ

ðiÞ

pðzk jxk Þpðxk jxk1 Þ

ð6Þ

ðiÞ ðiÞ ðxk jx0:k1 ; z0:k Þ

p

~ ðiÞ w ðiÞ wk ¼ P k ðiÞ ~ i wk

ð7Þ

As visible from (6), the formulation of the particle weights needs the transition density function pðxk jxk1 Þ linking two subsequent samples of the system state, the importance density function pðxk jx0:k1 ; z0:k Þ used to draw the samples, and the observation likelihood pðzk jxk Þ. Once the posterior PDF of the current condition is available, the RUL distribution can be estimated as described in Section 3.2. Then, Sections 3.3 and 3.4 provides some details on the PDFs and the likelihood involved in the algorithm assuming that the system state is an unidimensional variable representing the crack length, so xk ! xk 2 ½0; Rþ . 3.2. Prediction and RUL estimation Eqs. (5)–(7), describe how to compute the posterior distribution of the system conditioned to the observation sequence. The prediction of the future system evolution initiates with the posterior estimation in (5), and requires the estimation of the system state l step-ahead the current time step k (8), [26].

pðxkþl jz0:k Þ ¼

Z

pðxk jz0:k Þ

X

kþl Y

pðxj jxj1 Þdxk:kþl1

ð8Þ

j¼kþ1 ðiÞ

The prediction equation is rewritten using the approximation of the posterior distribution obtained with the samples xk , (9).

pðxkþl jz0:k Þ ¼

Z Ns kþl X Y ðiÞ ðiÞ wk pðxkþ1 jxk Þ pðxj jxj1 Þdxkþ1:kþl1 i¼1

X

ð9Þ

j¼kþ2

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

6

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

The estimation of the RUL is grounded on the prediction Eq. (9), assuming a pre-determined threshold of the system state xCR . Each sample is propagated with time using the transition density function pðxk jxk1 Þ until all of the particles reach the threshold state [26]. The time (or number of load cycles) required to reach the threshold is the RUL prediction associated to the sample. The algorithm therefore stops when all the particles reach the threshold xðiÞ ¼ xCR . In line with this approach, the evaluation of the RUL of the system has the same form of the posterior probability of the state (10). The subscript k in (10) indicates that the RUL has been calculated with the available information at the current time-step k.

pðRULk jz0:k Þ ¼

X ðiÞ wk dRUL i

ðiÞ k ;RULk

ð10Þ

3.3. Transition density function The equation in (1) defines the evolution of the system in a discrete time step k, from xk1 to xk . Since the random process

m affects the model f, the transition from step k  1 to step k becomes a stochastic transition defined through the transition density function pðxk jxk1 Þ. The transition density function is based on the FCG rate model, and the random process m required to make the evolution equation f a stochastic equation. In this work, the Walker’s law has been adopted as the FCG rate model. It consists of a modification of the well-known Paris’ law in order to take the load ratio R, by means of an additional empirical parameter usually defined as c ¼ cðRÞ2 , into account (11). The term DK is the stress intensity factor range, and its value depends on the applied stress and the crack length as well (12).

dx ¼ CðDKÞm dN

DK ¼ F DS

ð11Þ

pffiffiffiffiffiffi pffiffiffiffiffiffi px ¼ FSmax ð1  RÞc px

ð12Þ

The term DS refers to the definition of the stress range using the Walker’s law, thus considering the effect of different load ratios. The stress range is the input of the system defined as uk in Eq. (1), i.e. uk ¼ DS. The empirical parameters C and m depend on the material. The crack shape function F usually depends on the geometry and the crack length x, and Smax is the maximum stress applied within a single load cycle. Thus, the FCG rate model can be specified using (11) and (12). Once the FCG rate model is available, the hypothesis of linear damage accumulation makes the definition of the evolution  dx  equation in (1), which assumes the form in (13), feasible. The term dN represents the value of the FCG rate (11) calculated k1 with a crack length equal to xk1 .

xk ¼ xk1 þ DN

 dx  emk1 dNk1

ð13Þ

The dependence of the current state xk on the previous value xk1 is clearly visible from (13). The term DN defines the number of load cycles between two subsequent steps. In case of a constant-amplitude fatigue, the formulation is valid if the term DN is relatively small (theoretically, DN ! 1Þ. However, it must be chosen equal to one ðDN ¼ 1Þ in case of random loading conditions, since the value of the fatigue load randomly changes. So, the quantity DN is omitted for brevity henceforth. The term em is the random process used to transform the linear damage accumulation model into a stochastic model and is therefore fundamental for the definition of the transition density function, pðxk jxk1 Þ. If m is chosen as a Gaussian process, the term em becomes a lognormal process, and every sample of such a process cannot be smaller than 0. The reason for this choice is given by the crack growth phenomenon. The crack length can only increase over time (i.e. a crack cannot decrease its size). Then, the crack growth rate added to the current crack length must be larger than zero. The formulation  dx  em satisfies this condition, which is xk P xk1 . Using the evolution equation in (13), the system state xk is driven by (14). dN k1

xk ¼ xk1 þ

 dx  ln Nðl; r2 Þ dNk1

ð14Þ

It has to be noted that the term ln Nðl; r2 Þ represents a lognormal distribution with a mean l and a variance r2 . The expected value of the lognormal random process em must be equal to 1 to produce an unbiased swarm of particles [25]. Therefore, a proper tuning of the mean and the variance of the process is required. The relationship between these two parameters can be easily evaluated imposing that E½em  ¼ 1, (15).

l¼

r2

ð15Þ

2

Assuming such a model for the estimation of the system evolution, the transition density function assumes a shifted lognormal distribution form with a shift parameter xk1 , (16).

8 h   i2 9 > > dx  < lnðxk  xk1 Þ   r22 þ ln dN = k1 pffiffiffiffiffiffiffi exp  pðxk jxk1 Þ ¼ 2 > > 2r ðxk  xk1 Þr 2p : ; 1

ð16Þ

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

7

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

3.4. Importance density function and observation likelihood The definition of the weights in (6) requires also the importance density function pðxk jxk1 ; z0:k Þ, and the likelihood of the observation, pðzk jxk Þ. The importance density function is the PDF used to draw the samples of the system state xk , and in the most general case, it depends on the previous values of the particle, and the observations, as expressed in (6). The mathematical formulation of the filter can be simplified if the importance density function is chosen equal to the transition density, i.e. pðxk jxk1 ; z0:k Þ ¼ pðxk jxk1 Þ. Such an algorithm is called bootstrap particle filter, the most widespread particle filtering forðiÞ

mulation adopted in several papers. Therefore, it is possible to draw the samples of the system state xk using the transition ðiÞ

density function in (16) assuming that the value at the previous time step xk1 is known and a sample from the lognormal random process em (17). ðiÞ

ðiÞ

xk ¼ xk1 þ

ðiÞ   dx  r2 exp  þ r r dNk1 2

ð17Þ

The term r is a random number drawn from the standard normal distribution Nð0; 1Þ. According to the hypothesis of a bootstrap particle filter, the formulation of the weights in (6) can be simplified (18). ðiÞ ~ ðiÞ ~ ðiÞ w k ¼ wk1 pðzk jxk Þ

ð18Þ

As visible in (18), the likelihood of the observation is the term driving the weights of the particles. The explicit formulation of the likelihood depends on the measurement system and its uncertainty as well. For example, the authors in [4,28] use a logit model to simulate an ultrasonic inspection device (19). Following the approach in [4,28], the likelihood linked to the logit model is expressed in (20).

ln

zk xk ¼ b0 þ b1 ln þ gk d  zk d  xk

8 ! 9 < 1 ln zk  xk 2 = 1 d dzk pðzk jxk Þ ¼ pffiffiffiffiffiffiffi exp  : 2 ; zk ðd  zk Þ rg 2prg

ð19Þ

ð20Þ

zk where the value xk is equal to ln dz without considering any uncertainty of the measurement. This measurement model is k

used for all the simulations and all the case studies of this work, and its parameters are reported in Table 1. 3.5. Considerations on particle filters in the presence of random loads All the equations provided in the previous subsections assume a constant-amplitude input represented by the stress range DS. Yet, if the load history on the structure is a random process, the stress range DS changes at each discrete time step. As a consequence, the stress intensity factor range DK changes at each time step because of two random variables: the crack length xk1 and the stress range DSk1 (21).

  pffiffiffiffiffiffiffiffiffiffiffiffi m pffiffiffiffiffiffiffiffiffiffiffiffi m dx  ¼ C ðF DSk1 pxk1 Þ ¼ C FSmax;k1 ð1  Rk1 Þck1 pxk1 dNk1

ð21Þ

The term DSk1 ¼ Smax;k1 ð1  Rk1 Þck1 is defined according to the Walker’s law. Its value depends on the current maximum stress combined with a certain load ratio and then a certain value of the parameter c. It can be obtained by a real-time usage monitoring system or estimated according to the statistical features of the load spectrum. As a result, the RUL prediction using PF becomes more complex. First, the continuous monitoring of the crack growth must take into account that the load range is not constant. Otherwise, the particle swarm would be unable to track the damage growth correctly. Secondly, the estimation of the RUL distribution depends on the estimation of the expected loads, and the load sequence must be approximated using statistical methods. In the next section, the PF has been applied to random load cases proposing a solution for the case of a narrow-band stress history (case study 1), and for the case of a wide-band stress history (case study 2). The stress ranges DS is described in statistical terms, so that the monitoring of the crack and the RUL prediction take the variability of the loads into account. 4. On the applicability of particle filters for random load conditions The literature on FCG in the presence of random loading is extensive, and several methods have been developed to assess the time to failure of structures subjected to non-deterministic loads. Thus, the present section shows the potential application of Bayesian filtering algorithms. Two case studies are presented; the first is valid for narrow-band stress histories with a zero mean value, while the second allows the estimation of the RUL in the presence of wide-band stress histories. Each approach is described within this section highlighting the initial hypotheses, the potential application, and the performance. Both case studies withstand with the hypotheses stated in Section 2, and the following simplification for the real scenario: ðiÞ Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

8

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx Table 1 Parameters of the simulated measurement system (Eqs. (19) and (20)). Parameter

Value

d b0 b1

100 0.06 1.25 0

lg r2g

0.05

neglecting the load sequence effect and (ii) neglecting the measurement noise affecting the stress histories. Supposing a perfect filtering process of the input load signal, the uncertainty of the provided load measurement is not discussed. Moreover, the sample frequency for the recording of the stress history is supposed to be large enough to capture each significant variation of the stresses affecting the crack growth (a reasonable assumption considering modern load monitoring systems). The target crack growth used as a virtual test is simulated under the condition of stationary random loads. It is simulated by the following steps: (i) the extraction of the stress cycles from a pre-defined stress history using the rainflow method [20] and (ii) the simulation of the crack growth using the DS calculated at step (i). This procedure is identical in both case studies. Hence, the stress histories shown hereafter are representative of the stress histories affecting the simulated damage. Instead, the stress range used to propagate the particles in the algorithm is statistical defined according to two specific approaches. The PF algorithm adopted in this work is a sequential importance resampling with systematic resampling [27]. Nevertheless, the methodology explained below can be implemented with any kind of PF method and any kind of resampling approach. The number of particles used in all the simulation is Ns ¼ 3000. 4.1. Target crack growth The test structure is a simulated plate with a central crack with a known analytical solution of the stress intensity factor as well as the analytical relationship between the load and the stress, see Fig. 1. Since the common FCG models for central cracks are based on the semi-crack length, the variable x indicates the semi-crack length. The Walker’s law presented in Section 3 is employed as the FCG rate model. The statistical features of the random process m influencing the target crack growth are not provided for the PF algorithm, in order to reproduce a realistic environment in which the actual uncertainty concerning the crack growth is unknown. A linear damage accumulation model as the one expressed in Section 3 is used to produce the target crack growth (22).

 pffiffiffiffiffiffiffiffiffiffiffiffi m xk ¼ xk1 þ C FSmax;k1 ð1  Rk1 Þck1 pxk1 emk1

ð22Þ

Every time a new method is tested, a new stress history is simulated according to the hypothesis of the method itself, and the rainflow method is used to extract the mean and the range of the stresses from the stress history. This simple approach helps to understand the feasibility of each method. Table 2 shows the parameter values employed for the target crack growth. These parameters remain the same for all case studies. As a consequence, the target RUL changes according to the stress history acting on the simulated plate. Initial and critical semi-crack lengths are empirically selected. The crack shape function F is driven by the geometry, assuming that the stress field at the crack tip is not affected by the dimensions of the plate. The Walker’s parameter c changes according to positive or negative stress ratios, while the empirical parameters of the Walker’s law are reasonable values for aluminum alloys. The variance of the process noise m has been empirically selected and changes according to the case study (see Table 2). 4.2. Case study 1: particle filtering for narrow-band stress histories This case assumes several hypotheses about the load acting on the structure. The proposed method follows the considerations on narrow-band load histories made in [29], which is summarized below. Let us consider a stress history SðtÞ with a zero-mean described by one main frequency, and altered by a Gaussian random noise (23). The stress ranges DS can be modelled by the Rayleigh distribution (24) with the parameter r2DS equal to the variance of the stochastic process. An example of the simulated stress history is presented in Fig. 2, highlighting the signal in the time domain. The parameters of the stress history used to reproduce the graphs are reported in Table 3.

SðtÞ ¼ S0 þ Sa sinð2pf r tÞ þ Nð0; r2D Þ

ð23Þ

"

2 # DS 1 DS exp  2 2rDS 4r2DS

ð24Þ

pDS ðDSÞ ¼

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

9

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

Fig. 1. Simulated thin plate with central crack.

Table 2 Target crack growth parameters. Feature

Symbol

Unit

Value

Initial semi-crack length Critical semi-crack length Crack shape function Intercept Paris’ law parameter Slope Paris’ law parameter Walker’s parameter for R P 0 Walker’s parameter for R < 0 Load cycle increment Variance of process noise

x0 xCR F C m

(mm) (mm) (–) pffiffiffiffiffiffiffiffiffi (mm/cycle (MPa mm)m) (–) (–) (–) (cycles) (mm/cycle)2

1 25 1 2.382e12 3.2 0.68 0 1 6/ 1

c c DN

r2m

It has to be noted that the random process affecting the signal Nð0; r2D Þ is intended to reproduce a source of variability of the load conditions, and it is called stress disturbance process here. It does not reproduce the uncertainty related to the measurement system; the uncertainty is supposed to be filtered in a previous signal processing procedure, as stated at the beginning of the section. Moreover, its variance r2D is different from the variance of the stress history r2DS driving the Rayleigh PDF. The probability plot of Fig. 3 shows as the Rayleigh distribution is a proper assumption for the equivalent stress ranges DS extracted by the stress history in (23). Under the assumption of random loading conditions, the PF algorithm requires a proper technique to propagate the swarm of particles with time. Therefore, the following paragraphs show the proposed approach for an effective implementation of the algorithm. Since random samples of the stress range can be extracted by the inverse transformation of the Rayleigh Cumulative Distribution Function (CDF), the short-term prediction can be solved assuming that each stress cycle affecting the crack has a stress range belonging to the Rayleigh distribution. Using the Rayleigh CDF (25), a random sample of the stress range can easily be extracted (26) starting from a uniformly distributed random number u  Uð½0; 1Þ. 2

 DS2 2r

PðDSÞ ¼ 1  e

DS

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi DS ¼ rDS 2 ln u

ð25Þ ð26Þ

Assuming that the applied load is stationary, the last equation remains the same during the entire life of the structure, and can be used within the algorithm for the propagation of the particles at each time step. Therefore, the propagation of the particles from step k  1 to step k expressed in (17) changes according to the samples extracted using (26), and is defined in (27).

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

10

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

Fig. 2. Simulated narrow-band stress history in the time domain.

Table 3 Stress history features according to the Rayleigh distribution method. Stress history feature Mean stress Signal amplitude Signal main frequency Standard deviation of the stress disturbance process

S0 (MPa) Sa (MPa) f r (Hz)

rD

0 60 10 30

Fig. 3. Probability plot of the stress ranges DS extracted from the stress history.

ðiÞ

ðiÞ

xk ¼ xk1 þ

  

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffim  dx r2 ðiÞ ðiÞ ðDSðiÞ Þ emk1 ¼ xk1 þ C F rDS 2pxk1 ln u exp  þ rr dN 2 k1

ð27Þ

ðiÞ

where r is a random number from a standard normal distribution. So, the particle xk is drawn from the value at the previous ðiÞ xk1 ,

time step and two random numbers from a uniform and a standard normal distribution (u and rÞ, respectively. A new sample DS is extracted for each simulated load cycle, supposing that the actual stress range changes continuously during the real crack propagation. Thus, the transition density depends on the previous value of the sample, and the stress range ðiÞ

ðiÞ

extracted from the Rayleigh PDF, pðxk jxk1 ; DSðiÞ Þ. The equation to generate random samples of the crack growth evolution using the Rayleigh PDF has been derived. The next step is the estimation of the RUL using the hypothesis of a narrow-band random load with a null mean value. The RUL prediction using PF consists of the propagation of the particles in the discrete time domain until all of the particles themselves reach a pre-defined critical state or threshold. The approach used in the case of a constant-amplitude fatigue can be easily applied to the case under discussion, provided that samples of the stress ranges from the Rayleigh PDF are used to propagate the particles as in (27). However, a faster computation of the RUL is available by means of the statistical properties of the Rayleigh PDF. So, the computational issue related to the RUL estimation can be mitigated. The method used here is grounded on the work in [29] and the stochastic integration approach used in [7] for a constant-amplitude fatigue. Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

11

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

After a brief summary of the stochastic integration method from [7], the technique is extended to the case of random load conditions hereafter. Thereby the RUL can be effectively estimated under the hypothesis of a narrow-band history with a null mean value without the time consuming step-by-step simulation of the particles. Let us consider a model h describing the FCG rate according to the applied stress range DS, the model parameters h and the current crack length x. In general terms, this model can be interpreted as a first order ordinary differential equation in which the load cycle variable N substitutes the time domain (28). Assuming the hypothesis that a closed from solution of the inteR gral hðx;dx exists, the RUL in terms of number of cycles to failure can be evaluated by means of the separation of variable DS;hÞ technique (29). Eq. (30) shows the closed form solution of the RUL using the Walker’s law as FCG rate model.2 On the other hand, if an analytical solution of the FCG integration is not feasible, the numerical evaluation of the integral in (29) allows the implementation of the method [30]. Eq. (30) can be implemented in the Bayesian filtering framework substituting the term x0 ðiÞ

with the current value of the particle xk , to obtain the RUL associated to the particle itself (31).

dx ¼ hðx; DS; hÞ dN RUL ¼ Nf  0 ¼

ð28Þ Z

Nf

dN ¼

Z

0 ð1mÞ

RUL ¼

xCR

x0

dx hðx; DS; hÞ

ð29Þ

ð1mÞ

xCR 2  x0 2 m CF m DSm p 2 1  m2 ð1mÞ

ðiÞ

RUL0 ¼

ð30Þ

ð1mÞ

xCR 2  xðiÞ 2 m CF m DSm p 2 1  m2

ð31Þ

where the subscript 0 specifies the deterministic value of the RUL, based on the value of the particle xðiÞ . However, the value ðiÞ

RUL0 does not account for the artificial noise m affecting the crack growth process from the current state until the critical state, which is fundamental for the correct operation of the algorithm. As a matter of fact, the FCG rate model for the PF algorithm becomes a stochastic model because of the presence of the noise. According to the theory of random processes, the integration of the stochastic differential equation in (32) depends on the deterministic value of the solution altered by samðiÞ

ples of the noise (33). It allows the estimation of the RUL linked to the ith particle, RULk .

dx ¼ em hðx; DS; hÞ dN

ðiÞ RULk

¼

Z

dx ¼ hðxðiÞ ; DS; hÞ

ð32Þ Z

ðiÞ k;0

RUL

0

ðiÞ RUL 

 n k;0 X X r2 mj e dN  e DN ¼ exp  þ rrj 2 j¼1 j¼1

m

ð33Þ

ðiÞ

where r j is the jth sample from a standard normal distribution. The term RULk;0 is the deterministic value of the RUL assoðiÞ

ciated to the ith particle at the time step k. So, the integral of the random process em in the domain ½0; RULk;0  can be approxðiÞ

imated using a series of samples of the random process itself, emj . The number of samples n required to estimate RULk is the ratio between the deterministic 

then n ¼

ðiÞ RULk;0 Þ.

ðiÞ RULk;0

in (31) and the delta cycles DN used to perform the simulation (since DN ¼ 1 here,

The repetition of the procedure for each particle allows the estimation of the RUL distribution using (10).

The equations from (28)–(33) are valid for a constant-amplitude fatigue. The approach is therefore extended to the case of a random loading fatigue below. Following the approach in [29] and the formulation of the Rayleigh PDF in (24), the average FCG rate can be estimated by the integration of the FCG model over the domain of the stress ranges (34). Using the Walker’s formulation for h, the closed form solution of the average FCG rate is feasible, as well as the estimation of the average RUL (35) [29].

Z dx ¼ hðx; h; DS; RÞpDS ðDSÞdDS E dN

ðiÞ E½RULk 

2

 1m m ðiÞ 1 xf 2  xk 2 ¼   pffiffiffiffim  C 23=2 F rDS p C m2 þ 1 1  m2

ð34Þ

ð35Þ

This approach supposes a constant crack shape function F typical of very simple geometries and structure dimensions much larger than the crack size.

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

12

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

The term Cðm=2 þ 1Þ represents a Gamma function evaluated in m=2 þ 1. Eq. (35) shows the expected value of the RUL ðiÞ

assuming a narrow-band random load that substitutes the deterministic remaining life RULk;0 in Eq. (33). The estimation of the RUL of the ith particle is eventually modified according to (36).

ðiÞ

RULk ¼

Z

ðiÞ  k;0

E½RUL

E½RUL

em dN ¼

0

ðiÞ

Xk;0 j¼1

ðiÞ

 Xk;0  r2 exp  þ rr j 2 j¼1

E½RUL



emj ¼



ð36Þ

4.2.1. Results of the proposed method The crack growth monitoring and the residual lifetime prediction are shown in Figs. 4 and 5, respectively. The large dispersion of the observations (black circles) simulated with the model in (19) is well-filtered by the algorithm. Since the first iteration, the 95% confidence band of the estimated crack length PDF includes the target crack growth, and the error of the estimated crack length is acceptable. Even the RUL estimation shows good results. The RUL confidence band includes the target RUL for the whole operation of the algorithm, and the expected value of the RUL correctly converges to the target RUL with time. This preliminary result suggests that the approach based on the Rayleigh distribution is suitable to be embedded in a particle filtering framework. A clear advantage of the method is the computational cost; the systematic simulation of each particle to estimate the RUL is not required, and the propagation of the particle between two subsequent time-steps is computationally efficient due to the fast sample generation permitted by the Rayleigh PDF (Eq. (27)). Nonetheless, the hypothesis of a narrow-band stress history with a zero mean value is a strong limitation preventing the applicability on most of the structures subjected to random loads. 4.3. Case study 2: equivalent stress range estimation This method is grounded on a different approach. It is based on the particle filtering algorithm for the combined state-parameter estimation. The mathematical technique was initially proposed by [31] for the estimation of unknown parameters of the model driving the particles. Now it is used to estimate the input of the system, which is the stress range affecting the crack. A detailed theoretical explanation of PF for combined state-parameter estimation is available in many papers ([32–34] to name just a few). Therefore, only the fundamental equations to implement the algorithm are summarized hereafter. The interested reader can refer to the dedicated literature for a detailed explanation of the method. As highlighted in (1), the evolution of the system state depends on one or several model parameters h. If the actual value of these parameters is unknown, they can be estimated within the particle filtering framework by the augmentation of the system state vector x. A transition density function of the model parameters is used to draw the samples needed to simulate the system evolution at the next time-step, supposing that h is independent from x (37). Therefore, each sample of the system state has a sample of the model parameters, and subsequently the evolution of the particle depends on the system state and the value of the model parameters as well (38). ðiÞ

ðiÞ

hk  pðhk jhk1 Þ ðiÞ

ðiÞ

ðiÞ

ð37Þ ðiÞ

xk  pðxk jxk1 ; hk1 Þ

ð38Þ

In this way, the weights associated to the particles represent the posterior distribution of the augmented state vector y ¼ ½x; h, i.e. the posterior distribution of the system state and the model parameters pðyjzÞ. The widespread artificial dynamics method is used to build the transition density function of the model parameters pðhk jhk1 Þ; an artificial random process is ðiÞ

added to the parameter sample hk1 at each time step. For example, if n is an artificial random process described by a multivariate Gaussian distribution with a zero mean and a covariance matrix Rnk1 (39), the evolution of the parameters h follows ðiÞ

(40), where the term nk1 is a random vector extracted from the multivariate Gaussian distribution with a zero mean and a covariance matrix Rnk1 .

pðhk jhk1 Þ ¼ MVNðhk1 ; Rnk1 Þ ðiÞ

ðiÞ

ðiÞ

hk ¼ hk1 þ nk1

ð39Þ ð40Þ

According to the method, the evolution equation for the FCG problem assumes the form in (41), where the stress range fS k . A lognormal distribution has becomes a model parameter. Thus, the parameter vector becomes a scalar variable hk ¼ ln D fS indicates an been chosen to describe the stress range, since it is intrinsically defined as a positive quantity. The symbol D equivalent stress range and has a different meaning with respect to the stress range DS defined in (12). Its value allows the monitoring of the crack growth from the observations of the semi-crack length, but it has no physical meaning. It becomes a simple model parameter used to predict the evolution of the process. Thus, the value of this stress range cannot be used to deduce the actual load acting on the structure, but it can be used to monitor and predict the evolution of the crack. A Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

13

Fig. 4. Crack growth monitoring according to the Rayleigh distribution approach. Entire crack growth simulation (a) and focus on the monitored crack length from 2 mm to approximately 8 mm (b).

Fig. 5. RUL estimation according to the Rayleigh PDF approach.

fS, since the stress range is defined as a positive quantity. The systematic evolution of lognormal distribution is assumed for D f DS is therefore implemented in the algorithm using a logarithmic form (41). Eq. (42) allows the estimation of the posterior distribution of the system state and the equivalent stress range [32].

2 ðiÞ yk

¼4

ðiÞ

xk



fS ðiÞ D k1

fS ðiÞ ln D k

8  ðiÞ < xðiÞ þ emk1 dx D fS ðiÞ  DN k1 k1 dN 5¼ k1 : f ðiÞ ðiÞ ln DS k1 þ nk1

3

ð41Þ

        fS k jzk / p zk jxk ; D fS k p xk j D fS k ; zk p D fS k jzk pðyk jzk Þ ¼ p xk ; D

ð42Þ

The last term of (42) can be neglected if the artificial dynamics approach is adopted, hence, the posterior PDF of the augfS k jzk Þ. The random process n has a decreasing variance to allow the mented state vector becomes independent of the term pð D fS. A reasonable prior distribution of D fS is used to draw convergence towards the most probable values of the estimated D fS able to fit the measures. samples at the first step of the algorithm. Then, the algorithm focuses on the best values of D These samples of the stress range are used for the estimation of the crack growth between two time steps, and they are also inserted in the algorithm for the fast computation of the RUL presented for the first case study. Similarly to the equations in fS ðiÞ is used to calculate the determin(31)–(33), the calculation of the RUL for each particle is straightforward. The sample D k

ðiÞ

istic RUL at time step k; RULk;0 (43). The random samples of the process em then allow the estimation of the RUL associated to ðiÞ

ðiÞ

the ith particle RULk , following (44). When all of the RULk are computed, the distribution of the remaining lifetime is evaluated using (10). Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

14

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

Fig. 6. Stress history simulated using (45) in the time domain.

Fig. 7. Stress history simulated using (46) in the time domain.

Table 4 Stress history properties used to simulate the crack propagation. Stress history feature Mean stress Signal amplitude First signal main frequency Second signal main frequency Variance of the random disturbance

S0 Sa f1 f2

(MPa) (MPa) (Hz) (Hz)

r2S

S1 ðtÞ

S2 ðtÞ

50 50 50 – 100

50 50 10 70 100

Table 5 Properties of the equivalent stress range for the artificial dynamics method. Equivalent stress range features

Value

Starting value

g D S0 (MPa)

30

Initial variance

r2 e log DS;0 r2 e

0.05

Decreasing variance function

log DS;k

r2log DS;0 pffiffiffiffi 2k

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

15

Fig. 8. Crack growth monitoring using the equivalent stress range for the stress history S1 ðtÞ. Entire crack growth simulation (a) and focus on the monitored crack propagation (b).

Fig. 9. RUL estimation using the equivalent stress range for the stress history S1 ðtÞ.

Fig. 10. Crack growth monitoring using the equivalent stress method for the stress history S2 ðtÞ. Entire crack growth simulation (a) and focus on the monitored crack (b).

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

16

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

Fig. 11. RUL estimation using the equivalent stress range method for the stress history s2 ðtÞ.

ð1mÞ

ðiÞ RULk;0

m ðiÞ ð1 Þ

2 x 2 x ¼ mCR ðiÞm mk fS p 2 1  m CF D

k

ðiÞ RULk

¼

Z

ð43Þ

2

ðiÞ RUL 

 k;0 X dx r2   exp  þ rr j 2 fS ðiÞ ; h h xðiÞ ; D j¼1

ð44Þ

k

5. Results of the proposed method The stress histories used for this case study assume the form expressed in (45) and (46). Figs. 6 and 7 show examples of the simulated stress histories, while the parameters of the stress signals are reported in Table 4. The features of the equivfS implemented in the algorithm are reported in Table 5, where the term k indicates the kth observation or alent stress range D time step of the algorithm.

S1 ðtÞ ¼ S0 þ Sa sinð2pf r;1 tÞ þ Nð0; r2D Þ

ð45Þ

S2 ðtÞ ¼ S0 þ Sa sinð2pf r;1 tÞ þ Sa sinð2pf r;2 tÞ þ Nð0; r2D Þ

ð46Þ

The damage monitoring and the RUL prediction of the algorithm are reported below. The algorithm is able to produce acceptable estimations of the crack length as well as an acceptable estimation of the residual lifetime. Figs. 8 and 9 are related to the case of a single-frequency stress history (S1 Þ, while Figs. 10 and 11 are related to the case of a two-frequency stress history ðS2 Þ. From the analysis of the results, the method seems to be unaffected by the stress history properties (one or more than one frequency). However, further studies on the effect of the input signal frequency are required. The proposed approach produces large fluctuations of the quantity of interest due to the high initial uncertainty fS. In particular, each particle has a different sample of D fS, and the algorithm requires several on the equivalent stress range D observations to focus on the most probable values. However, the basic algorithm adopted in this work can be substantially improved using advanced techniques for the parameter estimation like the kernel smoothing approach [32,34], or automatic tuning of the parameter variance [35]. Moreover, the initial distribution of the equivalent stress range as well as the decreasing function selected to allow the convergence of the algorithm play a key role in the overall performance of the method. Reasonable values for the initial variance and the decreasing function are heuristically selected in this work. The large dispersion of the particles at the beginning of the process could be mitigated by real-time algorithms for the tuning of the variance of the artificial dynamics, as the method presented in [35]. One of the main advantages of the approach is that neither a load monitoring systems nor a hypothesis on the distribution of the stress ranges are necessary. As a matter of fact, a reasonable approximation of the stress range value acting on the structure is enough to implement the algorithm, provided that the features of the artificial dynamics approach (initial variance and decreasing function) are correctly tuned using manual or automatic methods. Despite the simplicity of the method, promising results have been pointed out. The algorithm correctly achieves the crack monitoring and the RUL prediction. Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

17

6. Conclusions The work presented in this paper proposes a Bayesian framework based on particle filtering for the prediction of the remaining lifetime of cracked structures subjected to random loads. The application of PF to random loading fatigue constitutes a novelty in the field of real-time damage prognosis, for this reason several simplifying hypotheses have been introduced in this preliminary work. Two different case studies have been presented: both of them follow the hypotheses of a stationary random load, the absence of load sequence effects, and the constant proportionality between the length of the observed time window and the number of load cycles in that time window. The feasibility of real-time prognostics for random loaded structures is discussed, with the analysis of the results obtained in the simulated crack propagations. This is neither a comprehensive discussion of the matter, nor a quantitative comparison of the adopted methods. However, the present paper describes two promising techniques that can be implemented in a prognostics framework, and highlights key issues for future research on the topic. Both methods have shown good results. The suggested solutions allow an effective and reliable estimation of the remaining life of the structure if the hypotheses of stationary random load and the absence of load sequence effect are satisfied. Despite this encouraging outcome, the simulation of the target crack evolution with the mathematical model in (22) is a limitation of the work. The validity of the Bayesian filtering methods for random load conditions has to be proven in a relevant environment. Experimental activities on real structures are therefore the next step to verify the two techniques. Real tests in a controlled environment, followed by experiments on real structures subjected to real loads will help to understand the operation of the algorithm on real data. On the other hand, also the identification of load histories necessary to produce and generalize the results is a nontrivial problem. The computational costs of the algorithms are not quantitatively evaluated, however some conclusions can be drawn. The two algorithms are computationally inexpensive and they seem very suitable for real time applications especially when the fast computation of the RUL distribution is adopted. References [1] Orchard ME, Vachtsevanos GJ. A particle-filtering based framework for real-time fault diagnosis and failure prognosis in a turbine engine. In: Mediterranean conference on control and automation 2007, July 27–29, Athens, Greece. [2] Orchard ME, Vachtsevanos GJ. A particle filtering approach for on-line failure prognosis in a planetary carrier plate. Int J Fuzzy Logic Intell Syst 2007;7(4):221–7. [3] Orchard ME, Vachtsevanos GJ. A particle-filtering approach for on-line fault diagnosis and failure prognosis. Trans Inst Measur Control 2009;313(4):221–46. [4] Cadini F, Zio E, Avram D. Monte Carlo-based filtering for fatigue crack growth estimation. Probab Engng Mech 2009;24:367–73. [5] Cadini F, Zio E, Avram D. Model-based Monte Carlo state estimation for condition-based component replacement. Reliab Engng Syst Saf 2009;94(3):752–8. [6] Zio E, Peloni G. Particle filtering prognostics estimation of the remaining useful life of nonlinear components. Reliab Engng Syst Saf 2011;96(3):403–9. [7] Corbetta M, Sbarufatti C, Manes A, Giglio M. Sequential Monte Carlo sampling for crack growth prediction providing for several uncertainties. In: Proceedings of the second European conference of the prognostics and heath management society 2014, July 8–11, Nantes, France. [8] Paris PC, Erdogan F. A critical analysis of crack propagation laws. Trans ASME – J Basic Engng 1963;85:528–34. [9] NASA JS Centre and Southwest Research Institute. NASGRO reference manual, Version 4.02; 2002. [10] Corbetta M, Sbarufatti C, Manes A, Giglio M. Stochastic definition of state-space equation for particle filtering algorithms. Chem Engng Trans 2013;33:1075–80. [11] Corbetta M, Sbarufatti C, Manes A, Giglio M. On-line updating of dynamic state-space models for Bayesian filtering through Markov chain Monte Carlo techniques. Chem Engng Trans 2013;33:133–8. [12] Chiachio J, Chiachio M, Saxena A, Rus G, Goebel K. An energy-based prognostic framework to predict fatigue damage evolution in composites. In: Annual conference of the prognostics and health management society 2013, New Orleans, US. [13] Chiachio J, Chiachio M, Sankararaman S, Saxena A, Goebel K. Condition-based prediction of time-dependent reliability in composites. Reliab Engng Syst Saf 2015;142:134–47. [14] Chiachio M, Chiachio J, Saxena A, Rus G, Goebel K. An efficient simulation framework for prognostics of asymptotic processes – a case study in composite materials. In: Proceedings of the second European conference of the prognostics and health management society, 2014, July 8–11, Nantes, France. [15] Baraldi P, Cadini F, Mangili F, Zio E. Model-based and data-driven prognostics under different available information. Probab Engng Mech 2013;32:66–79. [16] Leem SH, An D, Lim CK, Hwang W, Choi JH. Experimental validation of crack growth prognosis under variable amplitude loads. In: 53rd AIAA/ASME/ ASCE/AHS/ASC structures, structural dynamics and materials conference 2012, Honolulu, Hawaii. [17] An D, Choi JH, Kim NH. Prognostics 101: a tutorial for particle filter-based prognostics algorithm using matlab. Reliab Engng Sys Saf 2013;115:161–9. [18] Sankararaman S, Ling Y, Mahadevan S. Uncertainty quantification and model validation of fatigue crack growth prediction. Engng Fract Mech 2011;78(7):1487–504. [19] Sankararaman S, Ling Y, Shantz C, Mahadevan S. Uncertainty quantification in fatigue damage prognosis. In: Annual conference of the prognostics and health management society, 2009. [20] ASTM International. Standard practice for cycle counting in fatigue analysis. In: Standard note E1049-85 (reapproved 2011). ASTM international 2011, West Conshohocken, PA 19428-2959, US. [21] Pascoe JA, Alderliesten RC, Benedictus R. Methods for the prediction of fatigue delamination growth in composites and adhesive bonds – a critical review. Engng Fract Mech 2013;112-113:72–96. [22] Ling Y, Shantz C, Mahadevan S, Sankararamn S. Stochastic prediction of fatigue loading using real-time monitoring data. Int J Fatigue 2011;33:868–79. [23] Sankararaman S. Significance, interpretation, and quantification of uncertainty in prognostics and remaining useful life prediction. Mech Syst Sig Proc 2015;52–53:228–47. [24] Corbetta M, Sbarufatti C, Manes A, Giglio M. On dynamic state-space models for fatigue-induced structural degradation. Int J Fatigue 2014;61:202–19. [25] Corbetta M, Sbarufatti C, Manes A, Giglio M. Real-time prognosis of crack growth evolution using sequential Monte Carlo methods and statistical model parameters. IEEE Trans Reliab 2015;64(2):736–53. [26] Doucet A, Godsill S, Andrieu C. On sequential Monte Carlo sampling methods for Bayesian filtering. Stat Comput 2000;10:197–208. [27] Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans Signal Process 2002;50(2):174–88.

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008

18

M. Corbetta et al. / Engineering Fracture Mechanics xxx (2015) xxx–xxx

[28] Myotyri E, Pulkkinen U, Simola K. Application of stochastic filtering for lifetime prediction. Reliab Engng Syst Saf 2006;91:200–8. [29] Zuccarello B, Adragna NF. A novel frequency domain method for predicting fatigue crack growth under wide band random loading. Int J Fatigue 2007;29(6):1065–79. [30] Sbarufatti C, Corbetta M, Manes A, Giglio M. Sequential Monte Carlo sampling based on a committee of artificial neural networks for posterior state estimation and residual lifetime prediction. Int J Fatigue http://dx.doi.org/10.1016/j.ijfatigue.2015.05.017). [31] Gordon NJ, Salmond DJ, Smith AFM. Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEEE Proc F Radar Signal Process 1993;140(2):107–13. [32] Liu J, West M. Combined parameter and state estimation in simulation-based filtering. In: Doucet A, de Freitas N, Gordon N, editors. Sequential Monte Carlo methods in practice. New York: Springer Science+Business Media; 2001. [33] Storvik G. Particle filters for state-space models with the presence of unknown static parameters. IEEE Trans Signal Process 2002;50(2):281–9. [34] Kantas N, Doucet A, Singh SS, Maciejowski J M. An overview of sequential Monte Carlo methods for parameter estimation in general state-space models. In: 15th IFAC symposium on system identification (SYSID), Saint-Malo, France (invited paper); 2009. p. 102. [35] Daigle MJ, Goebel K. Model-based prognostics with concurrent damage progression processes. IEEE Trans Syst Man Cybern. Part A Syst 2013;43(3):535–46.

Please cite this article in press as: Corbetta M et al. Real-time prognosis of random loaded structures via Bayesian filtering: A preliminary discussion. Engng Fract Mech (2015), http://dx.doi.org/10.1016/j.engfracmech.2015.07.008