Prognosis under uncertainty – An idealised computational case study

The object of this paper is to investigate computationally the possibility of damage prognosis in a structure under the most idealised circumstances. A simple isotropic material – Titanium alloy Ti-6Al-4V – is assumed. The structure is a simple finite plate under harmonic uniaxial loading and the damage is assumed to be a central mode 1 through crack for which various approximations to the stress intensity factor are known. The damage propagation model is the Paris-Erdogan law. Where the paper departs from complete simplicity is in the assumption that the parameters of the damage propagation law are uncertain. The paper investigates the effect of the parameter uncertainty on the estimated lifetime of the specimen. Two approaches are adopted for the uncertainty propagation, a statistical Monte Carlo scheme and one based on interval arithmetic.


Introduction
One of the most difficult problems in Structural Health Monitoring (SHM) is prognosis, i.e. given an accurate diagnosis of the damage state of a structure and a characterisation of the future loading to be experienced by the structure, what is its expected residual (safe) life?Under normal circumstances, this will depend critically on the soundness of the diagnosis and a good understanding of the underlying physics of damage progression.The existence of an accurate diagnosis is by no means assured; however, there is a great deal of accumulated experience in SHM within the engineering research community [1,2], and it is arguable that it has reached a point where it is possible to discuss prognosis.Similarly, while the existence of a high-fidelity damage accumulation model is not assured for complex materials and structures, considerable progress is being made in the development of such models using both theory and multi-scale physics modelling [3,4].The object of the current paper is to investigate damage propagation in a idealised situation in order to explore the feasibility of the approach.
The situation considered will be one where the structure is a finite rectangular metallic plate under uniaxial loading.The damage will take the form of a mode 1 opening crack situated centrally, perpendicular to the axis of loading.In the simple situation of a crack in a metallic plate, a useful approximation to the damage growth is given by the Paris-Erdogan (PE) Law.However, it is known that for metallic specimens, the parameters of the PE equation are subject to random variation between samples and even random variation at different points within a sample [5,6].This means that the computed lifetime of a given specimen will also be a random variable.Given experimental bounds on the variability of the PE equation parameters, the current study will investigate the extent of the variation in the predicted lifetime of the sample.In order to obtain the full statistics of the lifetime, it is possible to carry out a Monte Carlo analysis; however, this is often computationally intensive.It is made feasible here by the fact that the PE law only involves two parameters.For a general material or structure, it may well be that there are many parameters needed in order to characterise the damage propagation.It will also generally be the case that the damage propagation can not be summarised in a simple analytical formula, it may be necessary to execute a complex numerical model in order to investigate the damage growth.If the model is computationally complex and it takes a substantial amount of time to complete a single predictive run of the model, it will not generally be feasible to use a Monte Carlo approach.Alternative approaches to uncertainty propagation may prove advantageous in this case, and one approach which is notable for its speed and simplicity is interval arithmetic [7].
In the context of interval arithmetic, an object of the current paper is to take interval bounds for the parameters and solve the PE equation using interval arithmetic in order to obtain bounds on the lifetime.As interval arithmetic is known to be conservative, the approach will be compared with a Monte Carlo analysis to assess the degree of conservatism.It is found that conservatism is not a problem in the context of the PE equation and that an even simpler method of uncertainty analysis -vertex propagation -is sufficient to accurately bound the lifetimes.
The layout of the paper is as follows: Section Two will introduce the details of the computational model considered here.In Section Three, a Monte Carlo approach to uncertainty propagation is considered for PE parameters specified in terms of a probability distribution.In Section Four, a different specification of the PE parameters in terms of an interval range is considered, again using Monte Carlo analysis in order to facilitate comparison with the results of Section Five where an interval approach to the uncertainty propagation is adopted.The paper is completed with some discussion and conclusions.

