A complex-variable virtual crack extension finite element method for elastic-plastic fracture mechanics

Accepted Manuscript A Complex-Variable Virtual Crack Extension Finite Element Method for Elastic-Plastic Fracture Mechanics Arturo Montoya, Daniel Ramirez-Tamayo, Harry Millwater, Matthew Kirby PII: DOI: Reference:

S0013-7944(18)30677-5 https://doi.org/10.1016/j.engfracmech.2018.09.023 EFM 6158

To appear in:

Engineering Fracture Mechanics

Received Date: Accepted Date:

5 July 2018 17 September 2018

Please cite this article as: Montoya, A., Ramirez-Tamayo, D., Millwater, H., Kirby, M., A Complex-Variable Virtual Crack Extension Finite Element Method for Elastic-Plastic Fracture Mechanics, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/j.engfracmech.2018.09.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Complex-Variable Virtual Crack Extension Finite Element Method for Elastic-Plastic Fracture Mechanics Arturo Montoya1 , Daniel Ramirez-Tamayo2 , Harry Millwater3 , Matthew Kirby4 1

Associate Professor, Department of Civil and Environmental Engineering, The University of Texas at San Antonio, San Antonio, TX, 78249 (corresponding author). E-mail: [email protected]

2

Graduate Student, Department of Mechanical Engineering, The University of Texas at San Antonio, San Antonio, TX, 78249 3

Professor, Department of Mechanical Engineering, The University of Texas at San Antonio, San Antonio, TX, 78249

4

Graduate Student, Department of Mechanical Engineering, The University of Texas at San Antonio, San Antonio, TX, 78249

Abstract The virtual crack extension (VCE) approach for computing the energy release rate using a complex-variable finite element method (ZFEM) is extended to nonlinear materials undergoing plastic deformation. The method consists of performing a numerical derivative of the potential energy with respect to a crack extension via the complex Taylor series expansion method (CTSE). ZFEM does not depend on the existence of a strain energy density and does not require contour or domain integrals. As a result, ZFEM overcomes the contour path dependency and monotonic loading limitations of the Jintegral in elastic-plastic fracture problems. The nonlinear fracture analysis capability is demonstrated using material models based on the deformation and incremental theory of plasticity and using a loading/unloading cycle. Under monotonic loading, the ZFEM method remains robust at all levels of plasticity and its results are of the same accuracy as those provided by the ASTM standards and the saturated J-integral value. Moreover,

1

when unloading the elasto-plastic material, the J-integral contours provide erroneous energy release rate values, whereas ZFEM results are in excellent agreement with those obtained by finite differencing.

Keywords: nonlinear fracture, energy release rate, virtual crack extension method, complex variable finite element method, elasto-plastic fracture

1

Introduction

Elastic-plastic fracture mechanics (EPFM) evaluates the crack growth potential in materials that exhibit nonlinear elastic behavior or plastic deformation during fracture. The crack growth assessment is based on calculating the energy release rate (ERR) of the structure [1]. The ERR is a function of the stress-strain behavior of the material, which can be estimated through a nonlinear finite element analysis (FEA) at several load increments. Hellen and Blackburn [2] conducted a comprehensive review on the numerical techniques that can be used to calculate the ERR within a nonlinear FEA. They classified the methods into two groups, contour integrals and virtual crack extension (VCE) techniques. These techniques show numerical issues during loading, when the plastic zone size is larger than the singularity zone, and after unloading from a yield state. Numerical methods based on contour integrals provide accurate ERR estimates for linear elastic fracture mechanics (LEFM) problems. Contour integral evaluations remain robust for nonlinear structural problems as long as the material is assumed to be homogeneous and exhibits an elastic response to unloading, e.g. deformation theory of plasticity [3]. Under elastic conditions, contour integrals are known to be path independent, i.e., the ERR estimate does not change by varying the path of integration. The path independent J-integral proposed by Rice [4] is the most popular contour-based method. The J-integral is typically implemented in FEA codes by means of a domain-based approach [5]. The integration is performed with respect to the surface area (2D crack) or volume (3D crack) enclosed by the contour and evaluated by specifying a mathematical function 2

that can be interpreted as imposing a normalized virtual displacement in the direction of crack growth. In FEA practice, contour integrals are evaluated well away from the crack tip in order to reduce possible numerical inaccuracies at the crack tip. The J-integral becomes path dependent when modeling irreversible plastic deformation. The calculated J-values are not constant but increase monotonically as the size of the contour increases. Nonetheless, the J-integral has been employed to analyze cracks in elastic-plastic materials under monotonic loading since numerical studies have shown that the J-integral approaches a saturated, accurate value as the integration contour encloses the entire plastic zone [6]. However, the path is restricted to begin and end on a traction free surface at the crack plane [7]. Non-trivial correction terms may be appended to the standard J-integral to recover its path independent property. These correction terms are problem specific; for example, modifications have been proposed for nonhomogeneous materials [8, 9], functionally graded materials [10, 11], large deformation [12], loading/unloading [13], shells [14], residual stresses [15, 16] and curved cracks [17]. Not all of these corrections are available in commercial finite element software. Another restriction of the J-integral is that its evaluation depends on the existence of a strain energy density from which stresses can be uniquely derived. This assumption is valid for deformation plasticity, but not for modeling history-dependent elastic-plastic materials [18]. In practice, an equivalent elastic function that accounts for the elastic strain energy and the plastic dissipation has been used to evaluate the ERR in elastic plastic fracture problems under monotonic loading conditions [19]. However, no unique function can be adopted to evaluate the ERR during unloading post yielding. In LEFM, the J-integral represents the driving force causing the crack to propagate. On the other hand, for EPFM, the J-integral is a measure of the intensity of the stress field at the crack tip [20, 21]. In elasticity, the energy released is used to form a new crack surface area. However, in EPFM, crack progression is a two stage process. First, plasticity damages the material ahead of the crack tip and then elasticity causes a microstability that propagates the crack [22]. Thus, the J-integral can be used to measure 3

