Engineering Fracture Mechanics 73 (2006) 2662–2684

www.elsevier.com/locate/engfracmech

A weight function methodology for the assessment

of embedded and surface irregular plane cracks
Hugo López Montenegro *, Adrián Cisilino, José Luis Otegui
Institute of Materials Science and Technology (INTEMA), University of Mar del Plata, Welding and Fracture Division,
J. B. Justo 4302, 7600 Mar del Plata, Buenos Aires, Argentina

Received 20 September 2005; received in revised form 2 April 2006; accepted 12 April 2006
Available online 14 June 2006

Abstract

A low cost numerical tool for the calculation of mode I stress intensity factors K in embedded and surface irregular
cracks is presented in this paper. The proposed tool is an extension of the O-integral algorithm due to Oore and Burns
for the assessment of embedded plane cracks using the weight function methodology. The performance of the O-integral
is assessed first by comparing its K results to exact solutions for embedded elliptical and rectangular cracks. From the ana-
lysis of this data it is found that the error in the K results systematically depends on the crack aspect ratio and the local
crack front curvature. Based on this evidence a corrective function is derived in order to remediate the limitations of the
O-integral. Solutions due to Newman and Raju are used to account for the effects of free surfaces and finite thickness. The
accuracy of the proposed procedure is assessed by solving a number of examples and by comparing the obtained results to
those available in the literature.
Ó 2006 Elsevier Ltd. All rights reserved.

Keywords: Fracture mechanics; Stress intensity factor; Three-dimensional cracks; Numerical algorithms; Surface cracks

1. Introduction

Design standards based on theory and rules of good art allow to minimize the risk of failure of structures of
all types. Among those, best known are ASME PVPC Code [1], used by manufacturers and operators of pres-
sure vessels and pipes, and others applicable to specific industries. These standards mostly come from USA
(API, SEAL, STUMP, etc.), and from European countries (DIN, ISO, etc.). Another interesting aspect refers
to the evaluation of equipment in operation that, due to fabrication errors or accumulated damage in service,
show discontinuities or crack-like defects that are not allowed by the manufacturing codes. Nowadays there is
an advanced theoretical understanding and an array of experimental and numerical tools to address this engi-
neering problem; such as failure analyses, root cause assessments, testing of laboratory samples, numerical
analysis for design and re-rating, etc. A set of documents and recommended practices standardize the analysis

*
Corresponding author. Tel.: +54 223 4816600; fax: +54 223 4810046.
E-mail address: hmonte@fi.mdp.edu.ar (H.L. Montenegro).

0013-7944/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved.
doi:10.1016/j.engfracmech.2006.04.007
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2663

Nomenclature

a, c crack dimensions: half depth and half width, respectively

r, r* crack aspect ratio
b half width of a plate
t plate thickness
K stress intensity factor
Kcorr corrected stress intensity factor
KI mode I stress intensity factor
K Q0 stress intensity factor at a point Q 0 on the crack front
r opening stress
r0 remote uniform stress
rQ opening force intensity (pressure) at a point Q
A crack surface area
s crack front
ds elemental length of the crack front s
lQQ0 distance from the load point Q to the point Q 0 on the crack front
qQ distance from point Q to the center of the elemental length ds
W QQ0 weight function
vn non-dimensional normalized curvature of the crack front
/ parametric angle for the crack geometry
r(h), h polar coordinates
fc corrective function
P, R, G defined functions of r
f1, f2 defined functions of vn
p, q, t coefficients in f2, depending on the parameter r
Y configuration function
Yb configuration function for the crack in an infinite body
YC b Yb solution obtained by using the Calyf program
ks free surface correction factor
f/ Newman and Raju solution function for an ambedded elliptical crack
Q parameter approximating the square of the elliptical integral of second type
FS, H Newman and Raju auxiliary functions
rm, rb remote membrane and pure bending stresses, respectively
tm, tb normalized values of rm and rb, respectively
fS auxiliary function for FS
gS, fw auxiliary functions for fS
MS auxiliary parameters for fS
q parameter related to the deformation of a circumference of radius one

of structures in service, which purpose is to prevent accidents and at the same time to reduce repair costs. The
British documents BSI PD 6493 [2], CEGB R6 [3]; and the American documents EPRI (Electric Power
Research Institute) of 1981 and 1990 [4,5] were among the first bodies to address the problem, leading to pres-
ent day procedures such as the British Standard BS 7910 [6], which replaces the PD 6493, the European stan-
dard procedure for structural integrity assessment SINTAP [7] and the American Petroleum Institute
Recommended Practice API RP 579 [8].
Assessing the engineering integrity and life expectancy of a cracked component or structure, either under
service conditions or during the design stage requires the determination of fracture parameters as the stress
intensity factor K. In this sense, most of the above-mentioned codes and standard procedures use assessment
rules with a certain degree of conservatism. For example, API RP 579 models embedded or surface cracks as
2664 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