The computational model
As described in the last section, the situation in this paper is highly idealised.First of all, the structure under consideration is a thin rectangular plate under a uniaxial tension as illustrated in Fig. 1, the stress applied will be denoted σ.The damage under consideration is a central fatigue crack in mode 1 opening, i.e. with its length running perpendicular to the loading axis.The width of the crack is denoted by 2a and the width of the plate by 2w.
An extremely important quantity governing the evolution of the crack is the Stress Intensity Factor (SIF) K which governs the magnification factor for the imposed stress σ in the neighbourhood of the crack tip, this is defined by, Y (ā) is a geometrical factor which is needed in order to accommodate the fact that the plate has finite width -it is a function of the dimensionless parameter ā = a/(2w).In fact, for an infinite plate, an analytical solution yields, One of the reasons for choosing the plate geometry assumed here, is that there are a number of well-validated approximations for the geometry factor in the SIF for a plate of finite width, the reference [8] cites examples, The effect of using these different approximations will be investigated later.
The law of damage propagation appropriate for a fatigue crack in a metallic plate is the Paris-Erdogan (PE) law [8].This relates the crack growth rate da/dn to the crack tip SIF range ∆K as follows, where C and m are constants of the material called respectively, the Paris coefficient and the Paris exponent.(The units of C are MPa −m m m/2+1 and m is dimensionless.)The parameters depend on the microstructure of the material and therefore can vary significantly from sample to sample.In fact, they can vary with position within a single sample.The variation of these coefficients, and its effect on the lifetime of the specimen is the driving force for this paper.It will be assumed that the material of the plate is the Titanium alloy Ti-6Al-4V, because this can be considered as an approximately isotropic metal and because it has been the subject of research to find values of the Paris coefficient and exponent [5,6].(In reality, the α phase grain size does vary with direction; this property is neglected here.) The SIF range ∆K in the situation considered here is generated simply by the fact that a time-varying stress σ(t) will have a corresponding stress range ∆σ, i.e., ∆K = Y (ā).∆σ √ πa (7) In order to completely specify the computation, it is necessary to specify the bounds of the stress range; this is usually done by defining the load ratio R by, Under conditions of cyclic loading the crack length will grow monotonically.One could specify two termination criteria for the computation i.e. the point where the plate fails.The first condition is simply that a = w i.e. the crack has reached the edge of the plate.This is a little unsatisfactory as the approximations Eqs (2) to (5) begin to lose their applicability as a ≈ w.A better termination criterion is to take the point when the crack-tip SIF exceeds the fracture toughness K Ic of the material.K Ic can be regarded as a critical value of the stress intensity beyond which the material cannot withstand fracture.The point at which the crack-tip SIF exceeds the fracture toughness will occur at a critical crack length a c .This value has to be obtained using numerical methods for expressions Eqs (3) to (5).In the case of the infinite-plate approximation Eq. ( 2), one can readily show that, The calculations which follow are all terminated when a = a c and failure is then assumed to occur immediately.
The numerical integration scheme used to step forward Eq. ( 6) in time is also chosen to be as simple as possible in order to maximise the possibility that it will generalise appropriately to the interval arithmetic case.A simple 'rectangle-rule' (more formally a forward-Euler explicit) recursion is used so that, The initial conditions are n 0 = 0 and (usually) a = 10 mm (corresponding to a crack length of 20 mm).In order to ensure accuracy, an appropriately small step-size dn is chosen for each integration on the basis of the size of (da/dn) 0 .During the integration, steps are taken in groups of 10 5 and after each group the step size is multiplied by a factor of 10.Note that da/dn is usually given in mm/cycles, while the SIF is usually given in MPa/ √ m (note that 'm' here denotes metres and not the Paris exponent).This means that a little care is needed in handling the various units.
The dimensions of the simulated plate were chosen to give correspondence with the experiments conducted in [5] in order to measure the Paris parameters.An appropriate choice of dimensions and load range to give a ∆K of 10 MPa/ √ m (the order of the quantity in [5]) is width w of 200 mm, thickness t of 2 mm and ∆σ of 40 MPa.The latter corresponds to a maximum applied load of 16 kN.

