Fatigue sensitivity analysis using complex variable methods

International Journal of Fatigue 40 (2012) 61–73

Contents lists available at SciVerse ScienceDirect

International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue

Fatigue sensitivity analysis using complex variable methods A. Voorhees a,1, H. Millwater a,⇑, R. Bagley a,1, P. Golden b,2 a b

The University of Texas at San Antonio, San Antonio, TX 78249, United States Air Force Research Lab, RXLMN, Wright Patterson Air Force Base, Dayton, OH 45433, United States

a r t i c l e

i n f o

Article history: Received 14 August 2011 Received in revised form 30 December 2011 Accepted 1 January 2012 Available online 24 January 2012 Keywords: Fatigue life prediction Sensitivity methods Complex Taylor series expansion

a b s t r a c t The sensitivity of the computed cycles-to-failure and other lifing estimates to the various input parameters is a valuable, yet largely unexploited, aspect of a fatigue lifing analysis. Two complex variable sensitivity methods, complex Taylor series expansion (CTSE) and Fourier differentiation (FD), are adapted and applied to fatigue analysis through the development of a complex variable fatigue analysis software (CVGROW). The software computes the cycles-to-failure and the sensitivities of the computed cyclesto-failure with respect to parameters of interest such as the initial crack size, material properties, geometry, and loading. The complex variable methods are shown to have advantages over traditional numerical differentiation in that more accurate and stable first and second order derivatives are obtained using CTSE and more accurate and stable higher order derivatives are obtained using FD. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The fatigue of structural components due to repeated loading cycles is a detrimental and dangerous problem. Structural failure due to fatigue can lead to costly and time consuming repairs, early retirement or catastrophic failure of structural and mechanical systems. Therefore, the calculation of the estimated ‘‘life’’ of a component or a structure is a critical element in a fracture control plan. As a result, there are a number of in-house and commercial computer codes such as AFGROW [1], NASGRO [2], and DARWIN [3] for fatigue crack growth evaluation. These codes compute an estimate of the cycles-to-failure of a cracked geometry under load by integrating the crack growth rate equation until failure occurs, with failure typically defined as net-section-yield or unstable fracture. Other ancillary outputs such as crack size and the stress intensity factor as a function of loading cycles, the residual strength, and the critical crack size are also provided. The inputs that are traditionally considered are the initial crack size, the applied loading, the geometry, and the material properties, and, in some cases, inspection schedules and probability-of-detection curves. The sensitivity of the computed cycles-to-failure to the various input parameters is a valuable, yet largely unexploited, aspect of a fatigue analysis. The fatigue analysis often contains: (a) a

⇑ Corresponding author. Address: Department of Mechanical Engineering, 1 UTSA Circle, San Antonio, TX 78249, United States. Tel.: +1 210 458 4481; fax: +1 210 458 6504. E-mail address: [email protected] (H. Millwater). 1 Department of Mechanical Engineering, 1 UTSA Circle, San Antonio, TX 78249, United States 2 Metals, Ceramics & NDE Division, 2230 Tenth St, Ste 1, United States 0142-1123/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2012.01.016

significant uncertainty or variability in the input parameters; (b) simplifications in the modeling such as a simplified geometry and loading; and (c) numerical approximations, e.g., when computing the geometry correction factor using handbook solutions or curve fits to numerical results such as weight functions. Sensitivity analysis can convey the significance of the various inputs on the computed cycles-to-failure. For example, the sensitivities express the expected change in the computed cycles-to-failure given a small change in the initial crack size, the part geometry, the loading cycles (gust, maneuver, ground-air-ground), the material properties (crack growth rate, fracture toughness), etc. If a sensitivity is relatively large, then a more thorough data collection effort or thorough analysis may be warranted. As a result, methods that determine the sensitivity of the cycles-to-failure to the inputs are valuable additions to a fatigue analysis. Both deterministic and probabilistic approaches can be applied for sensitivity analysis. Deterministic studies to determine the partial derivative of the cycles-to-failure with respect to each input of interest can be obtained using numerical finite differencing. McGinty discussed the sensitivities of damage tolerance analysis (DTA) elements by examining the deterministic ratio of the percent change of output to the percent change of input, i.e., the nondimensionalized derivative of the governing equation [4]. The following studies were developed: sensitivity of the stress intensity factor to the stress intensity geometry correction factor (b), the applied stress, and the crack size; sensitivity of the crack growth rate to b, Paris coefficient, Paris exponent, and crack size; sensitivity of the fatigue life to b, Paris coefficient, Paris exponent, and initial and final crack sizes. The results indicated that the geometry correction factor (b) and the applied stress are the important factors with respect to fatigue life.

62

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

Nomenclature an a0 c0 C1 C2 CD CTSE da/dN DTA FD

Fourier series coefficients initial crack depth initial crack length Paris constant for region I Paris constant for region II central differencing complex Taylor series expansion crack propagation rate damage tolerance analysis Fourier differentiation