the singular stress at the crack tip for stationary cracks, but loses its physical meaning for growing cracks [23]. In addition, during crack propagation, the monotonic loading assumption of the J-integral is not met as there are small regions where the material is unloaded [24]. Virtual crack extension techniques for evaluating the ERR involve determining the energy differences in the system with respect to a small crack tip extension. These can be performed by obtaining the difference in energies from two separate finite element analyses with the crack slightly extended in the second analysis or by using virtual crack extension (VCE) methods. The former, also known as finite differencing, has the advantage of being applicable for any plastic zone size and for cases involving loading and unloading. The disadvantage of this approach is that it is subject to subtraction cancellation errors and its accuracy is dependent on the size of the crack extension. The latter approach is preferred under monotonic loading as it reduces the computational time as compared to at least two nonlinear analyses and avoids truncation errors. VCE consists of calculating the changes in finite element variables that result from shifting nodal coordinates surrounding the crack tip in the direction of the crack advance. One example is the stiffness derivative finite element approach presented by Parks for linear [25] and nonlinear [26] applications. Hellen [19] evaluated the variation in the potential energy with respect to a VCE by transforming it into a volume integration over the elements surrounding the crack tip. DeLorenzi [27, 28] used the VCE scheme and continuum theory to derive an analytical expression for the ERR which can be used with any numerical technique, including FEA. Similar to the J-integral, the use of the current state of the art VCE methods in incremental elastic plastic materials is limited to small plastic zone sizes and monotonic loading conditions as their evaluation depends on a unique strain energy density function. Recently, a new approach based on implementing the complex Taylor series expansion (CTSE) [29] within a complex variable finite element code has emerged as an alternative VCE method to compute the ERR [30]. The complex finite element method, ZFEM, 4

numerically differentiates the finite element system of equations with respect to an extension of the crack length or crack area using CTSE. This method is analogous to the finite differencing method, but computes the difference in the strain energy during a single computational run. Most importantly, ZFEM calculates derivatives without subtracting similar numbers, thereby avoiding subtraction cancellation errors and circumventing the determination of the step size, as long as a sufficiently small step size is used. ZFEM has shown to have the same accuracy as the J-integral when computing the ERR for 2D and 3D linear elastic models [30]. Ramirez Tamayo et al. [31] demonstrated the use of ZFEM in thermal fracture problems. Moreover, ZFEM has been applied to compute very accurate weight functions by evaluating the derivative of the crack opening displacement with respect to crack length [32, 33, 34]. This paper extends the complex-variable VCE presented by Millwater et al. [30] for LEFM to EPFM. The computation of the ERR for non-linear materials requires to evaluate the energy difference between the systems with the original and perturbed cracks throughout the incremental solution procedure [19]. Due to the complex variable nature of ZFEM, this computation can be processed while solving the FEA system of equations. The real domain handles calculations related to the original crack position, while the imaginary domain calculations are based on the extended crack position. Montoya et al. [35] calculated sensitivities for bodies with nonlinear elastic and elastic-plastic materials using ZFEM, which was shown successful in providing accurate partial derivatives of stress dependent variables with respect to geometric, material, and load parameters. In this research, ZFEM computes the ERR for nonlinear fracture applications from the imaginary part of the global complex-variable strain energy that is accumulated over the loading increments. Due to its similarity to finite differencing method, ZFEM estimates are independent of the plastic zone size and loading conditions. The paper is organized as follows. First, a background on the nonlinear ERR calculations is provided. Then, the ZFEM approach for obtaining the ERR within an incremental-iterative finite element procedure is presented. In the numerical results sec5

tion, the approach is verified by comparing the ZFEM results to those obtained via the domain-based J-integral calculation and the ASTM approach for a compact tension specimen with Ramberg Osgood (deformation plasticity theory) and bilinear hardening (incremental plasticity theory) material models. The ZFEM capability of computing the nonlinear ERR during unloading is demonstrated by modeling a compact tension specimen and edge cracked plate under one loading/unloading cycle and verifying the results against finite differencing results. Lastly, the relevance of introducing ZFEM for nonlinear fracture mechanics is discussed, followed by concluding remarks.

2

Background

The definition of the energy release rate is the same for linear and nonlinear materials, except that G is replaced by J for nonlinear materials [1] , J =−

dΠ dA

(1)

where Π is the potential energy and A is the crack area. The potential energy for a linear or nonlinear material is defined as, (2)

Π=U −F

where U is the strain energy stored in a body and F is the work done by external forces. In a deformed linear or non-linear material, Π can be expressed as Z Z Z Π= W dV − Qi ui dV − Ti ui dS V

V

(3)

ST

where W is the strain energy density in a body, Qi and Ti are the prescribed body and surface forces, respectively, and ui is the displacement associated with these prescribed forces. V is the volume of the body and ST is surface area under mechanical loading. Rice [4] showed that J can be obtained by performing a line integral over a closed contour Γ as Z ∂ui J = W dx2 − Ti dS ∂x1 Γ

(4)

xi represents the rectangular coordinates with the x2 direction taken normal to the crack 6

line. The relationship between strain energy density and stresses, σij , is σij =

∂W ∂ij

(5)

where ij are the total strain components in the material. Shih et al. [5] presented a domain integral approach to facilitate the evaluation of the contour line integral. This integral requires the incorporation of a function, denoted commonly as “q”, that is used to evaluate the integral numerically. This function can be interpreted as imposing a normalized virtual displacement in the direction of crack growth [1]. For a linear or nonlinear elastic material under quasistatic conditions and in the absence of thermal strains, body forces and crack face tractions, Equation 4 becomes:  Z  ∂ui ∂q J= σij − W δ1i dA ∂x1 ∂xi A

(6)

where A is the area enclosed by Γ. Hellen [19] evaluated Equation 1 by implementing the VCE technique in a nonlinear FEA code. In the nonlinear finite element method, the equivalent expression to Equation 3 at a load level λ is given as: Π(λ) = U (λ) − uTλ Fλ

(7)

where the strain energy scalar U (λ) and the vectors uλ and Fλ correspond to the load and displacement for that particular load level. The strain energy represents the summation of the strain energy of all the elements, nel. Thus Π(λ) =

nel X j=1

Uj (λ) − uTλ Fλ

(8)

Under a VCE scheme, the approach consists of finding the variation in potential energy δΠ associated with a virtual crack extension δa. Nodes on and within an interior region surrounding the crack tip are virtually extended by δa. Hellen [19] showed that for two dimensional problems, dΠ 1 dΠ 1 δΠ = = lim dA t da t δa→0 δa

(9)

7

where t is the thickness of the structure and δΠ/δa simplifies to δΠ X ∂Uj (λ) = δa ∂a j

Equation 10 was transformed to the following equivalent integral, ! Z ˆ ∂Uj ∂W ˆ ∂|J| = |J| + W n dV ∂xni ∂xni ∂xi V0

(10)

(11)

ˆ is the determinant of the Jacobian transformation matrix, xn stands for the where |J| i coordinate system component at node number n, and V0 is the volume in natural coordinates. Note that the J-integral formulation (Equation 6) and the VCE formulation (Equation 11) depend on the existence of a unique expression for the strain energy density W throughout the loading cycle.

3

Methodology

The Complex Taylor Series Expansion (CTSE) method is a first-order 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 without difference operations; hence, the step size used in CTSE can be smaller in magnitude as compared to the finite differencing step size. The derivatives are computed using a small perturbation along the imaginary axes. That is, variable X = xo is perturbed to X = xo + ih, where i denotes an imaginary number and h denotes the step size. The formulae for the derivatives can be derived from the Taylor series representation of the function evaluated at the complex sample point. f (xo + hi) = f (xo ) + f (1) (xo )

