A data-driven stochastic approach to model and analyze test data on fatigue response

Computers and Structures 76 (2000) 517±531

www.elsevier.com/locate/compstruc

A data-driven stochastic approach to model and analyze test data on fatigue response R. Ganesan* Department of Mechanical Engineering, Concordia University, Montreal, Que., Canada Received 26 February 1998; accepted 28 May 1999

Abstract A stochastic approach to model and analyze test data on the fatigue response of materials and laminated composites is developed. The developed approach is `data-driven' in nature. It has been customary to describe the fatigue response of metallic and laminated composite materials using a suitable parameter that can serve as the indicator and descriptor of damage accumulation. In the present methodology, the fatigue response of the material is quanti®ed by interpreting the corresponding material parameter to be an embedded Markov process. The true probability distributions of the fatigue response parameter are extracted from sample test data based on an analytical approach, and they are used in the formulation. To this end, the Maximum Entropy Method is incorporated into the formulation. A recursive stochastic matrix equation is developed based on the test data using the theory of reliability and Fokker±Plank±Kolmogorov equation. Application of the methodology to a composite laminate is demonstrated. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Fatigue behavior; Failure of materials; Laminated composites; Reliability; Stochastic process modeling; Markov chains

1. Introduction The fatigue behavior of metallic and laminated composite materials has been studied using micro-mechanics as well as macro-mechanics approaches. The continuum damage theories correspond to the micromechanics approach. In the case of macro-mechanics approach, the response of the material to fatigue loading has been measured and quanti®ed through a suitable material parameter that is physically observable and easily measurable through laboratory experiments. The same parameter is also involved in: (i) the analysis of fatigue failure in structures, (ii) the fatigue design of

* Corresponding author. Tel.: +1-514-848-3164; fax: +1514-848-3175. E-mail address: [email protected] (R. Ganesan).

structures, and (iii) safety and reliability studies. As an example, in the case of metallic materials, the fatigue response is described and quanti®ed through the reduction in the strength properties of the material (S ) as a function of load cycles (N ). The corresponding experiments are conducted to obtain the S±N curves which will in turn be used in the fatigue design, and the safety and reliability studies. The case of composite laminates has been treated in a similar manner. The fatigue response of laminated composites is a complex mechanical phenomenon [14]. Unlike the case of metals, in the case of composite laminates various damage mechanisms and failure modes including de-lamination, ®bre breaking, debonding between ®bre and matrix, matrix cracking, ®ber pullout etc., become active at di€erent load cycles. These modes as well as their relative occurrences and contributions to the laminate fatigue re-

0045-7949/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 9 ) 0 0 1 2 4 - 8

518

R. Ganesan / Computers and Structures 76 (2000) 517±531

sponse cannot precisely be quanti®ed. Other complexities that cannot precisely be quanti®ed, also arise: failure mode changing from one to another, crack jumping, simultaneous activation of many failure modes and viscoelastic response. In this circumstance, it has not been possible to directly apply the damage accumulation theories (that have been originally developed and used for metallic materials) due to these complexities. It has been customary to approach the problem at the macroscopic (i.e., at the laminate) level encompassing the structural response. Moreover, in the case of metallic materials, fatigue damage results in the reduction of strength properties of the materials. But in the case of composite laminates, it has been established that in addition to strength properties, the sti€ness properties, particularly the laminate elastic modulus, is more dominantly a€ected by fatigue. Therefore, a structural-level response parameter, which is sensitive to various damage modes has been established and is in use to describe and quantify the fatigue response. In particular, it has been established [14] that the compliance of the laminate, which is the reciprocal of the elastic modulus of the laminate, that corresponds to the type of loading such as axial or bending loading, is a reliable fatigue response parameter for laminated composites. As the fatigue loading is applied, the compliance of the laminate has been shown to increase. Therefore, experiment is conducted to obtain the compliance versus number of load cycles curve, for use in fatigue design, safety assessment, reliability studies, etc. 1.1. Randomness in fatigue response The behavior and engineering properties of metallic and composite materials are signi®cantly a€ected by the technological process defects, design errors at both material and structural level, and environmental e€ects, such as the e€ects of operating temperatures and hydrothermal e€ects, which are random in nature. The relative occurrences and the contribution to the overall failure of various damage mechanisms and failure modes cannot be precisely quanti®ed. For instance, the case of composite materials can be cited. Laminated composites display considerable variabilities in their mechanical properties. These are due to variations in properties in ®bers, in matrices, in the nonuniformity of bonding at interfaces, in void content, in ®bre orientation, distribution and volume fraction, in manufacturing processes etc. Various damage mechanisms and failure modes become active at load cycle values that vary from specimen to specimen. Signi®cant variations in other complexities (that are mentioned above) have also been observed, and these variations cannot precisely be quanti®ed. As a result, signi®cant randomness exists in fatigue

response of materials and laminated composites. Tests on a single specimen yield a de®nite value for each material parameter. But when a number of specimens are tested, the parameter values randomly ¯uctuate from specimen to specimen. Hence, ecient and accurate modeling of the fatigue response through the corresponding test data on the changes in the macro-level fatigue response parameter has not been possible. This has resulted in disagreement between theory and experimental data. 1.2. Existing approaches and their limitations For quantifying the fatigue behavior of materials, so far the deterministic approach has in many cases been employed. This deterministic approach has resulted in signi®cant disagreement between theory and practice, and hence their practical applications are very limited. Recent developments led to statistical and probabilistic techniques. Weibull [15] introduced into fatigue failure studies a probability distribution function of extreme value type, in order to account for the variations in the life data of mechanical components made of metallic materials. Reliability analyses have been developed for metallic structures based on this distribution function. The Weibull distribution function has also been used to model the fatigue life data of composites [9,11,14]. The Weibull distribution function frequently provides a reasonable ®t to the probability distributions of life times that have been obtained from fatigue tests or service operation. Some aspects of variability in lifetime can be assessed and accounted for. But the probability distributions of damage parameters that are related to the strength and sti€ness properties, could not be reasonably represented by Weibull distribution in many cases. Moreover, there are only two limiting states of the material in the Weibull distribution based fatigue model, that are: (i) satisfactory (no damage) and (ii) failure (completely damaged). This imposes severe limitations on the validity and practical applications of this model. For instance, one can not incorporate the initial damage distribution (for instance in a composite laminate, this refers to the ¯uctuations in the initial elastic modulus), nor it is possible to trace the damage accumulation during service using the Non-Destructive Testing (NDT) techniques. Considering the above mentioned needs, a number of cumulative fatigue models have been developed for metallic materials and for materials such as reinforced concrete during the recent past [2,3,6]. In these models, (i) the crack length was employed as fatigue response parameter, and (ii) these model parameters are estimated from test data using statistical parameter estimation procedures. The true distributions (both individual and joint distributions) of fatigue response

R. Ganesan / Computers and Structures 76 (2000) 517±531

parameter, as revealed by the sample test data, have not been considered and incorporated in these models. In some cases, the Gaussian distribution was assumed. However, fatigue response parameters such as the crack length (in the case of metallic materials) and the compliance (in the case of composite laminates) do not follow a Gaussian distribution. This has been shown in many past studies in the literature, and the earlier studies of the present author [7,8]. In the present work, a stochastic approach to model and analyze test data on fatigue response of materials and composite laminates is developed. This approach is `data-driven' in nature and further, it can take into account the true probability distribution of the fatigue response parameter that is relevant and appropriate to the material being studied. 2. Markov process modeling of fatigue response parameter The modeling developed here considers the `Cumulative Damage (CD)' due to fatigue, which is the irreversible accumulation of fatigue damage in the material. The central idea of the present modeling is that in a suitably de®ned `Duty Cycle (DC)' a random amount of fatigue damage is accumulated in the material. This fatigue damage in turn results in the corresponding change in the fatigue response parameter, such as the increase in the crack length in the case of metallic materials and the increase in the laminate compliance (i.e., the reciprocal of elastic modulus) in the case of composite laminates. Here, the term duty cycle refers to a repetitive period of operation in the life of a material in which fatigue damage can accumulate. Test data regarding the fatigue response parameter (that corresponds to the particular loading applied onto the test specimens) are collected using material specimens that are subjected to fatigue loading. This fatigue response parameter is denoted here by C and the number of load cycles by n. The values of the stress ratio and the stress amplitude are denoted by R and S respectively. It may be noted here that the selection of a suitable fatigue response parameter appropriate to a particular material is not the main concern of the present work. Studies on material properties and behavior deal with this question, and in these studies for a particular material, an appropriate fatigue response parameter will be selected, evaluated and validated. Considerations regarding the nature of loading such as tensile or compressive or a combination of both, e€ects of load orientation relative to the cracks, suitability of the parameter for experimental work, etc. will be used in these studies as the basis for this purpose. Having selected a particular fatigue response parameter, a

519

number of test coupons will be fabricated and tested in accordance with the relevant ASTM standard and test procedure. For instance, in the case of composite laminates the ASTM standard D3039/D3039-93 has to be followed. For each test coupon, the fatigue response parameter (for instance in the case of composite laminates, the compliance) will be determined at di€erent load cycles using load-displacement measurements. The fatigue response parameter is denoted here by C and the number of load cycles is denoted by n. For each test specimen, the C±n curve can be plotted. Due to the inherent variability associated with the fatigue response, this curve varies from specimen to specimen. Having obtained a set of curves (corresponding to a number of test coupons) relating the fatigue response parameter and the number of load cycles, the stochastic modeling for fatigue response can be performed. This modeling methodology is developed and described in the following. The modeling methodology is more general in nature. It can be used to model the fatigue response of any material; in that case the fatigue response parameter C will be the one that has been established through material studies to be appropriate for that particular material. For any particular value of the number of load cycles n1, the values of fatigue response parameter C…n1 † corresponding to each and every test specimen, constitute a random variable, denoted by C1. For any two values of the number of load cycles n1 and n2, the corresponding values of C…n1 † and C…n2 † constitute two correlated random variables, denoted by C1 and C2, respectively. These two random variables are related by their joint probability distribution function. In a similar manner, the sets of compliance values at various load cycles constitute correlated, multi-variate, random variables which are described by multi-variate joint probability distributions. Hence, the set of C±n curves corresponding to each and every test specimen, essentially constitutes a stochastic process. For another fatigue loading with di€erent values of stress ratio and stress amplitude, a similar stochastic behavior can be observed but the magnitudes of fatigue response parameter are di€erent. Hence, the stochastic process modeling the fatigue response parameter is denoted as fC…n,S,R†,n 2 N,S 2 A,R 2 B g, where N, A and B are the respective sample spaces of the number of load cycles, the stress amplitude and the stress ratio. The sets N, A and B are the index sets of the stochastic process fC…n,S,R†g: A discrete set of damage states is determined based on the following procedure. The entire range of fatigue response parameter values is calculated from the sample test data. This range is then divided into m number of equal or non-equal intervals, …C2 ÿ C1 †,…C3 ÿ C2 †, . . . ,…Ci‡1 ÿ Ci †, . . . ,…Cm‡1 ÿ Cm †, where Cm‡1 is the maximum value of the fatigue response parameter.

520

R. Ganesan / Computers and Structures 76 (2000) 517±531

If the fatigue response parameter of a test specimen at any load cycle r is within the range …Ci‡1 ÿ Ci ), then the material specimen is said to be in the damage state i when the number of load cycles is r. When the number of load cycles is further increased, the specimen performs a transition from one damage state to another until the ®nal damage state, i.e., the fatigue response parameter being within the range …Cm‡1 ÿ Cm ), is reached. Due to the inherent variability in the material and damage properties of the material, the transition paths vary from specimen to specimen. Accordingly, they are probabilistically quanti®ed by the `Transition Probabilities'. The probability of any damage state is mathematically represented using the Conditional Probability Density Function as  pij ˆ P C…nj † ˆ jjC…njÿ1 † ˆ i,C…njÿ2 †  ˆ i ÿ 1, . . . ,C…0 † ˆ 1 …1† where …n0 ,n1 ,n2 , . . . ,njÿ1 ,nj , . . . ,nmax † is the index set N and …1,2,3, . . . ,j ÿ 1,j, . . . ,m† is the state space of the fatigue response parameter (stochastic) process. Here, nmax is the maximum number of load cycles and m is the ®nal damage state that is pre-set by the analyst. 2.1. Finite Markov chain The value of fatigue response parameter at the end of any DC fatigue loading depends, more or less, on both, the value of the fatigue response parameter at the start of that DC and the change in fatigue response parameter caused by the damage induced during the DC fatigue loading. So, the fatigue response parameter process is a Markov process. The fatigue response parameter process is non-homogeneous because the fatigue response parameter is a continuous function of the number of load cycles. The index set N of this nonhomogenous Markov process is …n0 ,n1 ,n2 , . . . ,njÿ1 ,nj , . . . ,nmax † and the state space is …1,2,3, . . . ,j ÿ 1,j, . . . ,m). Both these index sets are discrete. A Markov process that has discrete index sets is called Markov Chain (MC). Further, since the index sets are ®nite, the fatigue response parameter process becomes a Finite Markov Chain.