Fawaz and Harter studied the impact of various parameters on DTA estimates using deterministic parameter studies [5]. The expected variation in the parameters was estimated using engineering judgment. The cycles-to-failure was used as a metric of importance. Five different cracking scenarios were studied of a transport aircraft fuselage crack under remote tension: a single through crack at a hole, a double through crack at a hole, a single corner crack at a hole, a double corner crack at a hole, and a double oblique through crack at a hole. The material properties and their variation were chosen as representative of 2024-T3 aluminum sheet. Eight parameters were deemed important and, per the analysis results, divided into first and second order effects. The first order parameters were the initial flaw size, the geometry correction factor, the load interaction model, the crack growth rate data, and the stress intensity factor. The secondary parameters were the yield stress, fracture toughness, and threshold stress intensity factor. The distinction between first and second order effects was based on whether the life-cycle costs could be reduced via a more appropriate inspection schedule or if flight safety was affected. The parameters within each category were not comparatively ranked. Millwater and Wieland considered the probabilistic sensitivities of a T-38 wing (corner crack growing from a fastener hole) with respect to initial crack size and aspect ratio, hole diameter, crack growth rate, hole edge distance, geometry correction factor, fracture toughness, retardation and loading spectrum [6]. The conclusion from their study was that the probability-of-failure was sensitive to the fastener hole diameter, b, the crack growth rate at higher DK’s, and the stress spectrum scale factor. Of these variables, only the geometry correction factor was sensitive to variation in both its mean and standard deviation. Probabilistic approaches, while more comprehensive than deterministic approaches, are more arduous since significantly more data are required to determine the probability distributions of the inputs and the analysis is more time consuming to execute. Also, the results typically identify, after the fact, unimportant variables for which the data collection effort was not warranted. Therefore, fast yet accurate deterministic sensitivity methods have an important role to play. Numerical finite differencing is a straightforward commonlyused method to evaluate the derivatives of implicit functions. The method is simple in concept: change a parameter, rerun the analysis and determine the ratio of change in cycles-to-failure to the change in the input parameter. However, an estimate of the derivative is only accurate when the step size is small. On the other hand, when subtracting near-equal numbers, machine round off can also introduce error. As a result, the method is sensitive to the step size, which cannot be too large or too small. This becomes even more of a problem when calculating higher order derivatives since more subtraction operations are required. In summary, finite differencing, while easy to do, is prone to numerical issues that

FFT h i m1 m2 N Nf Dc DK b

fast Fourier transform step size or sampling radius of numerical differentiation square root of negative one Paris exponent for region I Paris exponent for region II number of terms in truncated Fourier series number of cycles to failure stride of the differential equation solver stress intensity factor stress intensity geometry correction factor

may be difficult to discern. It is typical that a laborious trial-anderror effort is required to locate a step size that results in a satisfactory derivative. Alternatives to finite differencing methods are complex variable sensitivity methods, in particular, complex Taylor series expansion (CTSE) and Fourier differentiation (FD). CTSE was originally described by Lyness and Moler [7,8] and was brought to the attention of the engineering community by Squire and Trapp [9]. In CTSE, the first derivative is calculated by perturbing the parameter of interest along the imaginary axes. For example, the initial crack size is given an imaginary perturbation, e.g., a0 + ih, where i denotes an imaginary number and h the step size. As derived below, the sensitivity is then estimated by evaluating the imaginary component of the cycles-to-failure and dividing by the step size. Consequently, CTSE involves no difference operations, thus, allowing for the step size to be made arbitrarily small. Hence the issue of choosing an accurate step size is avoided, making CTSE an easy to implement and highly accurate method for the numerical calculation of first derivatives. In order to calculate the higher order derivatives using CTSE, additional sample points along the imaginary axis are needed, then, in this case, differencing operations are necessary and the step size must be chosen carefully. Fourier differentiation (FD) is analogous complex variable sensitivity method for determining higher order derivatives using sampling in the complex plane. This technique was first described by Lyness et al. [7,8] FD requires the evaluation of sample points along a circular contour around the initial point in the complex plane and a fast Fourier transform (FFT) of the evaluated sample points is used to calculate the derivatives. The use of the FFT as a method to calculate derivatives was described by Lyness et al. [10] and by Henrici [11] and more recently by Bagley [12] to compute of sensitivities of implicit functions. The advantage of FD is that it can compute higher order derivatives more accurately than CTSE or finite differencing. The number of derivatives that can be obtained from the FFT is related to the number of sample points chosen. CTSE has been applied in several engineering fields but as yet is not widely known nor applied. In fluid dynamics, CTSE has been used to find sensitivities for the solution of the Navier–Stokes equation [13]. Furthermore, researchers have been able to apply CTSE techniques to finite element methods in the field of aerodynamics and aero-structural analysis [14,15]. CTSE has also been applied in the study of heat transfer [16], dynamic system optimization [17], pseudospectral [18] and eigenvalue sensitivity methods [19]. CTSE has been compared with automatic differentiation and shown that it is equivalent to the forward mode with comparable accuracy and much simpler implementation [20]. CTSE has been implemented within a solid mechanics finite element program to compute shape sensitivities with good accuracy and it was shown that the accuracy of the derivatives were a function

63

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

of the finite element mesh [21]. Each of these applications uses the same mathematical CTSE approach, thus, substantiating the generality of the method. A review of the literature has not found instances of CTSE being applied to fatigue problems. Standard finite difference techniques have been used in fatigue problems [22,23], as well as probabilistic sensitivity methods [6,24]. However, by adopting complex variable numerical differentiation techniques, it is possible to obtain more accurate sensitivity estimates. It is worth noting that several alternatives to numerical differentiation exist including the equation-based methods of direct differentiation and adjoint differentiation, as well as the code-based method of automatic differentiation. These methods calculate derivates via explicit differentiation and can offer high accuracy and are highly efficient for systems with a large number of input or output variables. However, these methods require either the lengthy derivation of auxiliary equations or heavy modification of existing code and as such have a high cost-of-entry. In contrast, complex variable sensitivity methods are straightforward in concept and offer excellent accuracy; with the caveat that the fatigue crack growth evaluation must be computed using complex variables. For this reason, this paper compares complex variable sensitivity and finite differencing methods. The goal of this paper is to explore the use of complex variable methods as an alternative to the finite differencing method in the determination of the sensitivities of cycles-to-failure with respect to input quantities such as loading amplitude, crack growth parameters, initial crack size, and fracture toughness, for typical fatigue analyses.

approximation of derivatives as a difference between two nearly equal numbers leads to error due to the truncation of terms in the function’s Taylor series. This error can be eliminated by making the step size as small as possible. However, as the step size gets very small, a new source of error arises. This new error is roundoff error and it is due to the fact that a computer cannot accurately calculate a small difference between two near-equal numbers. This means that for finite differencing there is a lower limit on the step size and also a limit on the maximum achievable accuracy. Hence, trial and error is often required to determine a step size that is neither too large nor too small. Higher order derivatives can be calculated through CD by using additional sample points. For example, using three analyses, the formula for the second derivative is.