hi (hi)2 (hi)3 + f (2) (xo ) + f (3) (xo ) + ... 1! 2! 3!

(12)

where f (1) denotes the first derivative, f (2) the second, etc. The first order derivative can be obtained by taking the imaginary part of both sides of Eq. 12, solving for f (1) , and ignoring higher order terms.   Im f (xo + hi) (1) f (xo ) ≈ h

(13)

8

Im[ ] denotes the imaginary component. The step size can be made arbitrarily small with no concern about round-off error. This method can be applied to evaluate the change in the strain energy of a body due to a step size h = ∆A, with respect to the crack surface,   Im U ∗ (A + ∆Ai) ∂U ≈ (14) ∂A ∆A where the superscript ∗ denotes a complex variable. In an FEA code, the strain energy of the body is the sum of the strain energy evaluated at nel elements. Thus, h i ∗ nel Im U (A + ∆Ai) X j ∂U ≈ ∂A ∆A j=1

(15)

For 2D problems, ∆A = t∆a, where t is the thickness and ∆a represents the magnitude of ~ in the direction of crack growth. Hence, Equation 15 simplifies the translation vector, ∆a, to

h i nel Im U ∗ (a + ∆ai) X j 1

∂Uj ≈ ∂A t

j=1

(16)

∆a

The evaluation of the ERR using CTSE requires: a) a complex-variable finite element method that allows perturbations in the imaginary axes, and b) the evaluation of the complex-variable strain energy at the element level. These two items are discussed further in the following subsections.

3.1

Complex Finite Element Method

The complex variable finite element method, ZFEM, is a complex representation of the traditional finite element method; i.e., nodal coordinates, finite element variables, and output variables become complex (having a real and imaginary part). Derivatives with respect to the crack length (or area) can be obtained by perturbing the crack length (2D) or area (3D) through the imaginary nodal coordinates. For example, Figure 1 shows a crack tip region modeled with a spider web mesh, which is considered the best practice approach for 2D fracture problems. In ZFEM, a virtual crack extension is implemented by elongating the nodes within close proximity to the crack tip (shadowed region) along ~ the crack’s self similar direction (indicated by the translation vector, ∆a). The magnitude 9

for the virtual increment is typically 10−10 multiplied by the size of the smallest element at the crack tip [30, 35].

~ ∆a θ

Figure 1: Schematic of the perturbation of the crack length along the self~ similar direction, indicated by the translation vector, ∆a. The value of the imaginary part of the complex nodal coordinates is defined by the components of the perturbation vector, x∗n = xn + ∆ax i yn∗

(17)

= yn + ∆ay i

where ∆ax = ∆a cos θ and ∆ay = ∆a sin θ. Hence, for a horizontal crack, x∗n = xn + ∆ai √ √ and yn∗ = yn , and for a 45◦ angled crack x∗n = xi + ∆a/ 2i and yn∗ = yn + ∆a/ 2i. Although there is no formal proof, empirical evidence suggests that the size of the perturbed region be selected as the ring of 4 or 5 elements surrounding the crack tip. At the completion of the analysis, the real part of a complex output is the traditional finite element output, while the imaginary part contains the sensitivity information of the variables with respect to the crack extension.

3.2 3.2.1

Energy and Energy Release Rate Calculation Linear Materials

The finite element notation of strain energy for an element j with linear elastic material is: 1 Uj = uTj Kj uj 2

(18)

10

where uj is the displacement vector and Kj is the element stiffness matrix. In ZFEM, all of these FEA variables become complex. Thus, the ERR for an element with a linear material, denoted as G, is obtained as Im[Uj∗ ] Gj = − ∆A (19) h i 1 T∗ ∗ ∗ =− Im uj Kj uj 2∆A The ERR of the body is then obtained by summing the imaginary component of all the elements in the model, refer to Equation 16.

3.2.2

Nonlinear Materials

For nonlinear applications, the strain energy formulation no longer has the quadratic form of Equation 18, but it is accumulated through the incremental-iterative solution procedure. Here below, two approaches are presented for calculating the strain energy and the ERR, the strain energy density approach and the force-displacement approach. The latter approach does not require the existence of the strain energy density to evaluate the ERR.

3.2.2.1

Strain Energy Density Approach The strain energy, U , is accumulated

during the loading increments and can be obtained by integrating the strain energy density, W , over the volume of the element. The complex strain energy at the element level can be evaluated by employing the traditional Gauss-quadrature approach. For a two dimensional isoparametric element, Z ngp mgp X X ∗ ∗ Uj = W dV = wk wl W ∗ (ξk , ηl )t |Jˆ∗ (ξk , ηl )| Vj

(20)

k=1 l=1

∗

ˆ is the complex determinant of the Jacobian, the pair (ξk , ηl ) represents an integration |J| point in the natural coordinates, and wk and wl are weights of the integration scheme. For deformation plasticity theory, W is well defined. For a Ramberg Osgood material model, W =

1 + ν 2 3 1 − 2ν 2 n α σe + p + σen+1 n−1 3E 2 E n + 1 Eσo 11

(21)

E is the elastic modulus, ν is the Poisson’s ratio, n is the hardening exponent, α is the yield offset at a stress of σo , σe is the Mises equivalent stress, and p is the hydrostatic stress. In ZFEM, the state dependent variables (SDVs) of the strain energy density formulation become complex when the crack is extended in the imaginary direction. Hence, σe and p will become complex and Equation 21 will yield a complex valued W . For incremental plasticity, W can be treated as nonlinear elastic under the condition of monotonic loading. Hence, a strain energy density for an equivalent elastic material can be defined as the sum of the elastic strain energy density, W el , plus the plastic dissipated energy per unit volume, W pl , el

W =W +W

pl

1 = σ : el + 2

Z

t

σ : ˙pl dt

(22)

0

el is the elastic strain tensor, ˙pl is the plastic strain rate tensor, σ is the stress tensor and “:” is the double contracted product. Following the complex representation of Equation 22, the ERR can be obtained following the aforementioned procedure. However, note that this approach requires an analytical expression for W and is only suitable for monotonic loading of elastic-plastic materials as a unique W cannot be defined for history-dependent plasticity.

3.2.2.2

Force-Displacement Approach Alternatively, the complex strain energy

can be obtained without invoking a strain energy density formulation, but using the force and displacement vectors computed during an FEA. This procedure is analogous to the experimental tests performed by Landes and Bangley [36, 37], in which the nonlinear strain energy was obtained from the area under the load-displacement curves. In nonlinear ZFEM, the internal force vector at an element level is computed as Z ∗ fj = BT ∗ σ ∗ dV

(23)

Vj

where B∗ is the complex strain-displacement matrix and σ ∗ is the complex stress vector. A numerical integration scheme can be employed to evaluate the strain energy of each