specimen entering a new damage state j at loading cycle nk , given that it is currently in damage state i at loading cycle nkÿ1 : Based on the information provided above as to the classi®cation of damage states of the material specimen, Eq. (2) for the TPF can be recast in the following form: pij …nkÿ1 ,nk †    ˆ P Cj < C…nk †RCj‡1 j Ci < C…nk †RCi‡1

…3†

Based on the de®nition of the conditional probability density function, joint probability density function and Bayes' theorem [5], this equation can be recast into pij …nkÿ1 ,nk † ˆ

    P Cj < C…nk †RCj‡1 \ Ci < C…nk †RCi‡1   P Ci hC…nk †RCi‡1

…4†

Now, it can be observed that both the probability distribution of the fatigue response parameter at load cycle nkÿ1 as well as the joint probability distribution of fatigue response parameter at load cycles nkÿ1 and nk , are required as the input information, in order to evaluate the TPF between damage states i and j, pij …nkÿ1 ,nk †: Since the fatigue response parameter at any load cycle is a random variable, the fatigue response parameters at load cycles nkÿ1 and nk are essentially two correlated random variables Ckÿ1 and Ck , respectively. The individual probability density functions of random variables Ckÿ1 and Ck are denoted by fCkÿ1 …ckÿ1 † and fCk …ck †, respectively. The joint probability density function corresponding to these two correlated random variables is denoted by fCkÿ1 Ck …ckÿ1 ,ck †: In terms of these probability density functions, the TPF between damage states i and j, can be determined according to … Cj‡1 … Ci‡1 pij …nkÿ1 ,nk † ˆ