ellipses or semi-ellipses, respectively, based on their width and depth, as they could be defined by ultrasonic
non-destructive testing. Similar approaches are applied to through cracks in the center or edge of a plate,
where the irregular geometry of the crack front is idealized as a straight fronted crack. In the case of a corner
crack, it is idealized as a quarter ellipse. With respect to multiple cracks, all documents establish approaches
for crack interaction.
Another well-established
pffiffiffiffiffiffiffiffiffiffiffi criterion to estimate stress intensity factors for simple cracks of arbitrary geom-
etry, is in terms of Area. Following this approach, Murakami et al. [9,10] proposed simple formulas for the
estimation of K maximum values in surface or embedded cracks, with a reported accuracy of 10%. In a recent
work, Noda [11] checks the validity of these formulas in rectangular and elliptical cracks of different aspect
ratios, and extends them to the case of mixed mode loading in homogeneous and heterogeneous materials.
When a K solution is required for a crack of arbitrary shape in a complex structural configuration methods
based on finite elements (FEM) or boundary elements (BEM) are widely applicable [12]. Although very effec-
tive, versatile, and capable of delivering solutions with a high level of accuracy, the drawback of these methods
is their relatively high computing cost and time needed for the model preparation.
Since the introduction of the weight function concept to crack problems by Bueckner [13] and Rice [14],
many numerical methods have emerged for the calculation of stress intensity factors. Weight function meth-
odologies are simple and accurate, but they present the drawback that a particular weight function solution is
necessary for each model geometry. This feature limits the applicability of the methodology, especially in the
case of three-dimensional problems, for which available weight function solutions are limited to simple regular
geometries [15]. In an effort to extend the applicability of the weight function methodology, Oore and Burns
[16] proposed a general, geometrically defined weight function for calculating the opening mode stress inten-
sity factor at any point on the front of an irregular planar crack embedded in an infinite solid subjected to an
arbitrary stress field. This methodology known as the O-integral has proved general and versatile. However,
depending on the local crack-front curvature, it could deliver K results with errors of up to 20% for simple
geometries such as elliptical cracks with low aspect ratios.
A corrective procedure to extend the applicability of the O-integral methodology is introduced in this
paper. In a first step, the performance of the O-integral is assessed and a correction function proposed in order
to remediate accuracy problems referenced in the previous paragraph. In a second step, a corrective procedure
is presented to allow the O-integral for the analysis of surface cracks and to account for the effect of finite
thickness. The accuracy of the procedure is assessed by studying a number of examples and comparing their
results to solutions obtained using BEM models and available in the literature.

2. The O-integral

According to the weight function (WF) concept [13], for an opening force intensity (pressure) rQ acting on
the crack surface area A, the stress intensity factor K Q0 at a point Q 0 on the crack front is an integral calculated
on the crack surface area as
Z
K Q0 ¼ rQ  W QQ0  dA ð1Þ
A

where W QQ0 is the weight function, which depends on the problem geometry only. It is worth mentioning that
most of the available weight function solutions are limited to two-dimensional cracks and regular crack geom-
etries in three-dimensions [15], what certainly limits the scope of problems which can be solved using this
technique.
In an effort to extend the applicability of the weight function methodology, Oore and Burns [16] examined
the structure of the weight functions for three-dimensional cracks for which there are known solutions, and
surmised that in general these weight functions could be rewritten as, or approximated by
pffiffiffi
2 1
W QQ0 ¼   ð2Þ
p  l2QQ0 R ds 1=2
S q2Q
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2665

where lQQ0 is the distance from the load point Q to the point Q 0 on the crack front, and qQ is the distance from
Q to the center of the elemental length ds of the crack front s (see Fig. 1(a)). Thus, the O-integral is a geomet-
rically defined weight function for calculating opening mode stress intensity factors at any point on the front
of an irregular planar crack, embedded in an infinite solid and subjected to an arbitrary stress field [17].
It has been shown that this integral satisfactorily predicts KI for a variety of normally loaded, embedded
cracks. Among others, particular solutions have been derived for a circular crack with a pair of opening, sym-
metrical, point loads, an elliptical crack under uniform or linear varying stress field, an annular circular crack
in a uniform stress field, and a parabolic crack in a uniform stress field [16]. For some of the aforementioned
test cases, it was possible to obtain closed-form solutions for Eq. (1).
In the case of cracks of arbitrary shape and loading, the O-integral needs to be evaluated numerically. In
this sense, Desjardins [17] proposed a numerical algorithm that employs interior and border elements for the
discretization of the crack surface area (see Fig. 1(b)). For the interior elements the evaluation of the integral

Q´
A

dA lQQ´
Q
s ρQ

ds

(a)

(b)

Fig. 1. (a) Definition of geometric variables used in the O-integral, (b) discretization strategy. Note the utilization of boundary and
interior elements used by Desjardins (from Ref. [17]).
2666 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

in Eq. (1) is carried out using the standard Gauss integration formula, while for the border elements the solu-
tion is obtained by means of a semi-analytical procedure.
Although general and versatile, it is found that depending on the local crack-front curvature, K results are
systematically under or overestimated by the O-integral. In the following section the performance of the O-
integral is assessed and a corrective function proposed in order to remediate the above-mentioned limitation.
In the same way a corrective procedure is introduced in order to extend the applicability of the methodology
to surface cracks and finite-thickness specimens.

3. Assessment and improvement of the O-integral

3.1. Embedded cracks

Elliptical and rectangular cracks were selected to assess the precision of the O-integral. Exact solutions for
the elliptical crack in an infinite domain are documented in the literature; the Irwin solution for uniform
remote tension [18] and the Shah and Kobayashi solution for pure bending [19].
In the work by Desjardin [17], circular and elliptical cracks under tension and bending stresses were solved
using the O-integral, and the results compared with the exact solutions. Maximum errors were found to be
almost independent from the loading case, but dependent on crack geometry. Maximum error was found
0.65% for the circular crack, while for elliptical cracks the error monotonically increased with the diminution
of the crack aspect ratio (the ratio between the length of the minor and the major axes of the ellipse). Maxi-
mum error was found in the later case equal to 18.38% for an elliptical crack of aspect ratio 0.2. This behavior
allows us to conclude that in the circumferential crack the error can be entirely attributed to the numerical
algorithm, as for this geometry the O-integral solution matches exactly the analytical weight function. In con-
trast, for the elliptical crack the largest error corresponds to the portion of the crack front with largest
curvature. Convergence analyses done using different discretizations showed that this error is due to the O-
integral formulation itself, not a consequence of the numerical evaluation. The fact that the error changes
monotonically with the ellipse aspect ratio and the local crack front curvature suggests the possibility of
obtaining a corrective function for K which depends on these two parameters only.
The aspect ratio of the ellipse with half axes c and a, is defined as r = a/c. At the same time, the non-dimen-
sional normalized curvature vn (from now on, simply referred as ‘curvature’) at a generic front point Q 0 asso-
ciated with the parameter / (see Fig. 2) is given by Eq. (3):
2
ða=cÞ
vn ð/Þ ¼ h i32 ð3Þ
2
sin2 ð/Þ þ ða=cÞ  cos2 ð/Þ