2. Methodology

f ð2Þ ðx0 Þ 

2.1. Numerical differentiation Finite differencing is a well-known method to obtain an estimate of a function’s derivative that is simple in concept and application. A derivative is defined as the limit (provided the limit exists) of the change in a function’s value across two different points as the distance between the two points goes to zero.

f 0 ðx0 Þ ¼ lim x!x0

f ðxÞ  f ðx0 Þ x  x0

ð1Þ

Finite differencing methods estimate derivatives by approximating the limit in Eq. (1) as a difference between a function evaluated at two distinct points located a distance h apart divided by h .

f 0 ðx0 Þ 

f ðx0 þ hÞ  f ðx0 Þ h

ð2Þ

This distance, h is often called the step size. When h is positive, the method is referred to as forward differencing. When h is negative it is called backwards differencing. When the forward difference and the backwards difference results are averaged, the method is called central differencing (CD). The equation for central differencing is

f 0 ðx0 Þ 

f ðx0 þ hÞ  f ðx0  hÞ 2h

Fig. 1. Schematic of central differencing.

f ðx0 þ hÞ  2f ðx0 Þ þ f ðx0  hÞ 2

h

where the superscript (2) denotes the second derivative. One of the problems with CD is that the calculation of higher order derivatives requires more sample points and more difference operations. Each additional difference operation results in an increase in the round-off error, which further restricts the lower limit of h. This means that CD is often not a good choice for the calculation of higher order derivatives. CTSE is a numerical differentiation method similar in spirit and concept to finite differencing but with significant advantages. CTSE uses the orthogonality of the real and imaginary axes of the complex plane to calculate derivatives with fewer difference operations and in turn less round-off error when compared to CD. Similar to finite differencing, CTSE requires an additional analysis but with a small perturbation along the imaginary axis. That is, variable X = x0 is perturbed to X = x0 + ih, where i denotes an imaginary number and h denotes the step size, see Fig. 2. The formulae for the derivatives can be derived from the Taylor series representation of the function evaluated at the complex sample point.

ð3Þ

The required analyses for CD in the complex plane are shown in Fig. 1. The diamonds indicate that two analyses are required around the sample point X, that is, the real value of X is perturbed by plus and minus h and the function is evaluated at the perturbed values. Finite differencing is very straightforward as no source code changes are required; the only effort is to perform multiple analyses. Although simple and often effective, the finite difference method suffers from one significant drawback. A small step size is needed to accurately represent the derivative. However, the

ð4Þ

Fig. 2. Schematic of the complex Taylor series expansion.

64

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

f ðx0 þ ihÞ ¼ f ðx0 Þ þ f ð1Þ ðx0 Þ

ih ðihÞ2 ðihÞ3 þ f ð2Þ ðx0 Þ þ f ð3Þ ðx0 Þ 1! 2! 3!

þ H:O:T:

ð5Þ

(1)

(2)

where f denotes the first derivative, f the second, etc. Taking the imaginary part of both sides of Eq. (5), solving for the first derivative and ignoring terms of O(h2) and higher yields an estimate of the first derivative.

f ð1Þ ðx0 Þ 

Imðf ðx0 þ ihÞÞ 2 þ Oðh Þ h

ð6Þ

Since no difference operation is needed for the first derivative, the step size can be made arbitrarily small with no concern about increasing round-off error. As a result, highly accurate estimates of the first derivative are obtainable since the estimate is O(h2) accurate and h can be made as small as desired. The estimate of the function evaluated at x0 can be obtained from the real portion of the Taylor series. 2

f ðx0 Þ  Reðf ðx0 þ ihÞÞ þ Oðh Þ

ð7Þ

Since h can be made arbitrarily small, the function’s estimate is recovered numerically. Higher order derivatives can be computed using complex Taylor series expansion. Taking the real part of Eq. (5), the formula for the second derivative with error O(h2) can be derived as

f ð2Þ ðx0 Þ 

2ðf ðx0 Þ  Reðf ðx0 þ ihÞÞÞ h

ð8Þ

2

Note, however, that the second derivative contains a difference operation meaning that round-off error will be a problem if h is set too small. By using more sample points along the imaginary axis, it is possible to solve Eq. (5) to obtain the higher order derivatives For example, if two sample points are used along the imaginary axis in addition to the original sample, the real and imaginary terms of the Taylor series can be used to create a system of equations that can be solved for the desired derivatives. For example, considering two points along the imaginary axis, 3

5

h h h  f ð3Þ ðx0 Þ þ f ð5Þ ðx0 Þ þ . . . 1! 3! 5! 3 2h ð2hÞ ð2hÞ5  f ð3Þ ðx0 Þ þ f ð5Þ ðx0 Þ þ ... Imðf ðx0 þ 2ihÞÞ ¼ f ð1Þ ðx0 Þ 1! 3! 5! 2 4 h h Reðf ðx0 þ ihÞÞ ¼ f ðx0 Þ  f ð2Þ ðx0 Þ þ f ð4Þ ðx0 Þ þ . . . 2! 4! 2 ð2hÞ ð2hÞ4 Reðf ðx0 þ 2ihÞÞ ¼ f ðx0 Þ  f ð2Þ ðx0 Þ þ f ð4Þ ðx0 Þ þ ... 2! 4!

Imðf ðx0 þ ihÞÞ ¼ f ð1Þ ðx0 Þ

Solving the system of equations yields an estimate for the first four derivatives.

f ð2Þ ðx0 Þ  f ð3Þ ðx0 Þ  f ð4Þ ðx0 Þ 

1 2p

Z

1

f ðx0 þ ceih Þeinh dh ¼

1

f ðnÞ ðx0 Þcn n!