Cj

C

( …i

fCkÿ1 Ck …ckÿ1 ,ck † dckÿ1 dck )

Ci‡1 Ci

…5†

fCkÿ1 …ckÿ1 † dckÿ1

2.2. Transition probability function (TPF) As a result of the Markovian property and nonhomogeneity of the stochastic process, Eq. (1) is rewritten as   pij …nkÿ1 ,nk † ˆ P C…nj † ˆ jjC…njÿ1 † ˆ i ; nk > nkÿ1 …2† where pij …nkÿ1 ,nk † is the Transition Probability Function (TPF), which is the probability of the material

2.3. Transition probability matrix (TPM) Generalizing Eq. (2) for all possible damage states included in the state space, the Transition Probability Matrix (TPM) of the material specimen, denoted by ‰P…nkÿ1 ,nk †], when the load cycle is increased from nkÿ1 to nk , can be obtained. The TPM corresponding to the load cycles nkÿ1 and nk is thus given as

R. Ganesan / Computers and Structures 76 (2000) 517±531



P…nkÿ1 ,nk † 2



p11 …nkÿ1 ,nk † 6 p21 …nkÿ1 ,nk † 6 ˆ6 6
4
pm1 …nkÿ1 ,nk †

p12 …nkÿ1 ,nk † p22 …nkÿ1 ,nk †

pm2 …nkÿ1 ,nk †















3

p1m …nkÿ1 ,nk † p2m …nkÿ1 ,nk † 7 7 7
7 5
pmm …nkÿ1 ,nk † …6†

where m is the number of damage states in the state space. Since all damage states are mutually exclusive and collectively exhaustive, the sum of the probabilities in any row of the TPM that has non-zero probabilities is equal to unity: j X pij …nkÿ1 ,nk † ˆ 1 i,k ˆ 1,2,3, . . . ,m;nk > nkÿ1

…7†

When the load cycle is nkÿ1 , the fatigue response parameter is C…nkÿ1 † ˆ i: Now, if the load cycle is increased to nk , the fatigue response parameter is C…nk † ˆ j: It is also possible that C…nk † ˆ C…nkÿ1 † so that j ˆ i: However, the system can not make a transition from a higher state to a lower state, i.e., C…nk † < C…nkÿ1 † so that jhi: Also, the ®nal damage state m is the so-called `absorbing state'. The specimen which has entered this damage state can not make a transition to any other future state and so pmm …nkÿ1 ,nk † ˆ 1: Imposing these constraints, the TPM can be written as   P…nkÿ1 ,nk † 2 3 p11 …nkÿ1 ,nk † p12 …nkÿ1 ,nk †

p1m …nkÿ1 ,nk † 60 p22 …nkÿ1 ,nk †

p2m …nkÿ1 ,nk † 7 6 7 7 ˆ6 0 0


6 7 40 5 0


0 0

1 …8†

2.4. M-step transition probability matrix The transition probability matrices corresponding to di€erent sets of load cycles n0 and n1, n1 and n2 , . . . ,nmÿ1 and nm , can be determined following the procedure outlined above. They are denoted by ‰P…n0 ,n1 †Š,‰P…n1 ,n2 †Š,‰P…n2 ,n3 †Š, . . . ,‰P…nmÿ1 ,nm †Š, respectively. The transition probability matrix of the material specimen corresponding to any initial load cycle n0 and any ®nal load cycle nm can be determined in terms of these transition probability matrices, based on the stochastic process theory. The resulting matrix of order m  m is called as the `M-Step Transition Probability Matrix (M-Step TPM)' and is denoted by ‰P…n0 ,nm †Š: In order to determine this matrix, use is made of the Chapman±Kolmogorov Theorem [5].

521