12

element. The simplest scheme is the trapezoidal rule Z u Nsteps −1 1 X ∗ ∗ ∗ ∗ ∗ Uj = fj (u )du ≈ (u∗k+1 − u∗k ) · (fk+1 − fk∗ ) 2 u0 k=1

(24)

where u∗k and fk∗ are the complex-variable displacement and force vectors at loading step k, respectively, Nsteps represents the number of loading steps and “·” is the dot product. The complex displacement vector increment is solved at every load increment by the Newton procedure. This approach requires that the complex internal force vector from the previous step k becomes available in increment k + 1. Since u∗k+1 − u∗k varies throughout the loading increments, other numerical integration schemes for unequally spaced data may be considered to increase the accuracy of the integral evaluation. For instance, adapting Simpson’s rule for unequally spaced load increments in FEA follows the form,  u∗ − u∗k+1  ∗  (u∗k+2 − u∗k )2 1 h ∗ 2 − k+2 · f + 2 − · fk+1 k ∗ ∗ ∗ ∗ ∗ ∗ 6 u − u (u − u )(u − u ) k+1 k k+1 k k+2 k+1 k=1 (25)  i ∗ ∗  uk+1 − uk ∗ + 2− ∗ · fk+2 uk+2 − u∗k+1 Nsteps −2

Uj∗

4

≈

X

Implementation in a Commercial Software Program

The nonlinear fracture mechanics feature of ZFEM can be implemented in Abaqus through a complex-variable user element (UEL) routine, as described by Montoya et al. [35]. Although complex variable operations are performed within the UEL, the complex system of equations is solved using the Abaqus incremental-iterative method for a real valued problems as shown in Figure 2. Duplicate nodes are defined in order to apply the perturbations in the imaginary domain (see Figure 3). The UEL assembles the complex coordinate variables by appending the value of the corresponding duplicate node in the Im imaginary term, e.g. x∗n = xRe n + xn i. The value of the duplicate node is defined by the

13

perturbation vector, as specified in Equation 17.

Begin Analysis Perturb the nodes sourounding the crack tip (∈ S) in the crack growth direction Define Initial Conditions Perturbation vector: ~ = {∆ax , ∆ay , ∆az } ∆a

Start of Step

User Element (UEL) Start of Increment

Read Input: {x1 , y1 , z1 , ..., xn , yn , zn } Complexify Input:

