Deriving exponential and power law slow crack growth parameters from dynamic fatigue of indented or as-polished window material coupons

Abstract. Knowledge of slow crack growth rate in a window material is required to determine the design life of a window in which cracks grow under the influence of applied service stress. Parameters that characterize crack growth rate can be measured by dynamic fatigue in which coupons are taken to failure at several constant stress rates. The slower the stress rate, the more time is available for cracks to grow during the strength measurement, and the lower the stress at which the coupon fails. That is, coupons tested at a slower stress rate are weaker than coupons tested at a faster stress rate. The rate of slow crack growth is commonly fit to either a power law or an exponential law. The exponential law provides the more conservative estimate of design life and is prescribed for use by the U.S. National Aeronautics and Space Administration to design windows for manned vehicles. Well-established analytical equations are used to derive the power law crack growth parameters from dynamic fatigue measurements of as-polished (unindented) coupons. There are no exact analytical equations to derive power law parameters from indented coupons or to derive exponential law parameters from unindented or indented coupons. Approximate procedures have been used to derive crack growth parameters when there are no exact equations. We now describe a numerical method that gives the power and exponential laws for both unindented and indented coupons by least-squares fitting of dynamic fatigue measurements. A MATLAB code is provided to carry out the calculations.


Introduction
The lifetime of a brittle optical window subject to stress is limited by the rate at which microscopic cracks grow under the influence of stress until the window fractures. For many window materials, slow crack growth is promoted by atmospheric moisture, which participates in bond breaking chemical reactions at the tip of a growing crack. 1,2 It is necessary to know the upper limit for the rate of slow crack growth to establish the design life of a window.
Slow crack growth rate is commonly characterized by dynamic fatigue experiments in which test coupons are broken at a series of constant applied stress rates:σ a ¼ dσ a ∕dt, where σ a is applied stress, t is time, andσ a is the applied stress rate. The more rapidly stress is applied, the less time there is for cracks to grow and the higher the stress at which failure occurs. Coupons are tested either with their as-polished surface having a natural distribution of microscopic flaws, or they are tested with a deliberate indentation that creates a reproducible flaw size larger than natural flaws.
Dynamic fatigue experimental results are fit equally well in the measured range by two empirical crack growth rate laws: a power law and an exponential law. When extrapolated to low stress intensity factors (K I ) common to many window applications, the two laws diverge and can predict window lifetimes that differ by several orders of magnitude. For most materials, we do not have data at low K I to differentiate which crack growth law applies. In the absence of this knowledge, U.S. National Aeronautics and Space Administration design rules prescribe that the more conservative exponential law shall be used for window lifetime prediction. 3 Power law slow crack growth parameters can be derived from dynamic fatigue measurements of as-polished coupons by using analytical equations given in ASTM C1368. 4 Fuller et al. 5 gave an approximate analytical fit to a numerical solution to derive power law parameters from dynamic fatigue of indented coupons. There are no analytical equations to derive exponential law parameters from dynamic fatigue of as-polished or indented coupons. Jakus et al. 6 described a numerical method to obtain the exponential law from dynamic fatigue of as-polished (unindented) coupons. Choi et al. 7 presented an approximate method to obtain the exponential law from dynamic fatigue of as-polished (unindented) coupons. Ritter et al. 8 gave an approximate method to obtain the exponential law from dynamic fatigue of indented coupons. The purpose of our paper is to present a numerical method for least-squares fitting of all data points in dynamic fatigue of as-polished or indented coupons to find parameters for the exponential law and the power law.  Tables 1 and 2 for as-polished BK7 glass disks immersed in water and broken in biaxial flexure at stress rates ofσ a ¼ 0.39, 7.7, 77, or 620 MPa∕s. 5 Figure 1(b) shows results from Tables 3 and 4 for BK7 glass disks indented by a Vickers indenter prior to measuring the dynamic fatigue strength in water.