4=3Imðf ðx0 þ ihÞÞ  1=6Imðf ðx0 þ 2ihÞÞ h 9=4f ðx0 Þ  7=3Reðf ðx0 þ ihÞÞ  1=12Reðf ðx0 þ 2ihÞÞ

Although Fourier differentiation appears more intricate than CTSE, operationally it is similar in that multiple analyses are required in the complex plane evaluated at N sample points along a circular contour c in the complex plane centered on the initial point, see Fig. 3. The vector of the output data is run through an FFT routine and the output is the first n derivatives in the function’s Taylor series. The nth order derivative of the function can then be calculated from the Taylor series coefficients by using the following relationship.

f ðnÞ ðz0 Þ ¼

an n! cn

ð12Þ

where an is the nth Taylor series coefficient. Due to the orthogonality of the Fourier series terms, the order of the estimation error is related to the number of sample points chosen. For example, when the first eight terms of the Taylor series are computed, then the estimation error will be O(h8) for each derivative. This allows for the potential use of much bigger sampling radii for FD than the step size for finite differencing or CTSE, i.e. c > h, while maintaining a similar level of accuracy. In general, the sampling radius can be at least an order of magnitude larger than the step size for CD or CTSE when using eight sample points. The FFT also relies on difference operations, which unfortunately means that error develops when the sampling radius becomes too small. Note that, while additional analysis points for CD or CTSE are farther from the origin, in contrast, all points using FD are equal-distant from the origin.

A complex variable fatigue analysis software (CVGROW) was developed in Matlab that contains the functions of a typical crack growth fatigue analysis. In particular, the software contains weight functions for stress intensity factor calculations for surface, corner, and through cracks in a finite thickness plate as well as round bar solutions; Paris equation, bilinear Paris, and sigmoidal crack

2

h 2Imðf ðx0 þ ihÞÞ  Imðf ðx0 þ 2ihÞÞ 3

h 6f ðx0 Þ  8Reðf ðx0 þ ihÞÞ þ 2Reðf ðx0 þ 2ihÞÞ h

ð11Þ

2.2. Complex variable fatigue analysis

ð9Þ

f ð1Þ ðx0 Þ 

periodic complex function. This is accomplished by adding a periodic, oscillatory, complex component to the function’s real independent variables. The resulting periodic function now has a Taylor series representation that takes on the properties of a periodic Fourier series. Fast Fourier transform techniques are used to determine the coefficients of this series. The resulting coefficients contain the function’s derivatives. A central feature of Fourier differentiation is to make the individual terms in the Taylor series oscillate at different frequencies. To achieve these oscillations, the perturbation to the variable x is taken to be a complex number as Dx = ceih, where c is a sampling radius. This expression extracts from the series the nth coefficient that contains the nth derivative of the function,

4

ð10Þ However, the quality of the derivatives using this approach will depend on the value of h. FD is a superior method for obtaining higher order derivatives [12]. The heart of FD is making a real valued function become a

Fig. 3. Schematic of analyses required for Fourier differentiation.

65

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

growth laws; variable amplitude loading including rainflow analysis; and unstable fracture or net section yield failure modes. CVGROW has been benchmarked against AFGROW and NASGRO for a number of scenarios and obtains cycles-to-failure results consistently within ten percent (often less than one percent) of the AFGROW and NASGRO results (usually between the two). The crack size versus cycles is determined by solving the crack growth equation using a forward Eulerian differential equation solver with the crack size increment, herein denoted as the ‘‘stride’’, defined as the independent variable. The differential equation solver returns the number of cycles-to-failure and the crack growth histories. For each crack growth increment, CVGROW increments the fastest propagating crack degree of freedom by a distinct crack growth increment defined as the current crack size times Dc, computes the cycles (DN) required to grow this degree of freedom by the designated crack growth increment, then grows the other degrees of freedom by the computed cycles DN, and finally updates the stress intensity factors. A check is then done to see if a surface crack has transitioned to a corner or through crack, or if a corner crack has Table 1 Complex variable crack growth code that will compute the sensitivities of the crack size with respect to various parameters. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

PROGRAM cvgrow IMPLICIT NONE ! Define Constants REAL:: Pi = 4.⁄ATAN(1.0) REAL:: zero = 0.0 REAL:: delN ! crack growth increment REAL:: ncycle = 0.0 ! counter for # cycles COMPLEX:: COMPLEX:: COMPLEX:: COMPLEX:: COMPLEX::

a ! current crack size c ! Paris intercept m ! Paris exponent delsig ! loading beta ! geometry correction factor

REAL:: h = 1.0E-20 ! imaginary step size delN = 0.01 !increment in cycles ! — Define imaginary component for parameter of interest — a = CMPLX(0.01, h) ! computing derivatives wrt initial crack size c = CMPLX(7.0E-11, zero) m = CMPLX(4.0, zero) delsig = CMPLX(80.0, zero) beta = CMPLX(1.12, zero) ! — Open output file — OPEN(UNIT=4,FILE=’datafile.txt’) ! — Grow the crack until 1000 cycles.— DO WHILE (ncycle <= 1000.) ncycle = ncycle + delN CALL grow(a) ! grow the crack ! — Write the cycles, crack size and derivative WRITE(4,⁄) ncycle,REAL(a),IMAG(a)/h END DO CONTAINS ! — Define crack growth routine — SUBROUTINE grow(a) COMPLEX, INTENT(INOUT):: a a = a + (c⁄(delsig⁄beta⁄SQRT(Pi⁄a))⁄⁄m)⁄delN END SUBROUTINE grow END PROGRAM cvgrow PROGRAM cvgrow

transitioned to a through crack; if it has, the stress intensity factor is recalculated using the appropriate model. CVGROW then checks Table 2 Sensitivity comparison of crack size at 1000 cycles: finite difference and CTSE. Sensitivity

Exact

Finite difference (forward differencing)