( Im x∗n = xRe n + ixn 0 if ∈ /S ∗ Re Im Im Im Im yn = yn + iyn xn , yn , zn = ~ if ∈ S ∆a Im Re ∗ zn = zn + izn

Start of Iteration

Finite Element Operations (UEL/UMAT):

Begin User Element

Converged?

Uj∗ (λ) = f (u∗1 , v1∗ , w1∗ ..., u∗n , vn∗ , wn∗ ) Decomposition of UEL Variables:   Re[K∗ ] −Im[K∗ ] Stiffness Matrix: K = Im[K∗ ] Re[K∗ ]   Re[f ∗ ] Force Vector: f = ∗ Im[f ]

No

Yes

Strain Energy: Write Output

Uj (λ) =

No

End of Step?

Yes

( R Re[ Vj W ∗ dV ] R Re[ fj∗ (u∗ )du∗ ]

Energy Release Rate:

Jj (λ) =

Im[Uj∗ (λ)] ∆A

Figure 2: Implementation of ZFEM in a Commercial Finite Element Software.

Im

Real

Nodal Coordinates

Figure 3: 16-noded isoparametric element (8 real and 8 imaginary)

14

The complex element stiffness matrix is decomposed into an all-real valued matrix by invoking the Cauchy-Riemann matrix representation of a complex variable. 

∗ Re[Kj ]

Kj = 



−Im[K∗j ]

Im[K∗j ]

Re[K∗j ]

(26)



The complex force vector is also decomposed into an all real valued vector as    Re[f ∗ ] j fj =   Im[f ∗ ] j

  

(27)

 

The Abaqus solver assembles the system of equations and obtains the displacement vector, u, and its derivative with respect to the crack size, a.

u=

  

uRe

  uIm = h ·

∂uRe ∂a

  

(28)

 

The nodal displacements of the element are passed to the UEL in order to calculate the stresses and the complex strain energy density. This calculation is performed in an internal user material subroutine for each of the integration points of the element. The complex strain energy of the element is evaluated by using either of the approaches detailed in the former section, the strain energy density approach (Equation 20) or the force-displacement approach (Equation 24). At the end of the analysis, the sum of the imaginary component of the strain energy of each element is used to determine the ERR, Equation 16.

15

5

Numerical Results

5.1

Monotonic Loading Conditions

A compact tension (CT) specimen was simulated using the (1) Ramberg Osgood (deformation theory) and (2) bilinear isotropic hardening (incremental theory) material models. The ERR results obtained by (a) the American Society for Testing and Materials (ASTM) standards [38, 39] and (b) Abaqus J-integral were compared to (c) ZFEM results. The CT specimen had dimensions as specified in ASTM E1820-18. The specimen width, distance between the point of application of the load to the right edge, was 50 mm and the initial crack length was 5 mm as shown in Figure 4. The two dimensional mesh of the plane strain model generated in Abaqus consisted of eight noded quadratic quadrilateral elements. A total of 9210 nodes and 2393 elements were used to discretize the CT specimen. The crack tip was meshed using the swept √ meshing technique. In elastic plastic fracture mechanics, the 1/ r singularity, characteristic of the LEFM crack tip-strain field, does not exist. Thus, the nodes of the elements at the crack tip were collapsed to a triangle, but were left untied. The location of the midside nodes was unchanged. This nodal arrangement generated a 1/r strain singularity at the crack tip [1]. The q-vector was specified as a unit translation in the x-axis direction in order to evaluate the domain-based J-integral that is employed in Abaqus. The J-integral output was requested for 10 contours. The numerical estimate of the energy release rate using ZFEM converges at the same rate with respect to mesh density as the J-integral [31]. The specimen was loaded by applying an incremental vertical displacement to the center of the pins. The pins were modeled by a quarter circle of elastic elements. The ASTM standards use the load and load-line displacement curves data, which were obtained from the node specified as VL in Figure 4, as input to the ERR semi-analytical formulation. The ASTM equations are specified in the Appendix.

16

The ZFEM analysis had the same mesh discretization as Abaqus, but employed twodimensional 16-node ZFEM elements (8 real nodes, 8 imaginary nodes). In ZFEM, the virtual crack extension was implemented by elongating the highlighted region in Figure 4 using imaginary nodal coordinates and using a step size of ∆A = 10−10 multiplied by the size of the smallest length of the element at the crack tip in the imaginary x-axis direction. The ERR was calculated through the integration of the strain energy density, Equation 20. u

VL VL

a

u W

Figure 4: Finite element model of the CT specimen. The shadowed region highlights the nodes that were used to implement the virtual crack extension under the ZFEM approach.

5.1.1

Ramberg Osgood material model (deformation theory)

The Ramberg Osgood material model parameters are provided in Table 2. These constants correspond to a 2.25Cr-1Mo steel [40]. The ASTM standard estimates and the third J-integral contour results computed by Abaqus were used to evaluate the accuracy of the ZFEM approach. The J-integral solution is independent of the contours for this particular material model, aside of the first two J contours. In fact, the Abaqus manual does not recommend using the results of the first contour as the inaccuracies are larger for elements close to the crack tip. Table 2 was obtained by normalizing the ERR estimates by the third J-integral contour at a load-line displacement of 0.8 mm. A very 17

close agreement is observed between ZFEM and the J-integral, while the ASTM solution provided slightly higher estimates of the ERR. Material Parameter E ν σo α n

Value 175 GPa 0.3 140 MPa 2.2 5.0

Table 1: Ramberg Osgood Material Model Constants Load-Line Displacement (mm) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Normalized J Abaqus ASTM rd 3 contour 0.0611 0.0624 0.1627 0.1669 0.2842 0.2919 0.4073 0.4185 0.5513 0.566 0.6905 0.7098 0.8355 0.8590 1 1.0279

ZFEM 0.0611 0.1627 0.2842 0.4073 0.5513 0.6906 0.8357 0.9997

Table 2: Nonlinear ERR

5.1.2

Bilinear isotropic hardening (incremental theory)

The material parameters that were used for the bilinear isotropic hardening model are provided in Table 3. For numerical purposes, a constant plastic modulus was estimated from the experimental data provided by Brocks and Scheider [18] for a low alloy ferritic steel. The ERR results of the ASTM standard and Abaqus J-integral are compared to the ZFEM solution in Figure 5, in which the magnified curves illustrate a discrepancy between the results of different J-integral contours, where higher numbered contours converge to the ZFEM value. The approximation provided by the ASTM curve is slightly above the ZFEM curve and the ERR estimates provided by the J-integral. Notice that at load-line displacements larger than 0.5 mm, the lower numbered J-integral curves crossover the ASTM and ZFEM curves, providing higher ERR estimates. When a large contour value 18

is defined, such as J35 , a divergent result is obtained, as the J-integral contours touch the boundaries of the model. Figures 6-8 group the results according to the level of plasticity in the CT specimen, (a) contained or small scale yielding, (b) uncontained yielding, and (c) large scale yielding. The left figures show the yielded regions in the CT specimen at a specified load-line displacement (LLD). The right figures show the elements that are enclosed by the contours to calculate the J-integral and the estimates of the ERR values at five displacement magnitudes. Figure 6 represents the small yielding case, where the first 3 contours do not encompass the full extent of the plastic zone, resulting in an underestimate of the ERR for lower numbered contours. However, as the contours fully enclose the plastic zone, the J-integral value approaches the single value provided by ZFEM. Figure 7(a) shows uncontained yielding due to the presence of plastic strains at the back face and to the right of the loading pins. The path dependency of the J-integral is more pronounced for this case. Again, the J-integral values increase progressively as the contours increase until reaching a saturated value, which is in agreement with the ZFEM value. Finally, for the large scale yielding case (Figure 8), when a plastic band is formed between the region in the crack tip and the back face of the CT specimen, the lower numbered J-integral contour values overestimate the ERR. The J-integral converges to the single ZFEM value at higher numbered contours. Material Parameter E ν σy Ep

Value 213 GPa 0.3 715 MPa 3.45 GPa

Table 3: Bilinear Isotropic Material Model Constants

19

1.2 0.23 0.22

J3

0.21

J1

1

0.2

Normalized J

0.8 0.19

0.29

0.6

J35

Increasing J Contour

0.18

0.3

0.31

0.4

ASTM 0.2

0

ZFEM

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Load-Line Displacement (mm)

Figure 5: J-integral and ZFEM results for the incremental plasticity model.

20

0.09 LLD*=0.1834 mm 0.08 PE, Max. In−Plane Principal (Avg: 75%) +1.299e−01 +2.704e−02 +2.602e−02 +2.501e−02 +2.399e−02 +2.297e−02 +2.195e−02 +2.093e−02 +1.991e−02 +1.889e−02 +1.788e−02 +1.686e−02 +1.584e−02 +1.482e−02 +1.380e−02 +1.278e−02 +1.177e−02 +1.075e−02 +9.729e−03 +8.710e−03 +7.692e−03 +6.674e−03 +5.655e−03 +4.637e−03 +3.618e−03 +2.600e−03 +0.000e+00

LLD=0.1679 mm 0.07

Normalized J

0.06 LLD=0.1447 mm 0.05 LLD=0.1294 mm 0.04 LLD=0.1064 mm

0.03 Plastic Energy Release Rate

J-integral ZFEM

2

0.02

1.5

0.01 1

0

LLD* = 0.1834 mm ODB: bilinear_job.odb

Y

Z

X

Abaqus/Standard 3DEXPERIENCE R2016x

J1 J2 ...

J7 0.5

1

2

3

4

5

6

7

8

Contour

(a)

(b)

Tue May 15 12:06:12 Central Daylight Time 2018

9

10

0

-0.5

Step: Step−1 Increment 30: Step Time = 0.3000 Primary Var: PE, Max. In−Plane Principal

Figure 6: Small scale yielding conditions. (a) Plastic regions in CT specimen. (b) J-integral contour results (points) versus ZFEM single ERR value (dashed) at different load line displacement (LLD) magnitudes. The depicted plastic zone corresponds to the LLD*=0.1834 mm.

0.38 PE, Max. In−Plane Principal (Avg: 75%) +3.005e−01 +2.704e−02 +2.602e−02 +2.501e−02 +2.399e−02 +2.297e−02 +2.195e−02 +2.093e−02 +1.991e−02 +1.889e−02 +1.788e−02 +1.686e−02 +1.584e−02 +1.482e−02 +1.380e−02 +1.278e−02 +1.177e−02 +1.075e−02 +9.729e−03 +8.710e−03 +7.692e−03 +6.674e−03 +5.655e−03 +4.637e−03 +3.618e−03 +2.600e−03 +0.000e+00

LLD*=0.4140 mm

0.36

LLD=0.3991 mm

Normalized J

0.34 0.32

LLD=0.3764 mm

0.3

LLD=0.3611 mm

0.28 LLD=0.3381 mm

0.26 Plastic Energy Release Rate

J-integral ZFEM 2

0.24

1.5

0.22 1

LLD* = 0.4140 mm

0.2

1

2

3

4

5

J1 J2 ... 6

J7 7

8

0.5

9

Contour 0

ODB: bilinear_job.odb

Y

Z

X

Abaqus/Standard 3DEXPERIENCE R2016x

(a)

(b)

Tue May 15 12:06:12 Central Daylight Time 2018

-0.5

Step: Step−1 Increment 55: Step Time = 0.5500 Primary Var: PE, Max. In−Plane Principal

Figure 7: Uncontained yielding conditions. (a) Plastic regions in CT specimen. (b) J-integral contour results (points) versus ZFEM single ERR value (dashed) at different load line displacement (LLD) magnitudes. The depicted plastic zone corresponds to the LLD*=0.4140 mm.

21

10

1.1 Plastic Energy Release Rate

2

1.05

PE, Max. In−Plane Principal (Avg: 75%) +6.313e−01 +2.704e−02 +2.602e−02 +2.501e−02 +2.399e−02 +2.297e−02 +2.195e−02 +2.093e−02 +1.991e−02 +1.889e−02 +1.788e−02 +1.686e−02 +1.584e−02 +1.482e−02 +1.380e−02 +1.278e−02 +1.177e−02 +1.075e−02 +9.729e−03 +8.710e−03 +7.692e−03 +6.674e−03 +5.655e−03 +4.637e−03 +3.618e−03 +2.600e−03 −9.583e−08

1.5

1

J-integral ZFEM

1

J1 J2 ...

J7

Normalized J

0.5

0.95

LLD*=0.7550 mm 0

LLD=0.7376 mm -0.5

0.9 LLD=0.7115 mm 0.85

LLD=0.6942 mm

0.8

LLD* = 0.7550 mm ODB: bilinear_job.odb

Y

Z

X

Abaqus/Standard 3DEXPERIENCE R2016x

0.75

LLD=0.6680 mm

1

2

3

4

5

6

7

8

9

10

Contour

(a)

(b)

Tue May 15 12:06:12 Central Daylight Time 2018

Step: Step−1 Increment 95: Step Time = 0.9500 Primary Var: PE, Max. In−Plane Principal

Figure 8: Large scale yielding conditions. (a) Plastic regions in CT specimen. (b) J-integral contour results (points) versus ZFEM single ERR value (dashed) at different load line displacement (LLD) magnitudes. The depicted plastic zone corresponds to the LLD*=0.7550 mm.

5.2

ZFEM Force-displacement approach under monotonic loading conditions

A ZFEM-based approach for estimating the ERR without utilizing a strain energy density, W , formulation is also proposed. The ERR is obtained from the imaginary part of the strain energy, which is calculated through a numerical integration of the forcedisplacement relationship. This section examines the accuracy of this approach in obtaining the ERR for the CT fracture problem with a) bilinear isotropic hardening and b) Ramberg Osgood material models as compared to the strain energy density approach. The a) trapezoidal and b) Simpson’s rules were employed to show the adaptability of the approach to different integration schemes.