Dynamic Fatigue Measurements for BK7 Glass
High scatter in the strength of as-polished coupons in Fig. 1(a) is typical for optically polished glass and ceramics. The standard deviation of the mean for the slope is 25% and 1.4% for the intercept. The intercept is the ordinate when logðstress rateÞ ¼ 0. Scatter is reduced by a factor of 10 for indented coupons in Fig. 1(b) because failure originates at the indentation crack, which is larger than natural flaws and has a reproducible size. Fewer indented coupons are required for testing because of the low scatter. The natural flaws in the as-polished specimens were identified by Quinn 5 as scratch/dig flaws, invisible on the polished surfaces, created during the specimen preparation that had all surface traces eliminated by final polishing. The near equivalence of the slopes in Figs. 1(a) and 1(b) is fortuitous because of the 25% uncertainty in the slope of Fig. 1(a). The theoretical slope in Fig. 1(a) should be 1∕ðn þ 1Þ, where n is the exponent in the power law Eq. (2) in the next section. The theoretical slope in Fig. 1(b) should be 1∕ðn 0 0 þ 1Þ, where the relation between n and n 0 0 is given by Eq. (15) in the paper of Fuller et al. 5 n ¼ ð4n 0 0 − 2Þ∕3. From the observed slope in Fig. 1(b) = 0.06101, we find n 0 0 ¼ 15.391. From this value of n 0 0 , we calculate n ¼ 19.855 and the slope of Fig. 1(a) should be 1∕ðn þ 1Þ ¼ 0.04795. The observed slope in Fig. 1(a) is 0.06212 AE 25% ¼ 0.0466 to 0.0776, which includes the theoretical value 0.04795. The only reason why the slopes in Figs. 1(a) and 1(b) are nearly the same is because of the large scatter in the data in Fig. 1(a).

Power Law and Exponential Law for Crack Growth Rate
Two types of empirical equations are commonly used to fit observed slow crack growth rates (v). The exponential law equation is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 9 7 Exponential law∶ v ¼  where c is flaw size, t is time, v o and β are measured parameters, and K I is the stress intensity factor. The elastic stress field magnitude near the crack tip is proportional to K I . The power law slow crack growth rate is given by the equation where A Ã and n are measured parameters and K Ic is the critical stress intensity factor (fracture toughness) for catastrophic failure. Equations (1) and (2) generally fit observed failure strengths equally well in dynamic fatigue experiments in the measured range of stress rates. However, Fig. 2 shows that the two crack growth laws diverge when extrapolated to low stress intensity factor (K I ). The exponential law predicts higher crack growth rate at low stress intensity factor than the rate predicted by the power law.
The low K I region is where windows with low applied stress spend most of their lives. As cracks grow under the influence of a constant stress, crack size increases and the rate of crack  growth accelerates, eventually becoming catastrophic. Because most of a window's life is spent at low K I , window lifetime predicted by the power law can be up to several orders of magnitude longer than lifetime predicted by the exponential law.
For most materials, we do not have data at low K I to know which crack growth law applies. Static fatigue experiments at a constant low stress until failure occurs can, in principle, distinguish the two growth laws, but such experiments are slow to conduct. Median times to failure for unindented material in Table 1 range from 0.2 to 210 s at decreasing stress rates. For indented material, median times range from 0.8 to 4600 s. A reviewer of this paper believes that static fatigue testing is more realistic than dynamic fatigue because trends can change at very long lifetimes. That point is well taken, but static fatigue measurements took that reviewer up to 2 years to conduct. Dynamic fatigue is a compromise between obtaining meaningful measurements and concluding the work in days instead of years. In the absence of knowing the applicable crack growth law, U.S. National Aeronautics and Space Administration design rules prescribe that the more conservative exponential law Eq. (1) shall be used for window lifetime prediction. 3 Power law crack growth parameters A Ã and n for Eq. (2) for unindented coupons are derived from the slope and intercept of the dynamic fatigue graph in Fig. 1(a) by analytical equations given in ASTM C1368. 4 There are no analytical equations to derive exponential law parameters from dynamic fatigue of unindented or indented coupons because Eq. (1) must be integrated numerically. We now present a numerical method for least-squares fitting of all data points in dynamic fatigue of as-polished (unindented) or indented coupons in Fig. 1(a) or Fig. 1(b) to find parameters for the exponential law or the power law.