Considering the fact that the Markov chain considered here is non-homogeneous, and based on this theorem, the product equation for the M-Step TPM is obtained as         P…n0 ,nm † ˆ P…n0 ,n1 †  P…n1 ,n2 †  P…n2 ,n3 †   


 P…nmÿ1 ,nm † …9† It can now be observed that both the Transition Probability Matrix ‰P…nkÿ1 ,nk †Š as well as the M-Step Transition Probability Matrix ‰P…n0 ,nm †Š consist of conditional probabilities corresponding to various damage states that the specimen can assume. They depend on the initial damage state of the specimen from which the transition is taking place. This is so because, (i) the TPM, ‰P…nkÿ1 ,nk †], consists of the probabilities of ®nding the specimen in any of the damage states 1, 2, 3, . . . , m at load cycle nk , given the present state of the specimen at load cycle nkÿ1 , and (ii) the M-Step TPM, ‰P…n0 ,nm †Š consists of the probabilities of ®nding the specimen in any of the damage states 1, 2, 3, . . . , m at load cycle nm , given the present state of the specimen at load cycle n0. In order to determine the unconditional probabilities of the damage states of the specimen, given the initial damage state distribution, when the load cycle is increased from n0 to nm , the initial state vector of the specimen is considered. Based on this, and from Eq. (9), the following equation is derived. n o      PU …n0 ,nm † ˆ P…0 †  P…n0 ,n1 †  P…n1 ,n2 † …10†   


 P…nmÿ1 ,nm † In this equation, fP…0†g denotes the initial damage state (row) vector of the specimen which consists of the probabilities of the specimen in di€erent damage states at load cycle n0 : Further, fPU …n0 ,nm †g is the damage state (row) vector consisting of the unconditional probabilities of ®nding the specimen in di€erent damage states at loading cycle nm : 2.5. Analytical probability density functions The evaluation of transition probabilities according to Eq. (5) requires the probability density function of the fatigue response parameter at load cycle nkÿ1 , which is denoted as fCkÿ1 …ckÿ1 †: This is determined, based on the Maximum Entropy Method [13]. Based on this method, the probability density function fCkÿ1 …ckÿ1 † is expressed as # " p X i fCkÿ1 …ckÿ1 † ˆ exp l0 ‡ …11† li ckÿ1 iˆ1

522

R. Ganesan / Computers and Structures 76 (2000) 517±531

The values of the parameters li , i ˆ 0,1,2, . . . ,p are determined such that the entropy of the random variable ckÿ1 given by …   S ˆ ÿ fCkÿ1 …ckÿ1 †ln fCkÿ1 …ckÿ1 † dckÿ1 …12† R

where R denotes the range of the random variable ckÿ1 , is maximum, subject to the following constraints: (i) the total probability is unity, i.e., … fCkÿ1 …ckÿ1 † dckÿ1 ˆ 1 …13† R

and (ii) the probabilistic moments of the random variable ckÿ1 , de®ned in terms of the probability density function fCkÿ1 …ckÿ1 † are equal to the corresponding statistical moments obtained from the compliance data, i.e., … cikÿ1 fCkÿ1 …ckÿ1 † dckÿ1 ˆ Mi , i ˆ 1,2, . . . ,q …14†

2. Each set of fatigue response parameter values is represented as a variable and so, as an abscissa in a three-dimensional plot. Two variables are de®ned this way. Each abscissa is then divided into a number of intervals, equal or unequal. A number of rectangular areas are formed as a result, on the surface de®ned by the two abscissae. 3. The relative frequencies that correspond to each and every rectangular areas are calculated from the two sets of fatigue response parameter values. Over a particular rectangular area, the value of relative frequency is assumed to be constant. The number of rectangular areas are to be determined by judgement. Since the joint density functions are to be used later for the calculation of transition probabilities using Eq. (5), it is preferable to de®ne the areas using the bounding fatigue response parameter values of di€erent damage states.

R

where Mi is the ith statistical moment about the origin of the compliance values ckÿ1 : Further details regarding the Maximum Entropy Method can be found in [13]. The above problem should be solved based on the techniques of optimization. The optimization technique that is best suited to this problem is the Nonlinear Programming Technique proposed by Jacobson and Oksman [12]. 2.6. Joint probability density functions The evaluation of transition probabilities according to Eq. (5) also requires the joint probability density function of the fatigue response parameters at load cycles nkÿ1 and nk (which are basically two correlated random variables) denoted as fCkÿ1 Ck …ckÿ1 ,ck †: The concept of maximum entropy that has been deployed for the determination of individual probability density functions can also be applied for the determination of joint probability density functions. This would however involve the optimization of a two-dimensional surface subject to constraints given by the joint probability moments, and would require a large data set. The data set obtained here with 35 samples is not adequate to perform such an optimization exercise. So, an alternative way of determining an estimation of joint probability density function based on the relative frequency approach [1] is deployed. This procedure aims at constructing a two-variable histogram as per the following procedure: 1. The ranges of fatigue response parameter values that correspond to two load cycle values are determined from the test data.

3. Estimation of the distribution of fatigue response parameter Now, the procedure for estimating the probability that the material specimen will enter a certain damage state, for instance a certain range of values of compliance in the case of composite laminate, at a given load cycle is detailed. The initial and ®nal values of load cycles that correspond to the test data are ®rst considered. This duty cycle is expressed in terms of m number of smaller duty cycles. For each small duty cycle, the transition probability matrix is determined according the procedure given in the foregoing. This way, m number of transition matrices will be available. For each entry in the transition probability matrix, m number of values will be available. For instance, consider any transition probability pij : If say 10 steps (smaller duty cycles) are considered, 10 values of this transition probability viz., 1 pij ,2 pij ,3 pij , . . . ,10 pij will be available. These 10 values are considered as constituting a random variable. Further, it is assumed here that this random variable, denoted by pij has a Gaussian distribution. The true distribution of each transition probability should be determined using the Maximum Entropy Method described in the previous section. However, at least 35 sample values are required for applying this method. If only few values are available, the true distribution can not be obtained using the limited data. To circumvent this diculty, Gaussian distribution is assumed. The 10 sample values are used as a basis for predicting the value of transition probability pij , that corresponds to the duty cycle de®ned by the ®nal load cycle value of the 10th (smaller) duty cycle and the value of the load cycle at which the estimation is to be made.

R. Ganesan / Computers and Structures 76 (2000) 517±531