CTSE

da(1000)/da0 da(1000)/dC da(1000)/dm da(1000)/dr da(1000)/db

3.250 2.067E+08 0.0419 7.235E4 0.05168

3.250 2.084E+08 0.04215 7.410E4 0.05260

3.246 2.064E+08 0.04186 7.222E04 0.05159

Table 3 Source code to compute the sensitivity of cycles-to-failure using CTSE. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

PROGRAM CTSEgrow2 IMPLICIT NONE ! Define Constants REAL:: Pi = 4.⁄ATAN(1.0) REAL:: zero = 0.0 REAL:: Kc = 50.0 ! Fracture toughness COMPLEX:: COMPLEX:: COMPLEX:: COMPLEX:: COMPLEX::

a ! current crack size c ! Paris constant m ! Paris exponent delsig ! loading beta ! geometry correction factor

COMPLEX:: acr ! critical crack size COMPLEX:: delN, ncycle, dela REAL:: h = 1.0E-20 ! imaginary step size ! — Define imaginary component for parameter of interest — a = CMPLX(0.010,zero) ! computing derivatives wrt initial crack size c = CMPLX(7.0E-11,zero) m = CMPLX(4.0,zero) delsig = CMPLX(80.0,zero) beta = CMPLX(1.12,zero) ncycle = CMPLX(zero,zero) acr = (Kc/beta/delsig)⁄⁄2/Pi dela = 0.00001⁄a delN = dela/(c⁄(delsig⁄beta⁄SQRT(Pi⁄a))⁄⁄m) ! — Open output file — OPEN(UNIT=4,FILE=’datafile.txt’) ! — Grow the crack until a > acr — DO WHILE (REAL(a) < REAL(acr)) ncycle = ncycle + delN CALL grow(a) ! grow the crack END DO ! — Write the crack size (RE(a)) and the derivative (IM(a)) vs. cycles — WRITE(4,⁄) REAL(ncycle),IMAG(ncycle)/h CONTAINS ! — Define crack growth routine — SUBROUTINE grow(a) COMPLEX, INTENT(INOUT):: a a = a + (c⁄(delsig⁄beta⁄SQRT(Pi⁄a))⁄⁄m)⁄delN END SUBROUTINE grow END PROGRAM CTSEgrow2

66

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

if failure has occurred through either unstable fracture or net section yield. Remarkably, CVGROW is largely a standard fatigue crack growth software and did not require many adjustments to ensure compatibility with the complex variable methods. In fact the only significant change was to ensure that only the real part of a complex value was used for inequalities, such as when comparing the stress intensity factor against the fracture toughness for evaluating unstable fracture. For example, if the initial crack size, crack growth law, or loading is complex, the stress intensity factor becomes complex. As a result, only the real part of the stress intensity factor is compared against the real part of the fracture toughness. In all cases, the cycles-to-failure becomes a complex variable. The real portion contains the cycles-to-failure and the imaginary portion divided by the step size contains the first derivative of the cycles-to-failure with respect to the parameter of interest.

Table 4 Sensitivities of number of cycles to failure calculated by code in Table 3.

dNf/da0 dNf/dC dNf/dm dNf/dr dNf/db

Exact

Finite difference (forward difference)

CTSE

224,580 2.88E13 6333 106.6 7616

213,800 2.81E13 5001 77.2 4944

201,562 2.887E13 5587 101.0 7205

3. Numerical sensitivity results 3.1. Academic example A simple numerical example is presented first to fix ideas. Table 1 shows a complete Fortran90 program that will compute the derivative of the crack size at time t with respect to several inputs: the initial crack size (a0), Paris constant (C), Paris intercept (m), applied load (r), and geometry correction factor (b). The fracture model is a through crack in a semi-infinite plate under constant amplitude loading modeling using a Paris crack growth equation. The parameters for which to compute the sensitivities are declared as complex type in lines 11–15. Lines 22–26 define which sensitivity will be computed; the imaginary component of the parameter of interest is set to the step size h, defined in line 17, the imaginary component of all other parameters is set to zero. In the example source code in Table 1, the derivative with respect to the initial crack size will be computed. To compute other derivatives, say da(t)/dC, the imaginary component of C is set to h and the imaginary component of a (and all other parameters) is set to zero. The step size h is defined with magnitude 1E-20 in line 17. Lines 32–37 grow the crack. Line 36 writes the crack size, REAL(a), and the sensitivity, IMAG(a)/h, for each time step. Of particular note is that h can be very small, 1E-20 here. This value is approximately nine orders of magnitude smaller than the smallest parameter, the Paris constant. Smaller values can be used without degrading the results, up to the smallest single precision representation of approximately 1038. Smaller values can be used for double precision. As described for CTSE, since the

Table 5 Material properties and initial crack conditions. Property

Value

Initial crack length (c0) Initial crack depth (a0) Part width (W) Part thickness (T) Fracture toughness (KIC) Region I Paris constant (C1) Region I Paris exponent (n1) Region II Paris constant (C2) Region II Paris exponent (n2) Region I to region II transition

0.005 in. 0.005 in. 0.5 in. 1 in. 100 ksi in.0.5 1.45E22 ksi in.0.5 25 ksi in.0.5 7E11 ksi in.0.5 4 3.6 ksi in.0.5

Table 6 First order sensitivities of cycles-to-failure for Dc equal to 0.01. Sensitivity with respect to

CD

CTSE

FD

c0 a0 C1 n1 C2 n2 Max stress (ground air ground cycle)

7.0922E8 1.3059E9 7.7526E15 9.3632E5 5.7061E27 9.9283E5 6.3530E5

7.0795E8 1.3059E9 7.7526E15 9.3630E5 5.7061E27 9.9253E5 6.3406E5

7.0845E8 1.3059E9 7.7526E15 9.3631E5 5.7061E27 9.9268E5 6.3489E5

Fig. 4. Variable amplitude load history.