For ellipses with aspect ratio r 6 1, the curvature diminishes monotonically with the parameter /, and it pre-
sents extreme values given by

A
Q´
a
φ
C

2c

Fig. 2. Elliptical crack and its characteristic geometrical parameters.

 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2667

vn ð0Þ ¼ 1=r; vn ðp=2Þ ¼ r2 ð4Þ

where / = 0 corresponds to the major curvature of the ellipse (point C), and / = p/2 to the minor curvature
(point A) in Fig. 2. The corrective function fc is defined as
Kðreference solutionÞ
fc ðr; vn Þ ¼ ð5Þ
KðOoreÞ
Then, the corrected stress intensity factor value Kcorr results from the product of the K(Oore) value and the
corrective function fc
K corr ¼ KðOoreÞ  fc ð6Þ
When the corrective function is calculated as a function of the crack front curvature for elliptical cracks with
aspect ratios 0.2 6 r 6 1, the set of curves shown in Fig. 3 is obtained. Exact solutions due to Irwin [18] were
used as reference in this case. As seen in the figure, all the curves are approximately parallel for curvatures
vn P 1. At the same time, as the aspect ratio approximates to one, the curves converge to the limiting point
L, which corresponds to a circumferential crack with constant curvature vn = 1 and corrective function fc = 1.
The K solutions of rectangular cracks given by Isida et al. [20] were employed as reference solutions to
study the behavior of the corrective function for geometries with null curvature and different aspect ratios.
Isida solutions, obtained by means of the body force method, allow computing KI at point A on the crack
front (see subfigure in Fig. 4) with a maximum error less than 1%. The corrective function for the rectangular
cracks is plotted in Fig. 4. It can be appreciated that fc approximates to one as r tends to zero. This behavior is
coherent with the Oore and Isida solutions. It is worth to note here that the geometrical WF proposed by Oore
exactly matches the closed form solution of the WF for a 3-D tunnel crack in an infinite solid or the centered
crack in a 2-D infinite plate [16].
Based upon the analysis of the above results, an algorithm for the computation of the corrective function fc
is proposed next. The scheme in Fig. 5 illustrates the typical behavior of the curves plotted in Fig. 3, which
present the maximum value of fc in the position with minimum curvature, vn = r2. Two points are defined
next: point R, which is given by the intersection of the curve and the vertical line vn = 1; and point P, which
corresponds to the intersection of the curve with the horizontal line fc = 1. The abscissa of P and the ordinate
of R were computed for all the curves in Fig. 3 (0.2 6 r < 1) using data interpolation. Obtained results are
presented in Fig. 6 as a function of the crack aspect ratio r. Fig. 6(a) corresponds to the results for the ordinate

1.05
aspect ratio r
1.0
L
1.00 0.9
0.8
0.7
corrective function fc

0.6
0.95
0.5
0.4
0.3
0.90 0.2

0.85

0.80
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
curvature χn

Fig. 3. Corrective function for elliptical crack as a function of the curvature.

2668 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

1.15

Oore versus Isida solution

Exact solution for
1.10 3-D tunnel crack (r = 0)

corrective function fc

1.05

y r=a/c
A
1.00 x
2a

2c

0.95
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
crack aspect ratio r

Fig. 4. Corrective function at point A for a family of rectangular cracks embedded in an infinite volume.

fc

P L
1
R

r2 1 χn

Fig. 5. Scheme of a typical curve reported in Fig. 3 for a given aspect ratio 0.2 6 r < 1.

of point R, function R(r), while Fig. 6(b) illustrates the results for the abscissa of point P, function P(r). The
best fitting functions R(r) and P(r) are given in Eqs. (7) and (8), respectively

RðrÞ ¼ 0:8915 þ 0:2953  r  0:2773  r2 þ 0:0905  r3 ð7Þ

P ðrÞ ¼ 0:839  r; 0 6 r < 0:2 ð8aÞ
2 3
P ðrÞ ¼ 0:278554 þ 2:58316  r  1:78777  r þ 0:483164  r ; 0:2 6 r 6 1 ð8bÞ

Note that since fc = 1 for r = 0 for both the rectangular and elliptical cracks, the function P(r) was linearly
approximated in the range 0 6 r < 0.2 (see Fig. 6(b)). On the other hand, the function R(r) was extrapolated
for r < 0.2 (see Fig. 6(a)). It is also worth noting that for r = 1, both P(r) and R(r) tend to one, i.e., the limiting
point L (see Fig. 5).
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2669

1.02
1.01
1.00
0.99
0.98
0.97
0.96
R (r)

0.95
0.94
0.93
0.92
0.91
0.90
0.89
0.88
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(a) aspect ratio r

1.0

0.9

0.8

0.7

0.6
P (r)

0.5

0.4

0.3

0.2

0.1

0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(b) aspect ratio r

Fig. 6. Intersection functions for ellipses: (a) ordinate R(r), (b) abscissa P(r).

As it was mentioned before, curves in Fig. 3 result approximately parallel for curvatures vn > 1. Thus, data
with vn > 1 for all the curves in Fig. 3 was shifted to contain point L by using Eq. (7). Then, a single expo-
nential function f1(vn) was obtained by fitting all the data (see Eq. (9) and continuous line in Fig. 7)
 
ðvn þ 0:68767Þ
f1 ðvn Þ ¼ 0:85076 þ 0:23245  exp ð9Þ
3:80872
Similarly, data with abscissa vn < 1 and only for the curves in the range 0.7 6 r 6 0.9 was also shifted to con-
tain point L. As it can be seen from Fig. 7, this data is adequately represented by the extrapolation of the expo-
nential fit given in Eq. (9) (see dashed line in Fig. 7). On the other hand, data in the curves of Fig. 3 for r < 0.7
does not follow the same pattern, and thus, an alternative fitting scheme for these curves is necessary for
vn < 1. Best results in this case were obtained using a parabolic interpolation function f2(vn) with positive con-
cavity as follows:
f2 ðvn Þ ¼ p  v2n þ q  vn þ t ð10Þ
2670 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