Monte Carlo case study 1
For the first case study, the Paris parameters were taken from the experimental quantification in [5].This paper followed an interesting approach whereby the statistical properties of the parameters were estimated from a single material sample i.e. exploiting the within-sample variation of the parameters as a crack progressed through the single sample.The authors of [5] showed that the parameter C was well-described (i.e. at confidence level 97.5%) by a log-normal distribution.Alternatively, log C (and this is a logarithm to base 10), was well-described by a normal distribution with mean −12.3 and standard deviation 4.1.The parameter m was found to be described (with a 90% confidence level) by a normal distribution with mean 7.3 and standard deviation 4.1.To give some idea of the scatter in these quantities, the 95% confidence interval for C was (4.613 × 10 −21 , 5.445 × 10 −5 ) and for m was (−0.736, 15.336).This naively gives a range for da/dn of (4.613 × 10 −21 , 1.933 × 10 12 ) and shows why it is necessary to use care in the selection of the integration step-size.
As a check of the accuracy of the integration scheme, the Eq. ( 6) was integrated from a = 10 mm up to a = w using the mean values of C and m and then compared with the exact result, The result of this calculation is shown in Fig. 2. The calculation was terminated at a critical crack size of 79.5 mm.This was determined from Eq. ( 9) using a fracture toughness value of 100 MPa √ m from [10].The value of K max corresponds to an assumed load ratio R of 0.8.The two curves overlay as desired, the lifetime of the specimen was estimated as 4627800 cycles compared with the exact value of 4627400.
For the Monte Carlo computation, 1000 sample runs of the integration were carried out, with C and m sampled from the distributions described above.(Note that the Monte Carlo method is versatile enough to account for correlations between the parameters C and m, however the distributions given in Ref. [5] were effectively the marginal distributions, so no information about possible correlation was available.The simulation here therefore assumes that the parameters are statistically independent.)Note that as the accuracy of the Monte Carlo approach is dominated by a term of the form 1/ √ S where S is the number of samples [11], for a univariate simulation 100 runs would usually be sufficient to give accuracy to one decimal place and 10000 runs would be needed for accuracy to 2 decimal places; the accuracy here is intermediate between these values, the number of runs being chosen for speed of computation as much as accuracy.1000 runs were carried out for each of the possible SIF forms Eqs (2) to (5).The calculations were carried out in Matlab version 7 using a Windows XP computer with a 3 GHz Pentium 4 processor.When integrating systems with the SIFS specified by Eqs (4) and ( 5) it was necessary to use protected versions of the functions to avoid the possibility of complex roots.For each run, the lifetime of the sample was defined as the number of cycles needed for a to reach a c , the exact lifetime was extracted by interpolation over the last time-step.Figure 3 shows a sample of the crack growth curves from this simulation, it is clear that the curves are evenly spread over several orders of magnitude, giving a huge range for the estimated lifetime.In order to display the results, the Probability Density Function for the lifetime was extracted from each 1000-point data set using Kernel Density Estimation (KDE).This essentially builds up the distribution from a series of kernel distributions of equal spread, each centred on a data point [12].KDE with a Gaussian kernel was used and the  4) Eq. ( 5) Eq. ( 6) Fig. 4. Probability Density Functions for lifetime for each SIF expression (Eqs (3) to ( 6)).First Monte Carlo simulation.
smoothing was chosen by cross-vaildation as described in [12].The results are shown in Fig. 4.
It is clear from Fig. 4, that all the expressions for the SIF, i.e. different damage propagation laws, lead to broadly similar distributions for the lifetimes.One can see that the expression for the infinite plate (Eq.( 2)) leads to slightly higher estimates for the mode of the distribution than the others.This is to be expected as Eqs (3) to (5) specifically account for the higher values of the SIF as the crack-tip approaches the plate edge.The statistics of the Monte Carlo simulation are summarised in Table 1.
The estimates of the lifetime are clearly of no use as the lower bound in each case is of the order of fractions of a second (even assuming very low frequencies of excitation)).This is entirely due to the scatter in the estimated Paris coefficients.While the approach taken in [5] yields an advantage in giving the whole probability distribution for the Paris Parameters, one must call into question the precision of the estimates.Before blaming the estimated parameters, however, one should note that the basic premise of [5] is that the parameters vary as a random field across the sample and it is possibly inappropriate to use the distributions to fix the parameters as constants over a whole sample.This question will be addressed in further work.
Finally, as a matter of interest, the distributions for the lifetimes all look quite like Gaussian distributions.In order to investigate this, the estimated PDF for the data corresponding to Eq. ( 2) is plotted in Fig. 5 and compared with a Gaussian distribution with the same mean and standard deviation as the 1000-point Monte Carlo sample.
There is certaintly a broad similarity between the estimated PDF and the Gaussian.However, as the distribution obtained is not actually very useful for the purposes of lifetime prediction, the similarity will not be pursued here.