Stress Intensity Field for Indented Test Coupons
First, we introduce equations taken from Fuller, Lawn, Cook, and coworkers 5,[9][10][11][12][13] to analyze dynamic fatigue of indented coupons. Then we apply these equations to find the parameters v o and β for the exponential law Eq. (1) or A Ã and n for the power law Eq. (2) for indented or unindented coupons.
Indentation of a ceramic or glass surface with a Vickers pyramidal diamond indenter leaves a residual stress field due to an elastic/plastic mismatch resulting in a localized stress intensity factor K r that depends on material and indentation load and the equation is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 1 2 5 Residual stress intensity factor∶ K r ¼ χ r P c 3∕2 : Here, P is the indentation force, c is the radius ("crack size") of the half-penny radial or median crack produced by indentation, and χ r is a proportionality constant that depends on the Fig. 2 Comparison of exponential and power law behavior as a function of K I . Crack growth rate parameters are derived from coupons in Fig. 1 elastic modulus and hardness of the material. If stress σ a is applied to an indented test coupon in a dynamic fatigue experiment, the combined stress intensity factor is the sum of the residual stress intensity factor in Eq. (3) plus the applied stress intensity factor, which is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 9 9 Applied stress intensity factor∶ K a ¼ Yσ a ffiffi ffi c p ; (4) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 5 6 Combined stress intensity factor∶ In Eqs. (4) and (5), Y is a geometric factor for which we choose the value Y ¼ 2∕ ffiffiffi π p for a halfpenny (semicircular) crack.
When an indentation crack of initial size c i is first formed, as in an equilibrium condition, the residual stress intensity factor is equal to the critical stress intensity factor, i.e., K r ¼ K Ic . Substituting c i into the denominator of Eq. (3) and K Ic in place of K r gives the relation ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 5 6 0 Combined stress intensity factor∶ The inert strength of a material is the average strength measured when a coupon is taken to failure in an inert atmosphere (such as dry nitrogen) in which slow crack growth is negligible. We can find an expression for inert strength by replacing K I with K Ic at the left side of Eq. (5) and then solving for the applied stress needed to break the coupon: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 4 6 6 Stress at failure∶ The applied stress σ a in Eq. (7) passes through a maximum that we will call σ m when the crack size is c ¼ c m . The maximum stress σ m is the inert strength, which is the strength of indented coupons when there is no slow crack growth. To find an expression for σ m , set the derivative dσ a ∕dc for Eq. (7) equal to 0 and solve for the maximum stress σ m , and the equation is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 3 6 0 ffiffiffiffiffi ffi c m p : Substitute the expression for σ m from Eq. (8) back into Eq. (7) to find an expression for c m , which is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 2 9 0 Crack size at maximum stress∶ c m ¼ 4χ r P K Ic At failure in an inert atmosphere, the stress intensity factor is K I ¼ K Ic , the stress is σ m , and the crack size is c ¼ c m . Making these substitutions into Eq. (6) gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 2 2 1 From Eq.
which we can substitute for ffiffiffiffiffi ffi c m p in the last term in Eq. (10) to solve for the initial indentation flaw size c i , giving the equation In summary, we have the following information from dynamic fatigue of indented coupons: 1. The inert strength σ m can be equated to the average strength of indented coupons measured in dry nitrogen at a high stress rate. 2. Crack size c m for the inert strength is calculated from σ m with Eq. (11). 3. The initial crack size c i from indentation is calculated from c m with Eq. (12). 4. The dynamic fatigue experiment in Fig. 1(b) gives multiple measurements of failure stress σ f for indented coupons tested at different constant stress ratesσ a ¼ dσ a ∕dt.