5.2.1

Trapezoidal Rule

The trapezoidal integration method, Equation 24, was selected for estimating the ERR of the bilinear isotropic hardening model. Figure 9 shows the ERR values provided

22

by both the strain energy density (W-ZFEM) and force-displacement (FD-ZFEM) approaches when breaking the load into 25, 100 and 1000 increments. The ERR results were normalized by the 10th J-integral contour value computed by Abaqus at full loading when using 100 increments, J10100 . A slight difference is observed in Figure 9 between both ZFEM approaches, W-ZFEM and FD-ZEM, for the curves corresponding to the same number of load increments. The ERR estimates obtained by using a larger number of load increments are higher than the ERR estimates obtained using a smaller number of load increments, regardless of the approach. However, this behavior is related to the incremental-iterative method used to solve nonlinear problems rather than the ZFEM-based approaches used to evaluate the ERR. The blow up magnification of the curves indicates that the ERR estimate obtained by W-ZFEM is slightly higher than that obtained by FD-ZFEM. The relative discrepancy of the FD-ZFEM approach using the trapezoidal integration rule was computed using the l2 -norm as discrepancy =

kJW-ZFEM − JFD-ZFEM k2 kJW-ZFEM k2

(29)

where JW-ZFEM and JFD-ZFEM are vectors with the ERR values obtained at a predefined number of steps, Nsteps , using W-ZFEM and FD-ZFEM, respectively. Figure 10 plots the discrepancy, (discrepancy/Nsteps ), as a function of Nsteps , where it can be observed that the discrepancy between both approaches reduces considerably as the number of loading steps increases.

23

Figure 9: Bilinear Isotropic Hardening Material Model: Nonlinear ERR evaluated using both ZFEM approaches (W-ZFEM and FD-ZFEM using the trapezoidal integration rule) and different number of load increments.

f⇤

Nsteps

u⇤

Figure 10: Normalized norm of the difference between the ZFEM approaches (W-ZFEM and FD-ZFEM) as a function of the number of load increments.

24

5.2.2

Simpson’s Rule

The Simpon’s integration method, Equation 25, was used to compute the ERR for the Ramberg Osgood model since this material model has a more pronounced nonlinearity than the bilinear model. Figure 11 shows the ERR curves for both approaches when solving the problem with 25, 100 and 1000 load increments. Again, the ERR results were normalized by J10100 . No visible difference is observed between the results of both approaches in the blow up magnification of the curves. Figure 12, plotting the discrepancy mean versus Nsteps , shows that the discrepancy decreases dramatically with an increase in the number of load steps.

Figure 11: Ramberg Osgood Material Model: Nonlinear ERR evaluated using both ZFEM approaches and different number of load increments.

25

f⇤

Nsteps

u⇤

Figure 12: Normalized norm of the difference between the ZFEM approaches (W-ZFEM and FD-ZFEM using the Simpson’s integration rule) as a function of the number of load increments.

5.3

Performance of the ZFEM-based ERR Estimating Method under loading/unloading conditions

5.3.1

Compact Tension Specimen

By definition, the J-integral is only valid under monotonic loading. This limitation is associated to the dependency of the evaluation of the J-integral on a strain energy density (SED) formulation. In contrast, as shown in Section 3.2.2.2, ZFEM does not require of an SED; hence, it provides accurate ERR results for loading and unloading. The J-integral evaluation during a full load/unload cycle is shown in Figure 13 for the CT specimen modeled with bilinear isotropic hardening plasticity and the same material properties as those provided in Table 3. Results for 28 J-contours are shown as a function of the pseudoloading time. The specimen is monotonically loaded up to a peak load line displacement of 0.802 mm at T =1 sec and unloaded between T =1 sec and T =2 sec. All results in this section were normalized by the maximum J-integral value when the structure is fully

26

loaded (T =1 sec). The J-integral contours exhibit erratic behavior during the unloading phase, do not converge to a saturated value, as in the loading case, and cannot be used to evaluate the ERR.