1.05
1.04
ellipses with 0.2 < r < 0.9 ; χ n > 1
1.03
exponential fit ; χ n > 1
1.02 ellipse with r = 0.7 ; χ n < 1
1.01
L ellipse with r = 0.8 ; χ n < 1
1.00 ellipse with r = 0.9 ; χ n < 1
f1 ( χ n ) 0.99 extrapolation ; χ n < 1
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.90
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
crack front curvature χn

Fig. 7. Data from Fig. 3 transferred on the basis of Eq. (7) and the exponential function f1(vn).

with coefficients depending on the parameter r by

1 þ GðrÞ  RðrÞ  GðrÞ  P ðrÞ
p¼ 2
ðP ðrÞ  1Þ
GðrÞ  ðP 2 ðrÞ  1Þ þ 2RðrÞ  2
q¼
ðP ðrÞ  1Þ2
t ¼ RðrÞ  p  q
where
0:85076  RðrÞ
GðrÞ ¼
3:80872
This selection of f2(vn) ensures that the C1 continuity condition is satisfied at vn = 1 (i.e., both interpolation
functions, f1(vn) and f2(vn), have the same slope at vn = 1). In this way, the corrective function fc(r, vn) results
in a piece-wise representation. It is worth to note that the C1 condition and the requirement of positive con-
cavity for f2(vn) allow obtaining the exact value for r* = 0.6895, the limiting point between the two fitting
schemes. In this way, for cracks with aspect ratios r P r* the corrective function fc(r, vn) is computed by means
of Eqs. (7) and (9). On the other hand, if r < r* the corrective function is computed using Eqs. (7) and (9) only
when vn P 1, while for vn < 1 Eq. (10) is employed. Fig. 8 illustrates the later case.
Finally, the algorithm for the computation of the corrective function is as follows:

(a) Compute the crack aspect ratio, r.

(b) Compute the local curvature vn at the point of interest on the crack front. If vn < 0, set vn = 0.
(c) If r P r* = 0.6895, compute the corrective function using Eqs. (7) and (9) by doing
fc ¼ f1  ½1  RðrÞ
(d) If r < r*, compute the corrective function using the procedure given in (c) for vn P 1, or the function f2
given in Eq. (10) for vn < 1.

The O-integral numerical solver and the aforementioned algorithm were implemented in a program named
Calsyf. Preliminary tests with the program have shown that the algorithm provides K results for elliptical and
rectangular cracks under tension with errors less than 1.6% and 1.06%, respectively. This level of accuracy
results very competitive when compared to other simple formulas for the estimation of maximum K values
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2671

1.05

fc = f2 fc = f1 - [1 - R (r )]
aspect ratio r
0.6
1.00 0.5
0.4
corrective function fc 0.3
0.2

0.95

0.90

0.85
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
curvature χ n

Fig. 8. Data from Fig. 3 for r 6 0.6 and correcting function fc defined from the fitting functions f1(vn) and f2(vn).

[9–11], with reported errors of a few percent. Finally, it is worth to mention that the proposed algorithm is not
limited to regular crack geometries. As it will be shown later, the corrective function can be easily and success-
fully extended to solve irregular cracks.

3.2. Surface cracks

3.2.1. Introduction
The O-Integral was developed to solve embedded plane cracks in an infinite domain. However, practical
fracture mechanics problems deal with surface cracks, mostly located in finite thickness plates or shells. In
what follows an algorithm is proposed in order to extend the applicability of the O-integral methodology
for the assessment of surface cracks located in semi-infinite bodies and in finite thickness geometries.
The generic expression of K in a point of the crack front is
pffiffiffiffiffiffi
K ¼ Y r pa ð11Þ
where the non-dimensional factor ‘‘Y’’ is a normalized K, or configuration function, which depends on geo-
metrical parameters and on the loading configuration. When treating surface cracks in a semi-infinite body,
the Y factor can be split in two terms as follows:
Y ! ks  Y b
where factor ks accounts for the correction due to the free surface, and the function Yb corresponds to the
configuration function for the crack in an infinite body. Note that this approach allows using the program
Calsyf to obtain the solution for Yb, named Y Cb , in arbitrary surface cracks under arbitrary loads, provided
that ks for the given geometry is available. From Eqs. (6) and (11), the parameter Y Cb is then expressed by
K corr
Y Cb ¼ pffiffiffiffiffiffi ð12Þ
r pa
As it will be shown, this procedure can be extended to the treatment of finite width plates.
A usual procedure that simplifies the study of surface cracks subjected to arbitrary loads consists in com-
puting the K solution for an equivalent embedded crack in an infinite body. This is done by considering the
free surface as a plane of symmetry. The procedure is illustrated in Fig. 9 for an irregular surface crack. Note
that for an irregular geometry like the one shown in Fig. 9, a generalized aspect ratio can be defined as the
2672 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

2c
A

Q´
a
φ
Free surface C

Crack surface

Fig. 9. Irregular surface crack and its equivalent symmetrical crack.

maximum axis ratio 2a2c

. The O-integral applies the above procedure when dealing with arbitrary cracks and
loads, as shown in the works of Oore and Burns [21] and of Grueter et al. [22], as a means to compute the
Yb factor.

3.2.2. Correction due to free surface and finite size geometry

In the resolution of elliptical or semi-elliptical surface cracks in a finite plate, the solution by Newman and
Raju [23,24] is broadly used. This solution was obtained from finite element results, with reported accuracy of
±5%. Considering a function f/ defined as
 1=4
f/ ¼ sin2 ð/Þ þ ða=cÞ2 cos2 ð/Þ ð13Þ