Evolution of K I and Crack Size in Dynamic Fatigue
Modeling the evolution of crack size and mechanical strength during dynamic fatigue is needed to find A Ã and n for the power law or v o and β for the exponential law. The power law has an analytical solution, but the exponential law requires numerical integration. We developed a numerical method which can be modified for use with either as-polished (unindented) or indented coupons, and for either power or exponential law crack growth. The procedure to be described below allows us to compute curves in Fig. 3 showing the progression of stress intensity factor and crack size during dynamic fatigue. The abscissa is the fraction of time to failure. Actual times to failure decrease by four orders of magnitude as the stress rate increases by four orders of magnitude. In Fig. 3(a) the stress intensity factor decreases rapidly from its initial value K I ¼ K Ic down to a nearly constant level that persists through most of the dynamic fatigue measurement, and then climbs rapidly to failure. Crack size in Fig. 3(b) jumps from its initial value c i and then increases gradually until it rises rapidly when close to failure. It takes 90% of the measurement time for the crack to double in size from the almost level value of ∼20 μm early in the experiment up to ∼40 μm near the end of the experiment.  14 Initial stress intensity factor: K I ¼ 0,

Select time step Δt
A variable time step is used to provide sufficient modeling fidelity when crack velocity is high, but reducing processing time at low crack velocity. For a time step beginning at time t, the step size is chosen as where v t is the crack velocity at time t andσ a is the stress rate in MPa/s. Equation (14) advances the crack by 10 −8 m in one time step. However, if Δt t exceeds Δt max , then use Δt max , which limits the change in stress in one step to 0.01 MPa.

New crack size at time t
The new crack size is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 2 9 0 where c t−1 is the crack size in the previous step, v t−1 is the crack velocity from the previous step, and Δt t−1 is the time step selected from the previous step.

New stress intensity factor at time t
The new value of K It at time t is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 1 8 6 Indented coupon∶ þ Yσ a t ffiffiffiffi c t p ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 1 2 9 Unindented coupon∶ K It ¼ Yσ a t ffiffiffiffi c t p ; whereσ a is the stress rate and t is the elapsed time since the start of the experiment.

New crack velocity at time t
The new crack velocity is v ¼ v o e βK It for the exponential law or v ¼ A Ã ð K It K Ic Þ n for the power law.

Calculate failure stress
Repeat Secs. 5.2.2 to 5.2.5 to obtain a new time increment, new crack size, new stress intensity factor, and new crack velocity at each new time until K It ≥ K Ic . The failure stress σ f is the highest value of σ a before the stress intensity factor exceeds K Ic .

Calculate failure stress for every measured point
Each data point in a dynamic fatigue experiment could have its own stress rateσ a . Generally, the experiment is designed so that there are groups of coupons in which each member of the group has the same stress rate. For each stress rate, calculate the predicted failure stress for a given pair of crack growth parameters (v o and β or A Ã and n) for each experimental data point.

Find best crack growth parameters
For a given trial set of crack growth parameters (v o and β or A Ã and n) compute R, which is the sum of squares of differences of logarithms between the observed failure stress and predicted failure stress for every point and the related equation is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 6 ; 4 6 9 R ¼ X ðlogðσ calc Þ − logðσ obs ÞÞ 2 : The best crack growth parameters for each crack growth law are the pair that minimizes R. For indented BK7 data in Tables 1 and 2, the minimum R value is 0.0049 at v o ¼ 9.57 Ã 10 −15 m∕s and β ¼ 42.39. We observed that the minimum value of R lies in a shallow well. Optimum slow crack growth parameters are sensitive to the number of significant digits used for stress rate and failure stress. We found it necessary to retain four significant figures for stress rate and failure stress as input to the computation of crack growth parameters.