Figure 13: J-integral results for a CT specimen with bilinear isotropic hardening material properties for a full load/unload cycle.

The performance of FD-ZFEM under loading/unloading conditions was evaluated against the finite difference method. FD-ZFEM (here on referred as ZFEM) was considered suitable under cyclic loading conditions since no strain energy density formulation has to be known a priori. To obtain the finite differencing results, the global strain energy was obtained from two different runs, one with crack length a and the other with crack length a + ∆a. The nonlinear ERR, J, at each load step was approximated as (U (a + ∆a) − U (a))/∆a. Figure 14 compares ZFEM results against those provided by finite differencing and selected J-integral contours. ZFEM and finite differencing results are in excellent agreement throughout the loading and unloading phases. The J-integral values for the three different integration contours provide erroneous, non-converging approximations during unloading.

27

Figure 14: Methods comparison for a CT specimen

5.3.2

Edge Cracked Plate

The edge cracked plate geometry is shown in Figure 15(a). The dimensions considered for this analysis had a crack length to width ratio (a/W ) of 0.5, and length to width (L/W ) ratio of 2. The plate was discretized using 8-node plane strain elements for Abaqus and 16-noded plane strain elements for ZFEM. The problem was modeled using the same material properties as the CT specimen and the symmetric boundary conditions as shown in Figure 15(b). In ZFEM, the nodes in the red shadowed region surrounding the crack tip were perturbed in the self-similar crack growth direction.

28

u

2L

a

W Y

Z

X

u (a)

(b)

(c)

Figure 15: (a) Edged cracked plate and (b) its corresponding finite element mesh with shadowed element indicating the perturbed elements. (c) Plastic zone regions at full loading illustrated by the non-gray colors.

Figure 15(c) shows that yielding extends far away from the crack tip. The three methods are in excellent agreement during loading as observed in Figure 16. However, during unloading, ZFEM and finite differencing results remain in excellent agreement, whereas the J-integral behavior deviates from the ZFEM and finite differencing results. While there is no analytical solution, the agreement of numerical differentiation results using complex and real variables gives confidence in the accuracy of the ZFEM results.

29

Figure 16: Methods comparison for an edge cracked plate

6

Discussion

The numerical examples in this paper showed that ZFEM is a robust method for evaluating the nonlinear ERR at all levels of plasticity and is superior to the J-integral for unloading. ZFEM provides a single, accurate ERR estimate for small scale, contained, and large scale yielding. On the other hand, the J-integral results are contour dependent for incremental plasticity. A saturated J-integral value is reached at different contours for varying load levels. In addition, a ZFEM-based approach has been presented that does not depend on the existence of a strain energy density formulation and, therefore, provides accurate ERR estimates under loading/unloading conditions. The ERR is computed from the complexvalued force displacement relationship that is obtained throughout the loading steps. In contrast, the J-integral is limited to monotonic loading conditions as a strain energy density cannot be uniquely defined for history-dependent elastic-plastic materials. The proposed method offers significant implementation advantages as compared to the J-integral and finite differencing approaches. The ZFEM implementation is based 30

on the complex representation of the finite element code. The virtual crack extension required for fracture applications is introduced through a perturbation along imaginary coordinate axes. ZFEM computes the ERR from the imaginary part of the strain energy without the need of performing any additional calculations or domain-based integrations. Conversely, the J-integral requires appending corrective terms to recover its path independent property for different physics. It is recognized that complex finite element method increases the computational time of the simulations as it duplicates the degrees-of-freedom of the analysis. Run time comparisons performed in former studies indicate that ZFEM runs are between 2.5 and 4 times longer than equivalent real-variable runs [35, 30]. Ongoing ZFEM efforts include reducing computational run times in order to increase the efficiency of the method.

7

Conclusions

ZFEM has emerged as a new alternative VCE method to compute the ERR in Elastic Plastic Fracture Mechanics. The method overcomes well-known limitations of the J-integral, specifically, contour path dependency and its restriction to monotonic loading. Two ZFEM approaches were presented to obtain the nonlinear ERR: the first one depended on a strain energy density formulation and the second one used a forcedisplacement relationship. The first approach obtains the ERR by integrating the energy density over the volume, with the imaginary part containing the ERR. This ZFEM approach was found accurate and robust when solving problems under monotonic loading and employing material models based on the deformation and incremental theory of plasticity. A ZFEM single value was found to be of the same accuracy as the ASTM standards and the saturated value of the J-integral. The force-displacement approach for evaluating the ERR was presented that was found suitable for problems under loading/unloading conditions. This approach consists of integrating the complex-valued force-displacement relationship at the element level.

31

These force and displacement vectors are those available within any FEA model. The advantage of this approach is that there is no need of invoking a strain energy density formulation. Under monotonic loading, the results were in excellent agreement with J-integral results and finite differencing. During the unloading phase, the force-displacement ZFEM results were in excellent agreement with those obtained by finite differencing, while the J-integral were inaccurate. Hence, ZFEM can become a valuable tool for modeling many challenging EPFM problems such as cyclic fatigue and ductile crack propagation.

8

Acknowledgments

This project was funded in part by the University of Texas at San Antonio: Office of the Vice President for Research, the Nuclear Regulatory Commission (Grant No. NRC-HQ60-17-G-0024), and Department of Defense (Grant No. W911NF-15-1-0456).

32

9

Appendix

9.1

ASTM E1820-18 Equations

The calculations of J for a CT specimen with a single edge notched specimen J = J el + J pl =

(30)

K 2 (1 − ν 2 ) ηApl + E t(Ψ − a)

where K =

P √ f (a/Ψ). B Ψ

For a compact tension specimen,

"    2  3  4 # 2 + Ψa a a a a a f( ) = 0.886 + 4.64 − 13.32 + 14.72 − 5.60
3/2 Ψ Ψ Ψ Ψ Ψ 1 − Ψa (31) and η = 2 + 0.522(b/Ψ), b = Ψ − a is the uncracked ligament length, Apl is the plastic work found from the load vs line displacement curve and Ψ is the horizontal distance between the point of application of load to the right edge of the CT specimen, a is the distance between the load line and the crack tip, and t is the thickness of the specimen.

33

References [1] T.L Anderson. Fracture Mechanics Concepts and Applications. CRC Press, Boca Raton, FL, 1995. [2] T. K. Hellen and W. S. Blackburn. Non-linear fracture mechanics and finite elements. Engineering Computations, 4(1):2–14, 1987. [3] D.M. Stump and E. Zywicz. J-integral computations in the incremental and deformation plasticity analysis of small-scale yielding. Engineering Fracture Mechanics, 45(1):61 – 77, 1993. [4] J.Rice. A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics, 35:379–386, 1968. [5] C. F. Shih, B. Moran, and T. Nakamura.

Energy release rate along a three-