the Newman and Raju solution for an elliptical crack in an infinite body under unitary uniform stress can be
expressed as (see Appendix A for further details):
p
Y b ¼ f/ = Q ð14Þ
where
1:65
Q ¼ 1 þ 1:464ða=cÞ ð15Þ
pffiffiffiffi
and the root Q approximates the elliptical integral of second type.
When a semi-elliptical crack in a semi-infinite volume under uniform stress is considered, the configuration
function Y is given by
Y ¼ ks  Y b ð16Þ
with the free surface correction factor ks (see Appendix B):
2
ks ð/; a=cÞ ¼ ½1:13  0:09ða=cÞ  ½1 þ 0:1ð1  sinð/ÞÞ  ð17Þ
Therefore, the procedure proposed in this work for dealing with surface cracks in thick plates that can be
assimilated to a semi-infinite domain consists in solving the equivalent embedded crack (see Section 3.2.1)
to compute the factor Y Cb (Eq. (12)) and using the correction factor ks (Eq. (17)) to account for the effect
of the free surface. Then, the configuration function Y at a point in the crack front is expressed by
Y ¼ ks ð/; a=cÞ  Y Cb ð18Þ
Eq. (18) allows finding K for regular and irregular geometries; by defining the parameter / associated to a
point Q 0 of the crack front, as it is shown in Fig. 9. Point C in the surface corresponds to / = 0°, and for
the deepest point A is / = 90°. Since the surface correction factor ks is a decreasing function of /, the correc-
tion is larger for point C than for point A, as expected. Note that in the case of an elliptical crack, the defi-
nition given by the angular parameter coincides with the one given in Fig. 2.
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2673

For cases that involve finite geometries, Newman and Raju introduced auxiliary functions FS and H which
contemplates the proximity of borders. For a semi-elliptical crack of aspect ratio a/c located amid a finite plate
of width 2b and thickness t (see Fig. 18(a) and (b) in Appendix B), the K solution for a linear normal stress
field is expressed by
rffiffiffiffiffiffi
pa
K ¼ ðrm þ H  rb Þ F S ða=c; a=t; c=b; /Þ ð19Þ
Q
where rm and rb are the remote membrane and pure bending stresses, respectively. Functions H and FS both
depend on the geometry parameters a/c, a/t, c/b, and on the angular parameter /. The function FS can be
rewritten:
F S ¼ fS ða=t; a=c; c=b; /Þ  f/ ð/; a=cÞ ð20Þ
where f/(/, a/c) has been defined in Eq. (13) and the function fS is (see Appendix B for details)
  a 2 a4 
fS ¼ M S1 þ M S2 þ M S3  gS  fw ð21Þ
t t
If r0 is the reference stress for a given problem, such as the stress in the most stressed point in the plate thick-
ness, then

rm ¼ tm  r0
rb ¼ tb  r0

The configuration function Y can be expressed in terms of the non-dimensional magnitudes tm and tb as
follows
K f/
Y ¼ pffiffiffiffiffiffi ¼ ðtm þ Htb ÞfS pffiffiffiffi ð22Þ
r0 pa Q

Then, from Eq. (14), it follows

Y ¼ ðtm þ Htb ÞfS Y b ð23Þ

where Yb is the configuration function of the equivalent crack located in an infinite body under unit uniform
stress, and the pre-factors in the equation account for the corrections due to the free surface, loads and the
plate geometry.
The application of the above procedure to the case of arbitrary cracks in finite plates is immediate. By
substituting the factor Yb by the factor Y Cb given by Calsyf for the equivalent crack:
Y ¼ ðtm þ Htb ÞfS Y Cb ð24Þ
In this way, the stress intensity factor K for surface crack of arbitrary shape under remote membrane and pure
bending stresses can be solved by using the previously defined parameters and the Calsyf solution Y Cb :
pffiffiffiffiffiffi
K ¼ r0 paðtm þ Htb ÞfS Y Cb ð25Þ

4. Results and discussion

4.1. Embedded cracks of irregular front

In order to evaluate the performance of the proposed corrective algorithm, two examples consisting of
irregular crack geometries embedded in an infinite body are studied next.
The first example consists in the crack geometry referred as ‘‘curvan’’, which is generated using the
function:
rðhÞ ¼ 1 þ q  cosð4  hÞ ð26Þ
2674 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

where the parameter ‘‘q’’ is related to the ‘‘deformation’’ of a circumference of radius one (note from Eq. (26)
that q = 0 results in a circumference). Fig. 10 illustrates the crack geometry for q = 0.1. The parameter / is the
same defined in Fig. 2. The analytic definition of the crack geometry allows computing the crack front curva-

y
1.2

1.0

0.8
Q'
0.6

0.4

0.2
φ
θ
0.0
-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2
-0.2 x

-0.4

-0.6

-0.8

-1.0

-1.2

Fig. 10. Curvan crack geometry (q = 0.1).

y
1.5

1.0

Q'

0.5

φ
θ
0.0
-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
x

-0.5

-1.0

-1.5

Fig. 11. Irregular crack geometry ‘‘Guitel’’.

 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2675

ture exactly. Extreme values for the curvature are vnjmax = 2.4169 for h = 0° (/ = 0°), and vnjmin = 0.8895
for h = 45° (/ = 35.3483°).
The crack geometry for the second example is illustrated in Fig. 11. This crack geometry will be referred as
‘‘guitel’’ and it results after the combination of the function
rðhÞ ¼ 1 þ q  cosð2  hÞ ð27Þ
with q = 0.4 for p/2 < h 6 p/2 and an ellipse with aspect ratio r = 0.4286 for p/2 < h 6 3p/2. This geometry
allows studying crack geometries with different aspect ratios and positive and negative curvatures. The crack
aspect ratio is r = 0.5092 and the maximum and minimum curvatures are vnjmax = 2.6295 (h = / = 0°) and
vnjmin = 1.844 (h = 90° and / = 25.3769°).
K-results obtained using Calsyf and the original O-integral algorithm were compared to results from BEM
computations for the two examples. K-computations using BEM were performed following the procedure pre-
sented by Cisilino and Aliabadi [25]. For both examples two BEM models with different discretizations were
carried out in order to verify the independence of the K-results with the model discretization. They are referred