Monte Carlo case study 2
The second Monte Carlo simulation is based on the Paris parameters tabulated in [6].The first point to note is that this paper considers a different version of the alloy Ti-6Al-4V in that the authors use Ti-6Al-4V-0.1 Ru (extra-low interstitials, (ELI)).As opposed to reference [5], the authors only extract one set of Paris parameters from each sample, of which 16 are considered.The tests are carried out in air and in both aerated and deaerated sea water, but still show considerably less scatter in the Paris coefficients than those from [5], so the data from all 16 tests are included here.The Paris parameters in [6] are tabulated and there is no suggestion for an underlying probability distribution; the parameter uncertainty will therefore be considered here as specified by an interval.The interval corresponding to log C is (−15.0,−11.6) and the interval for m is (3.7, 6.2).As the uncertainty is specified in terms of an interval, it is appropriate to propagate the uncertainty through the lifetime calculation using interval arithmetic.However, to explore whether the interval calculation leads to tight bounds on the lifetime, it will be compared with a Monte Carlo calculation carried out as in the previous section.In order to carry out the Monte Carlo runs, a uniform distribution was assumed for each parameter covering the appropriate interval.For the sake of simplicity, the parameters were once again assumed to be statistically independent, although it would have been possible here to estimate the statistical correlation from the dataset.(Although possible, using the dataset might not have been appropriate as the 16 samples came from three different environments.)The geometry for the plate and the loading regime is exactly as in the previous section.Apart from the different Paris parameters, the only difference in the calculation is in the value of the fracture toughness, which is taken as 75 MPa √ m as this is the (minimum) value quoted in [6].For each expression of the SIF (Eqs (2) to ( 5)), 1000 samples of the crack growth curve were obtained leading to 1000 values of the lifetime.Figure 6 shows the 1000 crack growth curves corresponding to the infinite plate SIF in Eq. ( 2).
The lifetime PDFs for the various SIF expressions are given in Fig. 7.
The PDFs in Fig. 7 are all very consistent with each other; all are clearly non-Gaussian.The other point to note is that the distribution corresponding to Eq. ( 5) appears to be bimodal.This is probably worth investigating further.
The statistics of the second Monte Carlo simulation are given in Table 2.
As one might expect, the scatter in the lifetime estimates from this set of Paris parameters is considerably smaller and actually yields usable estimates for the specimen lifetime.For example, If one uses the minimum lifetime corresponding to Eq. ( 2), one obtains 11817000 cycles for the crack to propagate from the initial 20 mm to the point where the specimen fails.At a nominal fatigue test frequency of 20 Hz for example, this corresponds to approximately 164 hours.In practice, this would be certainly enough time to take corrective action.The scatter on the lifetimes is still about 5 orders of magnitude.2) to ( 5)).Second Monte Carlo simulation.