Parameter Search Strategy
Any method to search for pairs of crack growth parameters to minimize R in Eq. (19) can be used. We employed a constrained interior-point approach utilizing barrier functions in the MATLAB (Mathworks, Natick, Massachusetts) lsqnonlin optimization solver embedded in the function solve(). Solve() analyzes the problem structure to select specific algorithms to use. This nonlinear constrained optimization problem should use lsqnonlin, which is confirmed by reading the output from solve(). R in Eq. (19) is evaluated at an initial point in parameter space and at a few surrounding points. Then, either a direct step is taken based on a finite difference gradient in the parameter space or, failing that due to local non-convexity, a conjugate gradient step within a trust region is taken. This process is repeated to descend down the gradient of the objective function R until a local minimum is reached. Parameter search space is constrained within upper and lower limits specified by user input. We demonstrated that the global minimum was obtained by comparison of the lsqnonlin minimum with that found by another global optimization algorithm (pattern search). [17][18][19][20]

Finding Exponential Crack Growth Law for Dynamic Fatigue of Indented BK7 Glass Disks
For indented dynamic fatigue of BK7 glass disks in Tables 1 and 2, the procedure in Sec. 5.2 gives the following intermediate and final results. New stress intensity factor K It ¼ K Ic

Method Validation
To validate our results for unindented (as-polished) coupons, Table 5 compares the power law derived by our numerical method with the power law derived from analytical equations of ASTM C1368. Parameters derived by both methods for three different materials are almost identical.
To validate our method for indented coupons, we can compare our power law parameters with those from the approximate equations given by Fuller et al. 5,11 Table 6 shows good agreement between the two methods for three different glass materials. Raw data for indented BK7 glass and soda-lime glass have relatively low scatter (slope uncertainties of 2% and 6%). Higher scatter of 18% in the raw data for the specialty glass did not lead to much discrepancy in the parameters derived by our method and Fuller's method. Fuller et al. 11 state that their analytical equations are within 1% of those from numerical integration for the exponent n and within 10% for the multiplier (A Ã ). We observe the same level of agreement (1% and 10%) between our numerical integration and Fuller's equations.
We also verified that our numerical method is predicting close to the same failure stress calculated by the equations of ASTM C1368 for unindented coupons or by Fuller's approximate equations for indented coupons over the wide range of input parameters in Table 7. Over the  Fig. 1(a) parameter space with 6 million combinations in Table 7, 99% of results are within 0.2% of the ASTM method for unindented coupons, and within 0.84% of Fuller's approximate method for indented coupons. Tables 5 and 6 demonstrate that our numerical method finds the same global minimum expected from the ASTM or Fuller equations.