(a)

(b)

Fig. 12. Model discretizations of Curvan: (a) Calsyf (upper right quarter), (b) BEM model (fine mesh).
2676 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

in what follows as ‘‘coarse’’ and ‘‘fine’’ meshes. Error of the BEM solution is estimated in a few percent.
Model discretizations are depicted in Fig. 12.
Results for the Curvan crack (see Fig. 10) are reported in Fig. 13. Results in Fig. 13(a) corresponds to the
loading case r = r0 (uniform stress), while Fig. 13(b) illustrates the results for the loading case r = r0y (pure

1.5
1.2
BEM (fine mesh) 1.0 q = 0.1
BEM (coarse mesh) Q'
1.4 0.8
Oore
this work 0.6

y
1.3 0.4

0.2
φ
K / σw πa

1.2 0.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2
x

1.1

1.0

0.9
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
φ [rad]
(a)

K / σw π a
BEM
0.8
Oore
0.7
this work
0.6

0.5

0.4

0.3

0.2

0.1

0.0
-1.5 -1.0 -0.5 0.0 0.5 1.0 1. 5
-0.1 φ [rad]
-0.2

-0.3

-0.4

-0.5

-0.6

-0.7

-0.8
(b)

Fig. 13. Normalized K results for the Curvan crack geometry with q = 0.1: (a) tension, (b) bending.
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2677

Table 1
Normalized Oore, Kcorr and BEM results at the point of minimum curvature of the Curvan crack geometry under tension
q Minimum curvature Normalized K (BEM) Normalized K (Oore) Error (%) This work Error (%)
0.1 0.8895 1.2651 1.2261 3.08 1.2810 1.26
0.2 4.2754 1.3981 1.3334 4.63 1.3931 0.36
0.3 10.8776 1.5706 1.4058 10.49 1.4688 6.48

1.5
BEM (fine mesh)
1.4
Oore
1.3 this work

1.2
K / σw π a

1.1

1.0

0.9

0.8

0.7
-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
x

Fig. 14. Normalized K results for Guitel crack geometry with q = 0.4 and uniform stress.

free surface bending tension

C x
1
σ ZZ σ ZZ
1
crack φ
a

t A y y
y
2c

Fig. 15. Semi-elliptical crack under bending and tension. The semi-infinite volume case corresponds to a/t ! 0.

Table 2
Normalized K results for the semi-elliptical cracks
Loading case Tension Pure bending
a/c Point Isida [26] Calsyf e (%) Isida [26] Calsyf e (%)
1.0 A 0.659 0.662 0.46 0.173 0.186 7.5
C 0.745 0.728 2.28 – 0.602 –
0.6 A 0.832 0.845 1.56 0.274a 0.270 1.5
C – 0.714 – – 0.595 –
a
Value for thick plate with a/t = 0.01.

bending). Also included in Fig. 13 are the results obtained using the O-integral and those of the BEM
computations.
2678 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

Results for cracks with q = 0.2 and q = 0.3 are summarized in Table 1. In all cases the errors are calculated
with respect to the BEM results obtained using fine mesh in the point with maximum error.
Similarly, results for Guitel corresponding to the loading case r = r0 (uniform stress), are given in Fig. 14.
These results allow to conclude that the proposed corrective function improves the K results obtained using
O-integral algorithm when dealing with irregular embedded cracks.

4.2. Surface cracks

The performance of the proposed algorithm is demonstrated in this section for surface cracks in thick and
thin plates under tension and bending. The studied geometries are limited to the case of semi-elliptical cracks,
as to the author’s knowledge, these are the only results available in the literature. The first example consists in
a semi-elliptical crack in a semi-infinite volume (see Fig. 15). Two crack aspect ratios are considered, a/c = 1

1.0

0.9
SURFACE POINT (C)
0.8

0.7

0.6 DEEPEST POINT (A)

a

Calsyf (A)
0.5
K/

Calsyf (C)
0.4 Shiratori et al. [28]
0.3 Newman-Raju [23, 24]
Isida et al. [26]
0.2
Glinka [29]
0.1

0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
(a) a/t

1.4
1.3 Calsyf (A)
1.2 Calsyf (C)
1.1 Shiratori et al. [28]
1.0 Newman-Raju [23, 24]
0.9 Isida et al. [26]
0.8 Glinka [29]
a

0.7
K/

0.6
SURFACE POINT (C)
0.5
0.4
0.3 DEEPEST POINT (A)
0.2
0.1
0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
(b) a/t

Fig. 16. Normalized K results for the semi-circumference (aspect ratio a/c = 1) in a finite plate: (a) tension, (b) bending.
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2679

(semi-circular crack) and a/c = 0.6. Both geometries are solved subjected to tension and pure bending stress
fields. Computed results are reported in Table 2 together with those due to Isida [26] tabulated by Murakami
[27]. Results are reported and comparedp for
ffiffi two positions: points A (/ = 90°) and C (/ = 0°). In all cases
results are normalized with respect to r0 pa.
Results computed using Calsyf show an excellent agreement with those of the reference with differences
around 2%. The only exception is point A for the semi-circular crack in bending. For this case the difference
is 7.5%.
In the second example the effect of finite-thickness is considered. The same crack geometries and loading
cases of the first example are studied for plate thickness a/t = 0.01, 0.2, 0.4, 0.6 and 0.8. Results for the
semi-circular crack are reported in Fig. 16, while the results for the semi-elliptical crack are given in
Fig. 17. Results due to Shiratori [28], Isida [26], Glinka [29] and Newman and Raju [23,24] are used for com-
parison. In all cases the reference results report an accuracy of a few percent.

1.2
1.1
DEEPEST POINT (A)
1.0 (OPEN SYMBOLS)
0.9
0.8
0.7
a