Interval analysis
As stated previously, the Monte Carlo technique of propagating many crisp (or certain) input data through the integration formula with crisp parameters in order to assess the range of output values which can be returned by equations with uncertain parameters becomes unfeasible very quickly as the number of uncertain parameters increases.Fortunately, there exist some alternatives to the above approach which employ dedicated arithmetics to allow for an assessment (generally guaranteed to contain the true solution set) to be provided with a single run.The first, and most well known, of these methods uses interval arithmetic, introduced in [7], to consider the outcome of problems where the parameters are represented by ranges.
The basic principal of interval arithmetic is that it represents an interval number by its upper and lower bound values.There exist rules for the standard mathematical operations on interval numbers which are well-known, but are repeated here for completeness.If an interval number takes the form [low, high] then addition of two interval numbers is: Multiplication of two interval quantities is given by: which takes into account the effect of negative values.Thirdly, multiplication of an interval number by a positive scalar c is given by: Finally, if an interval number is the input to a monotonically increasing function F then this can be written as: The biggest drawback with interval arithmetic is the problem of dependency, which occurs if variables occur more than once within the arithmetic of some problem.Interval arithmetic treats all ranges as being entirely independent of each other, which, if this should not be the case, may result in overestimation of the true solution set.This problem of 'bound explosion' is potentially most serious in recursive algorithms and as the integration rule ( 11) is recursive, it was a concern here.The approach to interval arithmetic taken here is as simple as possible; the quantities in (11) are interpreted as interval quantities (except dn, which is allowed to stay crisp).The recursion is then implemented in interval arithmetic starting from a thin interval [a 0 , a 0 ] with a 0 = 10 mm.The calculation was carried out using the package INTLAB.The calculation is very straightforward except for one point.Because the PE equation raises the SIF range to a positive power m, the upper limit of a reaches the critical crack size considerably before the lower limit.This means that if the calculation is continued up to the point where the lower limit reaches the critical crack size, the upper limit overflows in the computer's floating point arithmetic.Experimenting with the step size for the integration did not give a solution to this problem.The approach adopted here was to 'freeze' the upper limit of the interval a at a c and continue the calculation.In general, this may not give the correct final interval bounds; however it does ensure that the lower bound on the lifetime, which is the quantity of greatest interest, is correct.The interval calculation was only carried out here using the infinite-plate approximation to the SIF (Eq.( 2)) and the Paris parameters from Section four.The resulting crack growth curves are given in Fig. 8.
The bounds on the (Log) lifetime provided by the interval calculation are given by (6.9858, 12.7817).The bounds from the corresponding Monte Carlo calculation were (7.0725, 12.4188).It is clear that the Monte Carlo bounds are within the interval bounds, which means that the estimated lower bound on the lifetime is conservative as it should be.Interpreting the lower bound in terms of a 20 Hz fatigue test as before, the calculation gives a warning that failure will occur in 134 hours -30 hours ahead of the 164 hours predicted by Monte Carlo.Another advantage of adopting the interval approach is that the calculation took approximately 50 seconds as opposed to the Monte Carlo computations which took over an hour on the PC used.
One might imagine that the uncertainty in the lifetime estimate will be less if the crack is discovered at a later time and is therefore longer.This is not the case as a result of the exponential form of the PE law. Figure 9 shows the interval estimates of the crack growth laws for initial crack lengths of 10 mm, 20 mm, 30 mm and 40 mm.The uncertainty still spans several orders of magnitude even though the amount of warning before failure falls to 117140 cycles for the crack which starts from 40 mm.(In the context of the 20 Hz fatigue test, this amounts to a warning of 1.63 hours.) Under certain circumstances, it is possible to simplify matters even further from interval analysis.In general, the possible values for the parameters in an interval calculation occupy the interior of a hyprercuboid.The vertices of  this hypercuboid are given by the various combinations of the lower and upper bounds of each interval parameter.It can be shown that if the propagation function is monotonic in the various interval parameters, then the bounds of the output set (output ranges of the propagation function) can be found by propagating the vertices of the parameter set and taking the extremes [13].As the PE equation is clearly monotonic in the parameters Cand m, it is worth investigating if vertex propagation gives good results in the current problem.As there are only two parameters, there are only four possible vertex combinations.When these are used to carry out the integrations for a, the results are as shown in Fig. 10.
The upper and lower bounds on the (Log) lifetime are (6.9848,12.7807).Compared with the results, (6.9858, 12.7817), of the interval calculation, this shows that (at least up to the accuracy of the numerical integration) the interval bounds are tight.This is very important as it shows that the Monte Carlo calculation has overestimated the minimum lifetime of the specimen.A larger Monte Carlo calculation would be needed in order to refine the lifetime estimate.