67

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73 Table 7 Second order sensitivities of cycles-to-failure for Dc equal to 0.01. Sensitivity with respect to

CD

CTSE

FD

c0 a0 C1 n1 C2 n2 Max stress (ground air ground cycle)

7.2638E11 1.9756E12 2.1929E26 1.6771E6 7.8188E49 1.1747E6 4.8061E5

2.1998E11 1.9756E12 2.1929E26 1.6771E6 7.8188E49 1.1745E6 4.1936E5

4.4059E11 1.9756E12 2.1929E26 1.6771E6 7.8188E49 1.1746E6 4.6214E5

Table 8 Comparison of first order sensitivities of cycles-to-failure with respect to crack dimension (a0) for three different strides.

Dc

CD

CTSE

FD

0.02 0.01 0.001

1.3717E9 1.3059E9 1.2490E9

1.3712E9 1.3059E9 1.2489E9

1.3712E9 1.3059E9 1.2489E9

method does not depend upon the subtraction of two near-equal numbers, as for finite differencing, very small values for h can be used and the quality of the results is unaffected by the value of h.

Table 9 Comparison of second order sensitivities of cycles-to-failure with respect to crack dimension (a0) for three different strides.

Dc

CD

CTSE

FD

0.02 0.01 0.001

2.110E12 1.9756E12 1.8188E12

2.1110E12 1.9756E12 1.8542E12

2.1110E12 1.9756E12 1.8471E12

The sensitivity results computed using the program shown in Table 1 are given in Table 2 for a time of 1000 cycles and compared with exact and finite difference estimates. The exact solution was obtained by solving the differential crack growth equation symbolically, differentiating the exact solution with respect to the parameter of interest, then substituting numerical values. As seen in Table 1, the finite difference and CTSE methods were based upon a numerical solution to the differential equation. The calculation of the derivatives using the finite difference method required an iterative trial-and-error approach to find a step size that gave good accuracy. Unlike finite differencing, CTSE is largely step size independent and did not require any iteration on step size. The CTSE and finite difference results agree with the exact solution with the CTSE results typically more accurate. Although the numerical results are only given for N = 1000 cycles, it should be clear that the full time history of the sensitivities is available.

Fig. 5. The First derivative of the number of cycles to failure with respect to a0 calculated by CD.

68

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

A slightly more complicated example is shown in Table 3 that computes the sensitivity of the cycles-to-failure with respect to a parameter of interest. The key for this example is that in line 17, the cycle increment (delN), and total number of cycles (ncycle) are declared complex. Line 30 defines the critical crack size. Line 34 defines the crack growth cycle increment in terms of the current crack size, which will be complex if any of a0, C, m, r or b are complex. The lines 40–43 grow the crack until failure, that is, when the current crack size, contained in the real component of variable a, exceeds the real component of the critical crack size. After failure, the cycles-to-failure and the derivative of the cycles-to-failure with respect to the parameter of interest are written to the output file. Table 4 shows the sensitivities of the cycles-to-failure with respect to various inputs for the exact, finite difference, and CTSE estimates. The accuracy of the finite difference results were susceptible to the step size used so trial and error was employed to get reasonable results. The CTSE results were step size independent.

3.2. Variable amplitude surface crack example A more realistic numerical example was analyzed with the properties shown in Table 5. The geometry was a semi-elliptical

surface crack (surface length 2c and depth a) growing in a finite width plate. The initial surface crack had a one-to-one aspect ratio. A bilinear Paris Law relationship was used to model the da/dN vs. DK curve. Material properties for the model were taken from experimental lab data for Ti–6Al–4V [25]. The variable amplitude tensile loading spectrum, shown in Fig. 4, consisted of 70 load pairs representative of loading on an aircraft engine. A weight function methodology was used to compute the stress intensity factors at the surface (c) and deepest (a) points of the crack. For this analysis, sensitivities of cycles-to-failure, Nf, were calculated with respect to the following input parameters: c0, a0 (initial crack size), C1, C2, n1, n2 (crack growth law), and the maximum load pair of the applied stress. CD required two real analyses with a slight perturbation in the parameter of interest; CTSE required one complex analysis with a small perturbation in the imaginary axes of the parameter of interest. For FD, 16 runs were conducted with perturbations in both the real and imaginary components. The first order sensitivities of the number of cycles-to-failure with respect to these parameters are shown in Table 6, while the second order sensitivities are shown in Table 7. It is readily seen in these tables that the sensitivities agree well between the three methods – there is a moderate discrepancy in the second order derivatives with respect to c0 due to numerical instabilities in the

Fig. 6. The second derivative of the number of cycles to failure with respect to a0 calculated by CD.

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

derivatives arising from the finite stride of the equation solver. As will be shown, this instability is more problematic for CD and FD than for CTSE. All three methods provide fairly stable results for low order sensitivities with respect to the Paris constants, Paris exponents, and the maximum applied stress when using a proper step size, but, as will be shown, may provide unstable results for sensitivities with respect to the initial crack dimensions. In order to examine the issue of stability in more depth, sensitivities with respect to a0 are documented below for three different strides. Table 8 shows the first order derivatives of the cycles-to-failure with respect to the initial crack size depth, for CD, CTSE and FD and Table 9 shows the second order derivatives. The numerical results indicate that all methods produce similar results. The sensitivities are given as a function of the stride used to solve the crack growth differential equation. The smoothness of the cycles-to-failure from the differential equation solver can be adjusted by setting the stride of the solver, defined symbolically as Dc. Decreasing Dc increases the amount of time required to complete an analysis since smaller steps are used but results in smoother and more accurate solutions. Although all three methods produce comparable results, it is informative to examine the stability of the sensitivities over a range of initial values. It is assumed in numerical differentiation

69