SURFACE POINT (C)

0.6
(FILLED SYMBOLS) Calsyf, (A)
K/

0.5 Calsyf, (C)

0.4 Shiratori et al. [28]
0.3 Newman-Raju [23, 24]
0.2 Glinka [29]
Isida et al. [26]
0.1
0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
a/t
(a)

1.4

1.3
Calsyf (A)
1.2 Calsyf (C)
Shiratori et al. [28]
1.1
Newman-Raju [23, 24]
1.0 Glinka [29]
0.9 Isida et al. [26]
a

0.8
K/

0.7

0.6 SURFACE POINT (C)

(FILLED SYMBOLS)
0.5
DEEPEST POINT (A)
0.4 (OPEN SYMBOLS)
0.3

0.2
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
a/t
(b)

Fig. 17. Normalized K results for the semi-ellipse of aspect ratio 0.6 in finite plate: (a) tension, (b) bending.
2680 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

Results reported in Figs. 16 and 17 show a very good agreement between the obtained results and those
reported in the literature. In particular obtained results are almost coincident to those from Newman and raju
[23,24], from which the correction functions for the free surface and finite thickness effects were taken.

5. Conclusions

A numerical tool for the assessment of mode I embedded and surface cracks has been presented in this
paper. The proposed tool is an extension of the O-integral algorithm due to Oore and Burns for embedded
cracks, which has been effectively extended to deal with cracks of arbitrary shape and surface cracks.
Irregular crack fronts are solved by introducing a corrective function which depends on the crack aspect
ratio and the local crack front curvature. Solutions due to Newman and Raju are used to account for the effect
of free surfaces and finite thickness. The accuracy of the proposed tool is assessed by solving a number of
examples and the results compared to BEM computations and results available in the literature. Obtained
results compare well to the reference solutions.
The devised tool constitutes a low-cost and versatile means for computing stress intensity factors for cracks
of arbitrary shape. It delivers results within a few percent error while when compared to FEM and BEM mod-
els it demands a relatively low computing cost and data preparation.
This work will be extended and customized for the analysis of partially closed cracks and fatigue crack
propagation.

Acknowledgements

This work was funded by grants from the Agencia Nacional de Promoción Cientı́fica y Técnica de la
República Argentina (PICT 12-12528) and the University of Mar del Plata (ING125/03).

Appendix A. Embedded elliptical crack

An empirical K equation for an embedded elliptical crack in a finite plate subjected to remote tension r,
shown in Fig. 18 (see parameter / in Fig. 2), was obtained by Newman and Raju [23]. The equation is
rffiffiffiffiffiffiffiffiffi 
pa a a c 
KI ¼ r  F ; ; ;/ ðA1Þ
Q c t b

where
c1:65
Q ¼ 1 þ 1:464 ðA2Þ
a

2b

t
a
2t

2c

Fig. 18. Embedded elliptical crack configuration.

 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2681

and
 a2 a4 
F ¼ M1 þ M2 þ M3  g  fw  f/ ðA3Þ
t t
The function f/ was taken from the exact solution for an embedded elliptical crack in an infinite solid [18], fw is
a finite-width correction factor, and the function g is a fine-tuning curve-fitting function. For a/c 6 1, the func-
tions are expressed by
M1 ¼ 1 ðA4Þ
0:05
M2 ¼  3=2 ðA5Þ
0:11 þ ac
0:29
M3 ¼  3=2 ðA6Þ
0:23 þ ac
a4
t  
g ¼1  j cos /j ðA7Þ
1 þ 4 ac
" #
pcraffiffiffi 1=2
fw ¼ sec ðA8Þ
2b t
  1=4
a 2 2 2
f/ ¼ cos / þ sin / ðA9Þ
c
From the above equations, the limiting case of the elliptical crack embedded in an infinite volume under re-
mote tension (i.e., when a/t ! 0) is expressed by
KI f/
Y b ¼ pffiffiffiffiffiffi ¼ pffiffiffiffi ðA10Þ
r pa Q

Appendix B. Semi-elliptical surface crack in a finite width plate

The case of a semi-elliptical surface crack in a finite plate subjected to tension and bending was solved in
another work of Newman and Raju [24] (see Fig. 19(a) and (b) for details). The K equation for combined ten-
sion and bending loads is
rffiffiffiffiffiffiffiffiffi a a c 
pa
K I ¼ ðrm þ H  rb Þ  FS ; ; ;/ ðB1Þ
Q c t b
where rm represents the remote uniform-tension stress, and rb the remote outer-fiber bending stress. The Q
factor was given in Appendix A, and the function FS has a similar expression:
 a2 a4 
F S ¼ M S1 þ M S2 þ M S3  gS  fw  f/ ðB2Þ
t t
The functions fw and f/ are the functions already given in Appendix A. The other functions in the above equa-
tion are expressed by
 a
M S1 ¼ 1:13  0:09 ðB3Þ
c
0:89
M S2 ¼ 0:54 þ  ðB4Þ
0:2 þ ac
1:0  a24
M S3 ¼ 0:5  a þ 14 1:0  ðB5Þ
0:65 þ c c
  a 2 
2
gS ¼ 1 þ 0:1 þ 0:35  ð1  sin /Þ ðB6Þ
t
2682 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

2h a

2c

t
2b
(a)

Sm M

3M
2h Sb =
bt2
2c a

2b
t

Sm M
(b)

Fig. 19. (a) Elliptical surface crack in a finite plate, (b) surface-cracked plate subjected to tension or bending loads.

The function H has the form

H ¼ H 1 þ ðH 2  H 1 Þ sinp / ðB7Þ

where
a a
p ¼ 0:2 þ þ 0:6 ðB8Þ
c t
a a  a
H 1 ¼ 1  0:34  0:11 ðB9Þ
t c t
a a2
H 2 ¼ 1 þ G1 þ G2 ðB10Þ
t t
 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684 2683

In this equation for H2