Uncertainty in Crack Growth Parameters
To estimate uncertainties in Eq. (1) crack growth parameters β and log v o , we employ the jackknife procedure. 21,22 In this procedure, we delete the first point (σ a ¼ 77. 17 MPa∕s, Table 3 and compute the best values of β and log v o from the remaining 33 points. The first row of Table 8 shows that the best values are β ¼ 42.74 and log v o ¼ −14.08. For the second row of Table 8, we replace the first point back into the set and delete the second point. We continue deleting one point at a time from the original data to generate 34 new values of β and log v o in columns 2 and 3 of Table 8. The standard error (u ¼ standard deviation of the mean) of a parameter is an estimate of the uncertainty of that parameter. At the bottom of Table 8 Fig. 1  whereβ is the mean value of β, and n is the total number of data points. We write analogous expressions for the standard error and variance of að¼ log v o Þ. The second term on the right in Eq. (20) is the population standard deviation (STDEV.P in Excel ® ).
Covariance is a measure of the magnitude and direction in which a change in one variable affects another variable. Columns 4 and 5 of Table 8 compute the jackknife estimates of the covariance of β and log v o as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 6 ; 6 6 3 Covariance ðu aβ Þ ¼ ðn − 1Þ Let us estimate the uncertainty in crack velocity u v from the exponential crack growth law v ¼ v o e βK I for best fit parameters β ¼ 42.39 and log v o ¼ −14.02. For a chosen stress intensity factor K I ¼ 0.40 MPa ffiffiffiffi m p , the crack velocity is v ¼ v o e βK I ¼ 2.238 × 10 −7 m∕s. To find the uncertainty in crack velocity, we will use standard errors at the bottom of Table 8: u a ¼ 0.2374, u β ¼ 1.1264, and covariance u aβ ¼ −0.2644. The uncertainty in crack velocity u v is found from the equation Combining the velocity with its uncertainty, we can say v ¼ ð2.238 AE 0.274Þ × 10 −7 ¼ ð1.96 to 2.51Þ × 10 −7 m∕s. In this example, the covariance term in Eq. (23) is almost as large as the sum of the two variance terms, but opposite in sign. If we had neglected covariance, the relative uncertainty in velocity would have been 71% instead of 12%, which emphasizes the importance of including covariance.
To assess the jackknife procedure applied to finding standard errors u slope , u intercept , and covariance u slope;intercept for fitting a straight line, we calculated the standard errors and covariance for the slope and intercept of the straight line in Fig. 1(b) using analytical equations 11, 13, and X1.11 in ASTM C1368. 4 Jackknife standard errors differed from those of the analytical equations by −1%, þ6%, and −27% for u a , u β , and u aβ , respectively.

Discussion and Conclusion
Equations (5)-(12) from indentation fracture mechanics 5,[9][10][11][12][13] predict that the initial size of the indentation flaw is c i ¼ ð0.3968Þc m , where c m is the crack size at failure in an inert atmosphere. Our calculations are predicated on this calculated value of c i . It was observed at NIST that the size of the unstable crack grows in air before the dynamic fatigue experiment begins. It is also assumed that residual stress from indentation continues to exist throughout the mechanical test, but diminishes according to Eq. (3) as crack size grows during the mechanical test. We could question what is the "right" initial crack size to use in simulating dynamic fatigue and does the residual stress relax by mechanisms other than just crack extension during the test? Notice in Fig. 3 that the stress intensity factor drops and the crack size increases in a fraction of a second to values that remain almost flat for much of the dynamic fatigue experiment. We find from simulations in Table 9 that the starting value of crack size [c t−1 in Eq. (16) for the first step of the simulation] of dynamic fatigue can be varied over a factor of 2.5 with little effect on the derived values of the slow crack growth parameters. We still use the same theoretical value of the indentation crack size c i ¼ ð0.3968Þc m , but we allow the crack to propagate after indentation and prior to testing in these simulations. For the data in Tables 3 and 4, we conclude that if the crack grows to less than 2.5 Ã c i between indentation and testing, then pre-test crack growth has little effect on the derived slow crack growth parameters.
The other question raised above is does residual stress relax over time by mechanisms other than crack extension? We cannot provide a universal answer to this question, but we can cite one study 23 of soda-lime glass in which there was no significant decline in residual stress in air or water even up to 40 days. Spontaneous, post-indentation crack growth was observed for more than a day, but the residual stress driving the crack growth appeared to be constant.
In conclusion, we presented a step-by-step recipe for numerical integration of the slow crack growth exponential rate law Eq. (1) and power rate law Eq. (2) to derive the least-squares best-fit parameters v o and β or A Ã and n from dynamic fatigue experiments such as those in Figs. 1(a) and 1(b) with as-polished or indented test coupons. A link to our MATLAB code is provided. This code has been optimized to handle data sets with hundreds of points in a period of minutes. Tables 5 and 6 validate the method by demonstrating that the code gives the same crack growth parameters as found with the analytical equations of ASTM C1368 or close to the approximate equations of Fuller et al. 5,11 Section 7 describes how to use the jackknife method to estimate the variance and covariance of the crack growth parameters and to estimate the uncertainty in crack growth rate from those parameters. Table 9 demonstrates that limited crack growth up to 2.5 times the initial crack size after indentation and before dynamic fatigue testing has little effect on the derived crack growth parameters.