that the function being differentiated and its derivatives are smooth functions. This is not necessarily the case for fatigue analysis since the function being differentiated is a numerical solution to a first order differential equation (through crack) or a set of first order coupled differential equations (surface crack). This lack of smoothness may prevent the accurate calculation of derivatives of an arbitrary order but it may be possible to calculate derivatives up to a certain order if the numerical solution is sufficiently smooth over a local region around the initial point of differentiation. In order to test the validity and performance of the three numerical differentiation techniques, first through fourth order derivatives were calculated at several points within a small neighborhood. Fig. 5 displays the first order sensitivity of the number of cyclesto-failure with respect to a0 as calculated by CD, with h = 0.001  a0 for three different values of Dc: 2%, 1%, and 0.1% and Fig. 6 displays the second order sensitivity calculated by CD for the same values of Dc. In these figures, all input parameters were held constant through the analysis except for a0, that is, multiple first order derivatives were computed over a local range of a0. Fig. 6 shows that there is some instability in the calculation of the second order derivative using CD unless Dc is sufficiently small. That is, small changes in the initial crack size can cause a step change in the cycles-to-failure resulting in a large error in the estimate of the

Fig. 7. The first derivative of the number of cycles to failure with respect to a0 calculated by CTSE.

70

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

second order derivative for particular initial crack size values. The errors reduce in magnitude if Dc is sufficiently small but may increase in frequency. This is seen in Fig. 6C. However, in practice, the correct value for Dc is not known a priori and a smaller value of the crack growth stride results in longer solution times. The first and second order sensitivities calculated by CTSE with h = 0.001  a0 for the three different values of Dc are shown in Figs. 7 and 8, respectively. The CTSE results agree closely with the CD results for the first order derivatives (less than 2% relative difference) as seen in Table 7 and by comparing Figs. 5 and 7. Fig. 8 shows that the second order sensitivities are very smooth over the entire domain, unlike CD; compare Figs. 6 and 8. This stability is seen at every value of Dc tested. Thus, CTSE is at least as good as CD for the first order derivative and far superior for the second order derivative. Derivatives were also calculated using FD with the sampling radius of h = 0.001a0 and 16 sampling points and are given in Tables 7 and 8. By utilizing points along both the real and imaginary axes, FD results in sensitivities with instability that falls somewhere between the instability seen in CD and CTSE. FD also comes with an increased computational cost as compared to CD and CTSE, which means that it may not be the best choice for calculating low order sensitivities.

First and second order sensitivities were also calculated with respect to several other parameters such as material properties and the maximum applied stress. In all of these cases, CTSE returned smoother first and second order derivatives than CD or FD. Based on these results, CTSE is the best choice for first and second order sensitivity calculations. The biggest advantage of FD over CTSE and CD is that it allows for the accurate calculation of higher order derivatives, e.g., third order and higher [12]. In the case of CD and CTSE, the number of difference operations required for high order numerical differentiation makes the method error prone and unstable. In addition, each additional sample point for CD and CTSE is located further away from the initial point, whereas for FD all additional sample points are located at the same distance; see Fig. 3. CD, CTSE and FD were used to calculate sensitivities of the cycles-to-failure with respect to the initial crack size through four orders. Fig. 9 shows the stability of the fourth derivative with respect to a0 as calculated by CD with h = 1% of the initial value (Fig. 9A), CTSE with h = 1% of the initial value (Fig. 9B), and FD with a sampling radius of h = 5% and 16 sample points (Fig. 9C). Note that for higher order derivatives the step size h has to be much larger than lower order derivatives in order to prevent the accumulation of error due to machine roundoff. In each case, the sensitivities

Fig. 8. The second derivative of the number of cycles to failure with respect to a0 calculated by CTSE.

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

71

Fig. 9. The fourth derivatives of the number of cycles to failure with respect to a0. A, B – central differencing; C, D – complex Taylor series expansion; E, F – Fourier differentiation.

were calculated for two different Dc’s, 1% and 0.1%. In general, as the order of differentiation increases, the stability and accuracy of the numerical differentiation is reduced.

As confirmed in Fig. 9, FD generates more stable fourth order derivatives than CD and CTSE since FD generates stable results at both sizes of Dc whereas CD and CTSE exhibit large point instabil-

72

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

Fig. 10. The second order mixed term sensitivity of the number of cycles to failure with respect to C2 and m2 (@N 2f =@C 2 @m2 ). A – central differencing; B – complex Taylor series expansion; E, F – Fourier differentiation.

ities, although CTSE is superior to CD. In fact using the smaller Dc = 0.001, FD returns a remarkably smooth curve over the range of initial crack depths. CD, CTSE, and FD can all be used to calculate mixed term derivatives and, hence, evaluate the importance of interaction terms. Fig. 10 shows the stability of the mixed term sensitivity of the number of cycles-to-failure with respect to the Paris constant and Paris exponent for region II of the bilinear crack growth curve, i.e., @N 2f =@C 2 @m2 . For CD and CTSE, the sampling step size was 0.2% of the initial value of each parameter and for FD the sampling radius was 0.2% and four sample points were used for each variable (16 total sample points). These parameters were chosen as they are typically a good starting point for calculating second through fourth order derivatives. For this case, each method produced a stable mixed term derivative.

4. Conclusions Complex variable numerical differentiation techniques provide a straightforward and accurate method to calculate sensitivities useful for fatigue analysis. CTSE is analogous to finite differencing however the parameter of interest is perturbed along the imaginary axis rather than the real axis. The advantage this provides is