Discussion and conclusions
The main conclusions of this paper are that, in a very idealised context, there are marked advantages in using an interval approach to the propagation of uncertainty through a lifetime estimate.The most important advantage is the conservative nature of the lifetime estimate.The Monte Carlo calculation suggested that the specimen could last another 164 hours before failure, whereas in fact the interval and vertex calculations showed that the lifetime could be as low as 134 hours, failure occurring over a day earlier.The reason for the overestimate of the Monte Carlo approach is the comparatively low number of samples used.The Monte Carlo results could be improved by increasing the number of samples.However, even with 1000 samples, the execution time for the Monte Carlo estimate was over an hour compared with 50 seconds for the interval calculation (and an even shorter time for the vertex calculation).This problem with the Monte Carlo approach would only become worse with more parameters to consider.The arguments are not entirely in favour of interval arithmetic; in this case the interval bounds were tight, largely because of the simplicity of the PE law and the 'rectangle rule' integration.In general, with arbitrarily complicated damage propagation laws and with more complex integration algorithms, the tightness of the bounds is not guaranteed and the calculation may even suffer from bound explosion.The same reasons why the interval calculation gave tight bounds, made the vertex propagation possible.In general this may not be the case.There does exist an alternative to interval arithmetic, termed affine arithmetic [14], which, in many situations, significantly reduces the bounds explosion problem through its ability to monitor dependencies between uncertain parameters through a calculation.Another advantage of affine arithmetic is the inclusion of an approximation error term which can be readily converted into upper and lower estimates of the conservatism within the calculation, something which is not possible with interval arithmetic.
It is worth stating exactly how idealised the situation in this paper is, lest the reader begin to think that damage prognosis is straightforward.There are a number of very important restrictions here: -The material under investigation here is an isotropic metal.In general, it will be necessary to make lifetime estimates for anisotropic composites.It may even be necessary to take into account the anisotropy of some metals.For the material considered here, estimates of the Paris law parameters were available.In one case (as discussed in Section Three) these estimates were subject to so much scatter that the lifetime estimates from them were unusable.Even if experimentally obtained Paris law coefficients are available, one must be very careful in considering the nature of the alloy in question.In this paper, the two variants of Ti-6Al-4V had markedly different parameters.For example, the alloy considered in Section Three had a mean value for mof 7.3 and this value was not even included in the range of values of m (3.7-6.2) for the alloy in Section Four.-The structure is simple; a rectangular plate with a central crack.This means that experimentally validated expressions for the SIFs needed for the PE law were available.This will not be true in general.A typical aircraft component may have many geometrical features like ribs and stringers which complicate matters.-The damage type (a fatigue crack) was the simplest imaginable and a validated damage propagation law -the Paris-Erdogan law -was available.In general, the damage law may well be complex and very time-consuming to evaluate.It may be the case that the propagation of the damage requires a large multi-scale physics calculation.Also, in composites for example, which have many different damage scenarios, a damage propagation law will be needed for each type.-The loading regime is very simple here -the load is uniaxial and harmonic.In general damage will evolve under multi-axial loading.This means that even in the context of a fatigue crack, it will experience multi-mode opening.Another vital consideration here is that the future loading was known exactly in this case as a harmonic function -a more realistic situation will demand that the algorithm can deal with uncertain broad band loading.
Given these caveats, it is worth stressing that prognosis has at least proved possible in the idealised case considered here.

Fig. 5 .
Fig.5.Comparison of PDF of lifetimes using simulation from Eq. (2) and a fitted Gaussian.

Table 1
Statistics of first Monte Carlo simulation

Table 2
Statistics of second Monte Carlo simulation