a
G1 ¼ 1:22  0:12 ðB11Þ
c
a0:75 a1:5
G2 ¼ 0:55  1:05 þ 0:47 ðB12Þ
c c
The limit case of a semi-elliptical surface crack in a semi-infinite volume under uniform remote tension gives
(i.e., rb = 0), from the above equations:
rffiffiffiffiffiffiffiffiffi  
pa
K I ¼ rm  ½1:13  0:09ða=cÞ 1 þ 0:1ð1  sin /Þ2  f/ ðB13Þ
Q

taking account of Eq. (A10) in Appendix A, the above equation can be expressed by
KI
Y ¼ pffiffiffiffiffiffi ¼ ks  Y b ðB14Þ
rm pa
with
 
2
ks ¼ ½1:13  0:09ða=cÞ 1 þ 0:1ð1  sin /Þ ðB15Þ

References

[1] American Society of Mechanical Engineers (ASME) Boiler and Pressure Vessel Code. Section IX: Rules for Inservice Inspection of
Nuclear Power Plant Components. New York: American Society of Mechanical Engineers, September 2004.
[2] PD 6493: 1991. Guidance on some methods for the derivation of acceptance levels for defects in fusion welded joints, British
Standards Institution, August 1991.
[3] Harrison RP, Loosemore K, Milne I, Dowling AR. Assessment of the integrity of structures containing defects. CEGB R6, Rev 2
1980.
[4] Kumar V, German MD, Shih CF. An engineering approach for elastic-plastic fracture analysis. EPRI Report NP-1931, Palo Alto,
CA: Electric Power Research Institute, 1981.
[5] EPRI Report GS-6724. Condition assessment guidelines for fossil fuel power plant components. USA: Electric Power Research
Institute, 1990.
[6] BS 7910. Guidance on methods for assessing the acceptability of flaws in metallic structures. London, UK: British Standard
Institution, 1999.
[7] Pisarski H, Wallin K. The SINTAP fracture toughness estimation procedure. The Welding Institute; 2000.
[8] API recommended Practice 579, Fitness for Service, 1st ed., American Petroleum Institute, August 2000.
[9] Murakami Y, Namat-Nasser S. Growth and stability of surface flaws of arbitrary shape. Engng Fract Mech 1983;17(3):193–210.
[10] Murakami Y. Analysis of stress intensity factors of modes I, II and III for inclined surface cracks of arbitrary shape. Engng Fract
Mech 1985;22(1):101–14.
[11] Noda N-A. Stress intensity formulas for three-dimensional cracks in homogeneous and bonded dissimilar materials. Engng Fract
Mech 2004;71:1–15.
[12] Aliabadi MH, Rooke DP. Numerical fracture mechanics. Kluwer Academic Publishers; 1991.
[13] Bueckner HF. A novel principle for the computation of stress intensity factors. Z Angew Math Mech 1970;50:529.
[14] Rice JR. Some remarks on elastic crack-tip stress fields. Int J Solid Struct 1972;8:751.
[15] Wu Xue-Ren, Carlsson AJ. Weight functions and stress intensity factor solutions. 1st ed. Pergamon Press; 1991.
[16] Oore M, Burns DJ. Estimation of stress intensity factors for embedded irregular cracks subjected to arbitrary normal stress fields.
Trans ASME 1980;102(May):202–11.
[17] Desjardins JL, Burnsy DJ, Thompson JC. A weight function technique for estimating stress intensity factors for cracks in high
pressure vessels. J Pressure Vessel Technol 1991;113(1):55–64.
[18] Irwin GR. The crack extension force for a part-through crack in a plate. ASME J Appl Mech 1962:651–4.
[19] Shah RC, Kobayashi AS. Stress intensity factor for an elliptical crack under arbitrary normal loading. Engng Fract Mech
1971;3:71–96.
[20] Isida M, Yoshida T, Noguchi H. Prelim. In: Proc Japan Soc Mech Engrs and Japan Soc Precision Engng, Mie District, No. 823-3,
1982. p. 15–7.
[21] Oore M, Burns DJ. Estimation of stress intensity factors for irregular cracks subjected to arbitrary normal stress fields. In: Proceeding
of the 4th international conference on pressure vessel technology, London, I Mech. E, vol. 1, 1980. p. 139–47.
[22] Grueter L, Huget W, Kylla H. Weight functions and stress intensity magnification factors for elliptical and semi-elliptical cracks
under variable normal stress. Interatom Report INTAT No. 32.05095.0, May 1981.
2684 H.L. Montenegro et al. / Engineering Fracture Mechanics 73 (2006) 2662–2684

[23] Newman JC, Raju IS. An empirical stress intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1):185–92.
[24] Newman JC, Raju IS. Stress intensity factor equations for cracks in three-dimensional finite bodies. Fracture mechanics: 14th
symposium – theory and analysis. In: Lewis JC, Sines G, editors. ASTM STP 791, vol. I. American Society for Testing and Materials;
1983. p. I-238–65.
[25] Cisilino AP, Aliabadi MH. Three-dimensional BEM analysis for fatigue crack growth in welded components. Int J Pressure Vessels
Piping 1997;70:135–44.
[26] Isida M, Noguchi H, Yoshida Y. Tension and bending of finite thickness plates with a semi-elliptical surface crack. Int J Fract
1984;26:157–88.
[27] Murakami Y. Stress intensity factors handbook. The Society of Materials Science. Japan: Pergamon Press; 1990.
[28] Shiratori M, Miyoshi T, Tanikawa K. Analysis of stress intensity factors for surface crack subjected to arbitrarily distributed surface
stresses (2nd report, Analysis and application of influence coefficients for flat plates with a semi-elliptical surface crack). Trans Jpn
Soc Mech Engrs 1986;52(474):390–8.
[29] Glinka G. AFGROW program, ver. 4.0004.12.10. AFGROW web site, public service by the Analytical Structural Mechanics Brand,
Air Vehicles Directorate, US Air Force Research Laboratory.