The relevant rationale and procedure are outlined as follows. 1. For the random variable pij , a Gaussian probability density function is ®tted using the 10 sample values available. The parameters of this Gaussian density function, mij and sij , can be determined following the standard parameter estimation techniques [1]. 2. A standardized random variable zij is de®ned according to the relationship zij ˆ …pij ÿ mij †=sij : This standardized random variable will also have Gaussian distribution because the Gaussian nature of probability distribution does not change when a linear transformation of normal random variable is performed. 3. Now, the concepts of reliability theory and failure probabilities [13] are deployed. It can be proved using the probability theory and system reliability theory that a value of pij that corresponds to a reliability R (say 90%) can be obtained by ®rst determining the value of z90 , such that the area bounded by the density function curve of zij between the range ÿINF to z90 is 0.9, and then calculating the value of pij according to the relationship …z90 sij † ‡ mij : For instance, when the reliability is 90%, the value of z90 is approximately equal to 1.3. Now consider a transition probability, say p33 : By ®tting a Gaussian distribution, the parameters of Gaussian

523

density function for p33 can be obtained and say, they …m33 and s33 † are respectively equal to, 0.4819 and 0.2843 (it is worthwhile to note the magnitude of randomness in the values of p33 , which is approximately 60% in the above example). The value of p33 which will have a reliability of 90% (i.e., the probability that this value will be exceeded in actual test results is 10%) is {[1.3 (0.2843)] +0.4819} which is equal to 0.8515. This is the predicted value of p33 : This procedure is followed to determine the predicted values of all the transition probabilities. This way, the complete transition probability matrix is determined. Now, the (10+1)-step duty cycle de®ned by the initial value of the load cycle in the test program and the value of the duty cycle at which the estimation is to be made, is considered. Ten transition probability matrices can be determined from test data and the 11th transition probability matrix is predicted. The unconditional probabilities of damage states corresponding to the load cycle at which estimation is to be made can be determined using Eqs. (9) and (10). The value of m is equal to 11. These unconditional probabilities de®ne the histogram of damage states corresponding to the load cycle where estimation is to be made.

Fig. 1. Typical sample realizations of specimen compliance.

524

R. Ganesan / Computers and Structures 76 (2000) 517±531

4. Application to a composite material The modeling and estimation procedure described in the above is now illustrated. As an example, the case of a composite material is considered. The composite material considered is AS4/3501-6 Graphite Epoxy cross-ply layup [0/90]4s. Experimental investigation on the fatigue behavior of this composite material has been performed as reported in [7,8]. Details regarding the manufacturing of samples [4,10], experimental setup, microstructure observations, measurement, etc. were reported in these works and hence, they are not repeated here. In essence, an on-axis T±T fatigue load at 5 Hz frequency and with a load ratio R of 0.1 was applied; the compliance of the laminate specimen, that corresponds to the loading direction, is taken to be the fatigue response parameter. The values of the compliance at di€erent load cycles were experimentally determined and were reported in these works. The relevant test data as reported in [7,8] are plotted here again as in Fig. 1. 4.1. Transition probabilities The transition probability matrix is evaluated considering a certain duty cycle that is de®ned by two values of load cycles, one at the beginning of the duty cycle and another at the end of the duty cycle. These two values of load cycles are for instance say 0 and 101. The transition probability matrix is determined using Eq. (5), based on the individual probability density function of the compliance corresponding to the initial value of load cycles i.e., zero, which is determined using the maximum entropy method, and the joint probability density function of compliances corresponding to the initial and ®nal values of load cycles, i.e., 0 and 101. Accordingly, nkÿ1 ˆ 0, nk ˆ 101, fCkÿ1 Ck …ckÿ1 ,ck † ˆ fC0 C101 …c0 ,c101 † and fCkÿ1 …ckÿ1 † ˆ fC0 …c0 † in Eq. (5). The equation for any entry in the transition probability matrix can now be written. For instance, the transition probabilities p11 and p12 can be written as p11 …n0 ,n101 † ˆ

… C2 … C2 C1 C1 … C2 C1

p12 …nkÿ1 ,nk † ˆ

fC0 C101 …c0 ,c101 † dc0 dc101

… C3 … C2 C2 C1 … C2 C1

…15†

fC0 …c0 † dc0

fC0 C101 …c0 ,c101 † dc0 dc101

…16†

fC0 …c0 † dc0

In the above two equations, C1 and C2 refer to the bounding values of compliances that de®ne the damage

state 1. In a similar manner, C2 and C3 de®ne the bounding values of compliances that de®ne the damage state 2. As can be seen from the above two equations, the numerator of the equation for any transition probability is obtained through the double integration of the joint density function corresponding to the compliances at initial and ®nal values of the duty cycle (which are 0 and 101 in the above example) over the bounding values of compliances that de®ne the initial and ®nal damage states. The denominator is obtained through the integration of the individual density function of the compliance corresponding to the initial value of the duty cycle (which is zero) over the bounding values of compliances that de®ne the initial damage state. This way, the equations for all the upper-diagonal entries of transition probability matrix can be written. A computer program is written to perform the calculation of the transition probabilities as per this procedure. It has modules to perform the numerical integration of individual density function and the double integration of the joint density function. Once the transition probability matrix for a certain duty cycle is determined, the calculation proceeds to the next duty cycle. After m number of duty cycles have been considered, m number of transition probability matrices will be available. Based on these matrices, the m-step transition probability matrix can be determined, which corresponds to the initial value of the ®rst duty cycle and the ®nal value of the last duty cycle, that are included in the m-step duty cycle. The determination of the m-step transition probability matrix corresponding to the 5-step duty cycle with initial value of 0 and ®nal value of 301 is now illustrated in several sequential calculation steps. The ®ve smaller duty cycles that de®ne the duty cycle (0, 301) are: (0, 101), (101, 151), (151, 201), (201, 251) and (251, 301). The fatigue damage induced during the duty cycle (0, 301) is represented by 17 damage states, that are equal in size. Five transition probability matrices are obtained using six sets of compliance values corresponding to the following six load cycle values: 0, 101, 151, 201, 251 and 301. The evaluation of transition probability matrices is performed in the following sequence. Step 1 (Determination of the damage states). The largest value in the state space i.e., the largest value of the compliances of all test specimens, is 1.7831  10ÿ11 Paÿ1 and the lowest value is 1.6114  10ÿ11 Paÿ1, as can be seen from Fig. 1. This range of compliance is then divided into 17 damage states of same size. Step 2 (Generate probability density functions for each compliance value set C…nkÿ1 †). The single-variable distributions are estimated using the maximum entropy method from the data sets C…nkÿ1 †, k ˆ 1,2, . . . ,5: The