that the first order derivative can be obtained without any differencing operations; thus, allowing the step size to be as small as desired and avoiding the troublesome step size issue that is problematic for finite differencing. FD requires real and imaginary perturbations to the parameter of interest such that samples are obtained in a circle within the complex plane. The complex cycles-to-failure results are then post processed through an FFT routine to obtain the derivatives. A complex variable fatigue analysis software was written that is emblematic of the solutions required for a typical fatigue analysis. The complex variable version required only minor alterations of a traditional real valued fatigue analysis code. The software is complex in that any input (initial crack size, material property, loading) can be a complex variable; as a result, the cycles-to-failure becomes complex. For CTSE, the real portion of the cycles-to-failure variable holds the traditional cycles-to-failure results and the imaginary portion holds the derivative times the imaginary step size. All three methods, CD, CTSE and FD were shown to give accurate first and second order derivatives in general; however, CTSE outperforms the traditional CD method for the calculation of first and second order sensitivities in that the results are more stable, particularly for larger crack growth strides (Dc). In terms of computational cost, CTSE requires half the number of sample points required by CD, albeit each complex sample point requires takes

A. Voorhees et al. / International Journal of Fatigue 40 (2012) 61–73

about three times longer to evaluate. This effectively makes the computational cost of CTSE approximately 1.5 times that of CD. FD clearly outperforms CD and CTSE for the calculation of derivatives above 2nd order. The calculation of high order derivatives through FD requires an increase in the number of sample points required, however this increase also leads to an increase in the stability and accuracy of the low order derivatives. By utilizing more sample points in FD, the stability of the low order derivatives approaches the stability of the low order derivatives as calculated by CTSE. In addition, the complex variable sensitivity methods can also be used to compute partial derivatives of any order. Acknowledgments This work was supported under Air Force Agreement FA865007-C-5060, Dr. Patrick J. Golden, AFRL/RXLMN, Project Monitor. The authors would like to thank Juan Ocampo for many improvements to the crack growth software. References [1] Harter JA. AFGROW users guide and technical manual, AFRL-VA-WP-TR-19993016; 1999. p. 11–2. [2] NASGROW fracture mechanics and fatigue crack growth analysis software, v6.0, NASA-JSC and Southwest Research Institute; 2009. [3] Leverant GR, Millwater HR, McClung RC, Enright MP. A new tool for design and certification of aircraft turbine rotors. J Eng Gas Turbines Power 2004;126:155–9. [4] McGinty RD. Sensitivity of fatigue crack growth rates to operating conditions, sample problem no. MERC-3. In: Miedlar PC, Berens AP, Gunderson AW, Gallagher JP, editors. USAF Damage Tolerant Design Handbook: Guidelines for the Analysis and Design of Damage Tolerant Aircraft Structures, AFRL-VA-WP-TR2003-3002,Wright-Patterson Air Force Base, OH; November 2002 . [5] Fawaz SA, Harter JA. Impact of parameter accuracy on aircraft structural integrity estimates. US Air Force Technical Report, AFRL-VA-WPTR-2001-3057; 2001. [6] Millwater HR, Wieland D. Probabilistic sensitivity-based ranking of damage tolerance analysis elements. AIAA J Aircraft 2010;47(1):161–71.

73

[7] Lyness JN, Moler CB. Numerical differentiation of analytic functions. SIAM J Numer Anal 1967;4(2):202–10. [8] Lyness JN. Differentiation formulas for analytic functions. Math Comput 1968;22(102):352–62. [9] Squire W, Trapp G. Using complex variables to estimate derivatives of real functions. SIAM Rev 1998;40(1):110–2. [10] Lyness JN, Sande G. ENTCAF and ENTCRE: evaluation of normalized taylor coefficients of an analytic function. Commun ACM 1971;14(10):669–75. [11] Henrici P. Fast Fourier methods in computational complex analysis. SIAM Rev 1979;4(11):481–527. [12] Bagley RL. On Fourier differentiation – a numerical tool for implicit functions. Int J Appl Math 2006;19(3):255–81. [13] Anderson WK, Newman JC, Whitfield DL, Nielsen EJ. Sensitivity analysis for Navier–Stokes equations on unstructured meshes using complex variables. AIAA J 2001;39(1):56–63. [14] Newman JC, Whitfield DL, Anderson WK. Step-size independent approach for multidisciplinary sensitivity analysis. J Aircraft 2003;40(3):566–73. [15] Burg COE, Newman JC, Efficient Computationally. Numerically exact design space derivatives via the complex Taylor’s series expansion method. Comput Fluids 2003;32(3):373–83. [16] Gao XW, He MC. A new inverse analysis for multi-region heat conduction BEM using complex-variable-differentiation method. Eng Anal Boundary Elem 2005;29(8):788–95. [17] Kim J, Bates DG, Postlethwaite I. Nonlinear Robust performance analysis using complex-step gradient approximation. Automatica 2006;42(1):177–82. [18] Cervino LI, Bewley TR. On the extension of the complex-step derivative technique to pseudospectral algorithms. J Comput Phys 2003;187(2):544–9. [19] Wang BP, Apte AP. Complex variable method for eigensolution sensitivity analysis. AIAA J 2006;44(12):2958–61. [20] Martins JRRA, Sturdza P, Alonso JJ. The complex-step derivative approximation. ACM Trans Math Softw 2003;29:245–63. [21] Voorhees A, Millwater HR, Bagley RL. Complex variable methods for shape sensitivity of finite element models. Finite Elem Anal Des 2011;47:1146–56. doi:10.1016/j.finel.2011.05.003. [22] Zeiler TA, Barkey ME. Design sensitivities of fatigue performance and structural dynamic response in an automotive application. Struct Multidiscip Optimizat 2001;21(4):309–15. [23] Zeiler TA. Use of structural dynamic and fatigue sensitivity derivatives in an automotive design optimization. Struct Multidiscip Optimiz 2002;23(5): 390–7. [24] Castillo E, Fernandez-Canteli A, Hadi AS, Lopez-Aenlle M. A fatigue model with local sensitivity analysis. Fatigue Fract Eng Mater Struct 2007;30(2):149–68. [25] Gallagher JP, van Stone RH, deLaneuville RE, Gravett P, Bellows RS, Slavik DC, et al. AFRL-ML-TR-2001-4159, Improved high-cycle fatigue (HCF) life prediction, Wright-Patterson Air Force Base, OH; 2001.