dimensional crack front in a thermally stressed body. International Journal of Fracture, 30(2):79–102, 1986. [6] R. M. McMeeking. Path dependence of the j-integral and the role of j as a parameter characterizing the near-tip field. pages 28–41, 1977. [7] R. Riddle, R. Streit, and I. Finnie. The failure of aluminum compact shear specimens under mixed-mode loading. Fracture Mechanics: Eighteenth Symposium, ASTM STP 945, pages 118 – 133, 1988. [8] A. Haddi and D. Weichert. Elastic-plastic j-integral in inhomogeneous materials. Computational Materials Science, 8(3):251–260, 1997. [9] J. W. Eischen. Fracture of nonhomogeneous materials. International Journal of Fracture, 34(1):3–22, 1987.

34

[10] J. Kim and G. Paulino. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. International Journal for Numerical Methods in Engineering, 53(8):1903–1935, 2002. [11] J. Kim and G. Paulino. The interaction integral for fracture of orthotropic functionally graded materials: evaluation of stress intensity factors. International Journal of Solids and Structures, 40(15):3967 – 4001, 2003. [12] R. Schapery. Correspondence principle and generalized J integral for large deformation and fracture analysis of viscoelastic media. Int J Fract, 25:194–223, 1984. [13] F.W. Brust, J.J. McGowan, and S.N. Atluri. A combined numerical/experimental study of ductile crack growth after a large unloading, using T, J and CTOA criteria. Engineering Fracture Mechanics, 23(3):537 – 550, 1986. [14] A. Sedmake, M. Berkovi´c, and J. Jari´c. J integral for thin shells. Mechanical Engineering Publications, pages 44 – 53, 1991. (Accepted). [15] Y. Lei, N.P. O’Dowd, and G.A. Webster. Fracture mechanics analysis of a crack in a residual stress field. International Journal of Fracture, 106(3):195–216, 2000. [16] J. Newman, S. Smith, B. Seshadri, M. James, R. Brazill, R. Schultz, J. Donald, and A. Blair. Characterization of Residual Stress Effects on Fatigue Crack Growth of a Friction Stir Welded Aluminum Alloy, 2015. [17] M. Lorentzon and K. Eriksson. A path independent integral for the crack extension force of the circular arc crack. Engineering Fracture Mechanics, 66(5):423 – 439, 2000. [18] I. Scheider, W. Brocks, and A. Cornec. Comprehensive Structural Integrity aspects of nonlinear fracture mechanics, In: de Borst, R., Mang H.A (eds) Comprehensive structural integrity numerical and computational methods, vol 3. Elsevier, New York, pp 127-209. 35

[19] T. K. Hellen. Virtual crack extension methods for non-linear materials. International Journal for Numerical Methods in Engineering, 28(4):929–942, 1989. [20] J. W. Hutchinson. Singular behaviour at the end of a tensile crack in a hardening material. Journal of the Mechanics and Physics of Solids, 16(1):13–31, 1968. [21] J. R. Rice and G. F. Rosengren. Plane strain deformation near a crack tip in a powerlaw hardening material. Journal of the Mechanics and Physics of Solids, 16(1):1–12, 1968. [22] C. E. Turner. A re-assessment of ductile tearing resistance. i. the geometry dependence of j − r curves in fully plastic bending. (retroactive coverage). pages 933–949, 1990. [23] O. Kolednik, R. Sch¨ongrundner, and F. D. Fischer. A new view on j-integrals in elastic–plastic materials. International Journal of Fracture, 187(1):77–107, 2014. [24] S. Goutianos and B. F. Sørensen. The application of j integral to measure cohesive laws under large-scale yielding. Engineering Fracture Mechanics, 155:145 – 165, 2016. [25] D. M. Parks. A stiffness derivative finite element technique for determination of crack tip stress intensity factors. International Journal of Fracture, 10(4):487–502, 1974. [26] D.M. Parks. The virtual crack extension method for nonlinear material behavior. Computer Methods in Applied Mechanics and Engineering, 12(3):353 – 364, 1977. [27] H. G. DeLorenzi. 3-d elastic-plastic fracture mechanics with adina. Computers and Structures, 13(5):613–621, 1981. [28] H. G. DeLorenzi. On the energy release rate and the j-integral for 3-d crack configurations. International Journal of Fracture, 19(3):183–193, 1982.

36

[29] W. Squire and G. Trapp. Using complex variables to estimate derivatives of real functions. SIAM Review, 40(1):110–112, 1998. [30] H. Millwater, D. Wagner, A. Baines, and A. Montoya. A virtual crack extension method to compute energy release rates using a complex variable finite element method. Engineering Fracture Mechanics, 162:95 – 111, 2016. [31] D. Ramirez Tamayo, A. Montoya, and H. Millwater. A virtual crack extension method for thermoelastic fracture using a complex-variable finite element method. Engineering Fracture Mechanics, 2017. [32] D. Wagner and H. Millwater. 2D weight function development using a complex Taylor series expansion method. Engineering Fracture Mechanics, 86:23–37, 2012. [33] H. Millwater, D. Wagner, A. Baines, and K. Lovelady. Improved WCTSE method for the generation of 2d weight functions through implementation into a commercial finite element code. Engineering Fracture Mechanics, 109:302 – 309, 2013. [34] Z. Jing and X.R. Wu. Wide-range weight functions and stress intensity factors for arbitrarily shaped crack geometries using complex taylor series expansion method. Engineering Fracture Mechanics, 138:215 – 232, 2015. [35] A.Montoya, R. Fielder, A. Gomez-Farias, and H. Millwater. Finite-element sensitivity for plasticity using complex variable methods. Journal of Engineering Mechanics, 141(2):04014118, 2015. [36] J.A. Begley and J.D. Landes. The Effect of Specimen Geometry on J Ic . ASTM STP 514, American Society for Testing and Materials, Philadelphia, PA, 1972. [37] J.A. Begley and J.D. Landes. The J-Integral as a Fracture Criterion. ASTM STP 514, American Society for Testing and Materials, Philadelphia, PA, 1972. [38] ASTM Standards, Standard Test Method for Measurement of Fracture Toughness, ASTM E1820-18 American Society for Testing and Materials, Philadelphia, 2018. 37

[39] ASTM Standards, Standard Test Method for Plane-Strain Fracture Toughness of Metallic Materials, ASTM E399-90, American Society for Testing and Materials, Philadelphia, 1997. [40] A. Saxena. Nonlinear fracture mechanics for engineers. CRC Press, Boca Raton, Fla, 1998.

38

Highlights • • • • •

This study presents a new approach to calculate the nonlinear energy release rate. The method evaluates the derivative of the potential energy using a complex-variable FEM. The method overcomes the path dependency and monotonic loading limitation of the J-integral. The method remains robust all levels of plasticity, i.e. small, contained, and large scale yielding. The method is superior to the J-integral when unloading from a yield state.