R. Ganesan / Computers and Structures 76 (2000) 517±531

525

Table 1 List of the values of li at ®ve load cycles

0 cycle 101 cycles 151 cycles 201 cycles 251 cycles

l0

l1

l2

l3

l4

ÿ59.75927 ÿ52.95307 ÿ10.20866 11.63882 ÿ1.021048

12.11768 4.758957 0.032665 ÿ1.913443 ÿ0.437083

ÿ0.976104 ÿ0.163753 0.0168172 0.0765457 0.0163917

0.0350279 0.0024078 ÿ0.000467 ÿ0.001223 ÿ0.000210

ÿ0.000468 ÿ0.000013 0.0000034 0.0000067 0.0000008

selection of bounds for the compliance values …Akÿ1 ,Bkÿ1 † in the density function fC kÿ1 has to be based on judgment. Also, the bounds must be selected such that the sum of probability densities in any row of the transition matrix that has non-zero entries must be theoretically unity, which however for practical purposes may be within the 10% tolerance range. The ®rst four statistical moments, i.e., M1, M2, M3 and M4, of each compliance value set C…nkÿ1 † are determined and are used as constraints in the determination of fC kÿ1 : The maximum entropy method yields the values of the ®ve parameters viz., li , i ˆ 0,1,2, . . . ,4, for each fC kÿ1 : The values of these parameters are listed in Table 1. Step 3 (Determination of joint probability density functions). The bounding values of compliances for the 17 damage states have to be used to de®ne the rectangular regions, for which the corresponding relative frequencies are determined from the test data. Five histograms that correspond to ®ve smaller duty cycles (0, 101), (101, 151), (151, 201), (201, 251) and (251, 301) are thus determined. Step 4 (Evaluation of transition probability matrices).    P…n0 ,n1 † ˆ P…0,101† ˆ Row1 Row2

The transition probability matrices can now be evaluated based on the single-variable distributions obtained in step 2 and the 3D histograms obtained in step 3, using Eq. (5). The single-variable integration involved in Eq. (5) is performed using numerical integration. The integral of the double integration in Eq. (5) is the volume directly above the rectangular area bounded by …Ci ,Ci‡1 † and …Cj ,Cj‡1 † in the histogram. At this stage, the sum of the non-zero transition probabilities in each row of the transition probability matrix is calculated. Theoretically, the sum must be equal to unity. However, for practical purposes, a tolerance margin of 210% is used in the present work. If the sum of the transition probabilities is more than 1.1 or less than 0.9, suitable modi®cations have to be made as to the bounding compliance values of the damage states Ci , i ˆ 1,2, . . . ,18 and/or the bounds used for the singlevariable distributions. An iterative (trial and error) procedure is hence called for. The ®ve transition probability matrices, each of order 17  17 corresponding to the ®ve duty cycles (0, 101), (101, 151), (151, 201), (201, 251) and (251, 301) are obtained using the above procedure. Parts of the rows of these ®ve matrices that have non-zero entries are given below:



0:0 0:4369 0:0 0:2003

0:2496 0:1872 0:1248 0:5342 0:1335 0:0668

0:0 0:0668

0:0 0:0

   P…n1 ,n2 † ˆ P…101,151† ˆ Row2 Row3 Row4 Row5 Row6

0:0 0:0 0:0 0:0 0:0

0:3064 0:0 0:0 0:0 0:0

0:7149 0:5243 0:0 0:0 0:0

0:0 0:5243 0:3935 0:0 0:0

0:0 0:0 0:5902 0:2792 0:0

0:0 0:0 0:0 0:5583 0:9089

0:0 0:0 0:0 0:0 0:0

   P…n2 ,n3 † ˆ P…151,201† ˆ Row2 Row3 Row4 Row5 Row6

0:0 0:0 0:0 0:0 0:0

0:8249 0:0 0:0 0:0 0:0

0:0 0:7325 0:0 0:0 0:0

0:0 0:4578 0:6279 0:0 0:0

0:0 0:0 0:2093 1:0321 0:0

0:0 0:0 0:0 0:0 0:9968

0:0 0:0 0:0 0:0 0:0





526

R. Ganesan / Computers and Structures 76 (2000) 517±531

   P…n3 ,n4 † ˆ P…201,251† ˆ Row2 Row3 Row4 Row5 Row6

0:0 0:0 0:0 0:0 0:0

1:001 0:0 0:0 0:0 0:0

   P…n4 ,n5 † ˆ P…251,301† ˆ Row2 Row3 Row4 Row5 Row6

0:0 0:0 0:0 0:0 0:0

0:9431 0:0 0:0 0:0 0:0



0:0 0:6418 0:0 0:0 0:0

0:0 0:3851 0:6094 0:0 0:0

0:0 0:0 0:3482 1:0036 0:0

0:0 0:0 0:0 0:0 1:0921

0:0 0:0 0:0 0:0 0:0



0:0 0:5112 0:0 0:0 0:0

0:0 0:3408 1:0379 0:0 0:0

At this stage, the sum of the non-zero transition probabilities in each row of the transition probability matrices is not exactly equal to unity, but within the tolerance range. The algebraic di€erence between the row sum and unity (total probability), which is the net computational error, is distributed to each entry in the row, and the transition probabilities are this way adjusted, so as to be proportional to the distributed error. Step 5 (Calculation of m-step transition probability matrices). Once the ®ve transition probability matrices are evaluated corresponding to the ®ve smaller duty cycles, the 5-step transition probability matrix ‰P…n0 ,n5 †Š, which is ‰P…0,301†Š can be determined using Eq. (9). For illustration purposes, the relevant equation is given as 

       P…n0 ,n5 † ˆ P…0,301† ˆ P…0,101†  P…101,151†        P…151,201†  P…201,251†  P…251,301†

0:0 0:0 0:0 0:2261 0:5762

0:0 0:0 0:0 0:0 0:2881

initial damage states of the specimens and is determined from the distribution of C(0). That is, the compliances of all specimens at the beginning of the test before any duty cycle loading is applied, are considered. The number of samples, the initial compliances of which belong to the ®rst damage state is determined. The ratio of this number of samples over the total number of samples is calculated, and this ratio is equal to the ®rst entry of fP…0†g: This procedure is repeated considering each and every damage states. This way, all entries of the row vector fP…0†g are determined. For the present case, the row vector fP…0†g is given by 

P…0 † ˆ f0:516129,0:483871,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0g

The unconditional probabilities of the compliances at n5 ˆ 301 is given below: n

The entries in the ®rst two rows and ®rst six columns of the transition probability matrices have non-zero values and they are given below:    P…n0 ,n5 † ˆ P…0,301† ˆ Row1 Row2

0:0 0:0 0:0 0:9045 0:0

o    PU …n0 ,n5 † ˆ P…0 † x P…0,301† ˆ f0:0, 0:0967 0:0967 0:3870 0:2580 0:1290 0:0322 0


0g



0:0 0:0

0:4369 0:2003

Step 6 (Calculation of the unconditional probabilities of ompliances at load cycle 301). This last step is carried out using Eq. (10) to determine the probability density function of compliances at the end of duty cycle (0, 301), which in turn is a function of both the initial distribution of compliances fP…0†g and the 5step transition probability matrix ‰P…0,301†Š: The row vector fP…0†g consists of the probability distribution of

0:2496 0:1872 0:1248 0:5342 0:1335 0:0668

0:0 0:0668

The histogram corresponding to these unconditional probabilities is presented in Fig. 2. Now, this result can be compared with the actual distribution of C…n5 † at n5 ˆ 301 cycles, which is plotted in Fig. 3. Probabilities pertinent to the damage states are given below: f0:0, 0:0967 0:0967 0:3870 0:2580 0:1290 0:0322 0


0g

Fig. 2. The unconditional probabilities of compliance at load cycle 301.

Fig. 3. Histogram obtained from Fig. 1 of compliance at load cycle 301.

528

R. Ganesan / Computers and Structures 76 (2000) 517±531

Fig. 4. Predicted probability distribution for the compliance at load cycle 351.

As can be seen from the above, the probability distribution of the compliances at load cycle 301 determined analytically by the proposed model compares very well with the distribution of the actual data.

The unconditional probability matrix fPU …0,351†g is determined according to the following equations: 

     P…0,351† ˆ P…0,101†  P…101,151†

4.2. Estimation of damage state probabilities at load cycle 351 Now, the duty cycle (301, 351) is considered. Using the transition probability matrices ‰P…0,101†Š, ‰P…101,151†Š, ‰P…151,201†Š, ‰P…201,251†Š, and ‰P…251,301†Š, the transition probability matrix ‰P…301,351†Š is predicted with a reliability of 90%. All the non-zero transition probabilities are located within the ®rst six rows and six columns of predicted matrix ‰P…251,301†Š, and they are given below:   P…301,351† ˆ Row1 Row2 Row3 Row4 Row5 Row6

0:0 0:0 0:0 0:0 0:0 0:0

0:3414 1:1415 0:0 0:0 0:0 0:0

0:1950 0:7022 0:8515 0:0 0:0 0:0

     P…151,201†  P…201,251†      P…251,301†  P…301,351† 

   PU …0,351† ˆ P…0 †  P…0,351†

The (predicted) unconditional probability matrix that corresponds to the load cycle 351 is given by fPU …0,351†g:

0:1463 0:1043 0:6060 1:0259 0:0 0:0

0:0975 0:0522 0:0 0:5545 1:2590 0:0

0:0 0:0522 0:0 0:0 0:4752 1:2924

R. Ganesan / Computers and Structures 76 (2000) 517±531



0:0 0:0709

0:0965 0:2989 0:3494

0:0899 0:0944 0 0


0

The corresponding histogram of compliance at load cycle 351 is plotted in Fig. 4. From the test data, the histogram of compliances at load cycle 351 is determined (Fig. 5) and the damage state probabilities are determined from this histogram. They are given below: 

0:0 0:0968

0:0968 0:3548 0:2581



fatigue damage in the composite laminate. From the transition probability matrix ‰P…0,101†Š, it can be seen that p11 ˆ 0:0, p12 ˆ 0:4369, p13 ˆ 0:2496, and so on. This means that, 43.69% of samples, the compliances of which lie in the range of 1:61  10ÿ11 ±1:62  10ÿ11 Pa before the application of duty cycle are so damaged

0:0968 0:0968 0 0


0

The damage state probabilities, that are obtained (i) by analytical estimation and (ii) from the test data, can now be compared. That the predicted values of damage state probabilities and the damage state probabilities obtained from test results are in reasonable agreement, can be observed. The transition probability matrices can be seen to both characterize as well as quantify the evolution of

529



that their compliance values increase after the application of 101 load cycles thereby reaching the range of 1:62  10ÿ11 to 1:63  10ÿ11 Pa. Out of the remaining samples, 24.96% of samples display a compliance increase to the range of 1:63  10ÿ11 to 1:64  10ÿ11 Pa, 18.72% of samples display a compliance increase to the range of 1:64  10ÿ11 to 1:65  10ÿ11 Pa, and 12.48% of samples display a compliance increase to

Fig. 5. Histogram of the compliance at load cycle 351.

530

R. Ganesan / Computers and Structures 76 (2000) 517±531

the range of 1:65  10ÿ11 to 1:66  10ÿ11 Pa. In the same manner, the transition probabilities in the second row can be interpreted. Further, it can be observed that: (i) none of the samples, the compliances of which lie in the range of 1:61  10ÿ11 to 1:62  10ÿ11 Pa at load cycle zero, retain the same compliance value after the application of duty cycle, thereby indicating that each and every sample gets damaged when 101 load cycles are applied, (ii) the same is not the case when the duty cycles (101, 151), (151, 201), (201, 251) and (251, 301) are considered, and in fact since p22 of ‰P…101,151†Š ˆ 30:64%, p22 of ‰P…151,201†Š ˆ 82:49%, and p22 of ‰P…201,251†Š ˆ 100%, it can be said that when the samples get `aged' due to `service', it is more likely that they resist vigorously, further damage (i.e., the laminated samples display more `inertial resistance' when the duty cycle `dynamic' loading through fatigue damage tries to move them from one damage state to another; or equivalently, the laminate sample behaves like the so-called `hard spring' which becomes sti€er as it gets deformed), (iii) 2.00% of samples the compliances of which lie in the range of 1:62  10ÿ11 to 1:63  10ÿ11 Pa at load cycle zero retain the same compliance even after the application of 101 load cycles, thereby indicating that the damage evolution process is sensitive to the initial conditions, i.e., the damage present before the application of load cycles, (iv) the maximum value of compliance encountered in any of the samples after the application of 101 load cycles is 1:67  10ÿ11 Pa. That the damage evolution is sensitive to the initial damage can be observed also from the transition probability matrices ‰P…101,151†Š, ‰P…151,201†Š, ‰P…201,251†Š, and ‰P…251,301†Š: The transition probabilities are distributed in a di€erent manner in di€erent rows of the transition probability matrices. Since each row corresponds to a particular range of compliance values, it can be observed that the magnitudes of fatigue damage setting in a set of samples, the initial compliances of which fall within the same range, are di€erent from that of another set of samples that belong to a di€erent range of initial compliance values. This trend emphasizes the sensitivity of damage evolution process during the application of a duty cycle to the initial damage. Any dynamic process that is sensitive to initial conditions is a non-linear dynamic process. Hence, the fatigue damage process in composite laminates is a non-linear process. This non-linearity hinders the application of Miner's rule which is more appropriate to linear damage processes, to the composite laminates.

5. Conclusions The fatigue response of metallic materials and laminated composites, is considered, and a stochastic approach to model and analyze the test data on fatigue response is developed. The following major sources of variability in cumulative fatigue are incorporated in the stochastic model developed in the present work: (i) initial state of the material and the variability in the initial values of the fatigue response parameters such as the compliance of composite laminates. Since the fatigue behavior is a non-linear process with respect to load cycles, changes in the initial values of the fatigue response parameter signi®cantly in¯uence the fatigue response process; (ii) severity and order of both the load cycles as well as the material failure (including the changes in material parameters); (iii) state of the material at the failure of the test specimen. The severity of load cycle and the material failure (including the changes in the material properties) is speci®ed in the model in terms of the probability transition matrix. The variations in the severity of repetitive load cycles and the material failure are incorporated through the non-homogeneous nature of the fatigue process. The randomness in the initial values of the fatigue response parameter due to manufacturing errors, technological process defects, design errors etc., is speci®ed by the initial damage state vector. The stochastic modeling and analysis developed in the present work can be employed to estimate the probability distribution of the fatigue response parameter corresponding to a given load cycle value. However, this estimation is limited to the particular material, the loading, the fatigue response parameter (that was pre-selected in the test program), the set of parameters and circumstances, the de®nition and the value of the stress ratio, the stress amplitude, the loading history etc. that correspond to the test data. In other words, the estimation corresponds to the probability distribution of the fatigue response parameter, that would have been calculated directly from the test data if the data (that is, the set of sample values of the fatigue response parameter) at that given load cycle value were experimentally obtained. Stochastic modeling and analysis of the sample test data on a CFRP composite material that was reported in [7,8], is then performed. Many observations as to the fatigue behavior of this CFRP laminated composite have been made. It may be noted that the model developed in this work is `data-driven' in nature, that is, the parameters of the model should be determined from test data. For instance, for a particular laminate con®guration, material system and fatigue loading, these parameters must be determined based on sample tests.

R. Ganesan / Computers and Structures 76 (2000) 517±531

Acknowledgements The support provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) for the research project is greatly appreciated. The author likes to thank Dr. Siyu Zhang and Dr. Moataz ElKarmalawy for their assistance in the work that is involved in Section 4. The author likes to thank also Professor Suong Van Hoa for his interest in and encouragement for the research project.

[8]

[9]

References [1] Bendat JS, Piersol AG. Random data: analysis and measurement procedures. New York: Wiley±Interscience, 1986. [2] Bogdano€ JL, Kozin F. Probabilistic models of cumulative damage. New York: Wiley, 1985. [3] Bolotin VV. Estimation of service life for machines and structures. New York: ASME Press, 1989. [4] Cardrick AW, Smith MA. An approach to the development of meaningful design rules for fatigue-loaded CFRP components. Composites 1974;5:96±100. [5] Cox DR, Miller HD. The theory of stochastic processes. London: Metheun, 1965. [6] Desayi P, Rao KB. Markov chain model for cracking behavior of reinforced concrete beams. Proceedings of ASCE: Journal of Structural Engineering 1989;115:2129± 43. [7] Ganesan R, Hoa SV, Zhang S, El Karmalawy M. A sto-

[10]

[11]

[12] [13] [14] [15]

531

chastic cumulative damage model for the fatigue response of laminated composites. In: The Proceedings of the Eleventh International Conference on Composite Materials (ICCM-11), July 14-18, Gold Coast, Australia. 1997. Ganesan R, El Karmalawy M, Hoa SV. Characterization of fatigue damage development and progression in composite laminates using stochastic modeling and analysis of test data. In: Proceedings of the 12th Annual Technical Conference of the American Society for Composites (ACM), October 6±8, Dearborn, Michigan. 1997. Hahn HT, Kim RY. Proof testing of composite materials. Journal of Composite Materials 1975;9:297±311. Ni RG, Adams RD. The Damping and Dynamic Moduli of Symmetric Laminated Composite Beams - Theoretical and Experimental Results. Journal of Composite Materials 1984;18:104±21. Ramani SV, Williams DP. Notched and unnotched fatigue behavior of angle-ply graphite/epoxy composites. In: Reifsnider KL, Lauraitis KN, editors. Fatigue of Filamentary Composite Materials, ASTM STP 636. 1977, pp. 27-46. Siddall JN. Optimal engineering design. New York: Marcel Dekker, 1992. Siddall JN. Probabilistic engineering design: principles and applications. New York: Marcel Dekker, 1983. Talreja R. Fatigue of composite materials. Lancaster: Technomic Publishing, 1987. Weibull W. A statistical distribution function of wide applicability. Transactions of the ASME: Journal of Applied Mechanics 1951;18:293±7.