Virtual Crack Closure Technique: History, Approach, And Applications

112Krueger : Virtual crack closure technique: History, approach, and applicationstwo-dimensional, four-noded elements, as shown in Fig. 4,the work E required to close the crack along one elementside therefore can be calculated as E 21 关 X i u ᐉ Z i w ᐉ 兴 ,(2)where X i and Z i are the shear and opening forces at nodalpoint i and u ᐉ and w ᐉ are the shear and opening displacements at node ᐉ as shown in Fig. 4. Thus, forces and displacements required to calculate the energy E to close thecrack may be obtained from one single finite element analysis. The details of calculating the energy release rate G E/ A, where A is the crack surface created, and theseparation into the individual mode components will be discussed in the following section.3 EQUATIONS FOR USING THE VIRTUALCRACK CLOSURE TECHNIQUEIn the following, equations are presented to calculate mixedmode strain energy release rates using two-dimensional finiteelement models such as plane stress or plane strain. Differentapproaches are also discussed for the cases where the crackor delamination is modeled with plate/shell elements or withthree-dimensional solids.Appl Mech Rev vol 57, no 2, March 2004modeled using single nodes, or two nodes with identical coordinates coupled through multipoint constraints if a crackpropagation analysis is desired. This is discussed in detail inthe Appendix, which explains specific modeling issues.For a crack propagation analysis, it is important to advance the crack in a kinematically compatible way. Nodewise opening/closing, where node after node is sequentiallyreleased along the crack, is possible for the four-noded element as shown in Fig. 6共a兲. It is identical to elementwiseopening in this case as the crack is opened over the entirelength of the element. Nodewise opening/closing, however,results in kinematically incompatible interpenetration for theeight-noded elements with quadratic shape functions asshown in Fig. 6共b兲, which caused initial problems wheneight-noded elements were used in connection with the virtual crack closure technique. Elementwise opening—whereedge and midside nodes are released—provides a kinematically compatible condition and yields reliable results, whichwas demonstrated in Refs. 关5兴,关54兴,关55兴 and later generalizedexpressions to achieve this were derived by Raju 关21兴.The mode I and mode II components of the strain energyrelease rate, G I and G II , are calculated for four-noded elements as shown in Fig. 7共a兲:G I 1Z 共 w w ᐉ * 兲 ,2 a i ᐉ(3)1X 共 u u ᐉ * 兲 ,2 a i ᐉ(4)G II 3.1Formulas for two-dimensional analysisIn a two-dimensional finite element plane stress, or planestrain model, the crack of length a is represented as a onedimensional discontinuity by a line of nodes as shown in Fig.5. Nodes at the top surface and the bottom surface of thediscontinuity have identical coordinates, however, and arenot connected with each other as shown in Fig. 5共a兲. Thislets the elements connected to the top surface of the crackdeform independently from those connected to the bottomsurface and allows the crack to open as shown in Fig. 5共b兲.The crack tip and the undamaged section, or the sectionwhere the crack is closed and the structure is still intact, isFig. 4Modified crack closure method 共one-step VCCT兲Fig. 5 Crack modeled as one-dimensional discontinuity. a兲 Initially modeled, undeformed finite element mesh and b兲 deformedfinite element mesh.

Appl Mech Rev vol 57, no 2, March 2004Krueger : Virtual crack closure technique: History, approach, and applicationswhere a is the length of the elements at the crack front andX i and Z i are the forces at the crack tip 共nodal point i兲. Therelative displacements behind the crack tip are calculatedfrom the nodal displacements at the upper crack face u ᐉ andw ᐉ 共nodal point ᐉ兲 and the nodal displacements u ᐉ * and w ᐉ *at the lower crack face 共nodal point ᐉ*兲, respectively. Thecrack surface A created is calculated as A a 1,where it is assumed that the two-dimensional model is ofunit thickness 1. While the original paper by Rybicki andKanninen is based on heurisitic arguments 关19兴, Raju provedthe validity of the equation 关21兴. He also showed that theequations are applicable if triangular elements, obtained bycollapsing the rectangular elements, are used at the crack tip.The mode I and mode II components of the strain energyrelease rate, G I , and G II , are calculated for eight-noded elements as shown in Fig. 7共b兲:G I 1131关 Z 共 w w ᐉ * 兲 Z j 共 w m w m * 兲兴 ,2 a i ᐉ(5)1关 X 共 u u ᐉ * 兲 X j 共 u m u m * 兲兴 ,2 a i ᐉ(6)G II where a is the length of the elements at the crack front asabove. In addition to the forces X i and Z i at the crack tip共nodal point i兲 the forces X j and Z j at the midside node infront of the crack 共nodal point j兲 are required. The relativesliding and opening behind the crack tip are calculated atnodal points ᐉ and ᐉ* from displacements at the upper crackface u ᐉ and w ᐉ and the displacements u ᐉ * and w ᐉ * at thelower crack face. In addition to the relative displacements atnodal points ᐉ and ᐉ* the relative displacements at nodalpoints m and m * are required, which are calculated fromFig. 6 Kinematic compatiblecrack opening/closure. a兲 Nodewise crack opening for four-nodedelement and b兲 crack opening foreight-noded element.

114Krueger : Virtual crack closure technique: History, approach, and applicationsdisplacements at the upper crack face u m and w m and thedisplacements u m * and w m * at the lower crack face 关21兴. Thecrack surface A created is calculated as A a 1,where it is assumed that the two-dimensional model is ofunit thickness 1. The equations are also applicable if triangular parabolic elements, obtained by collapsing the parabolic rectangular elements, are used at the crack tip 关21兴. Thetotal energy release rate G T is calculated from the individualmode components asG T G I G II G III ,(7)where G III 0 for the two-dimensional case discussed.The VCCT proposed by Rybicki and Kanninen did notmake any assumptions of the form of the stresses and displacements. Therefore, singularity elements are not requiredAppl Mech Rev vol 57, no 2, March 2004at the crack tip. However, special two-dimensional crack tipelements with quarter-point nodes as shown in Fig. 8 havebeen proposed in the literature 关21,56 –58兴. Based on thelocation of the nodal points at 0.0, 0.25, and 1.0, thesequarter-point elements accurately simulate the 1/冑r singularity of the stress field at the crack tip. Triangular quarter-pointelements are obtained by collapsing one side of the rectangular elements, as shown in Fig. 8共b兲. The mode I and modeII components of the strain energy release rate, G I , and G IIare calculated for eight-noded singularity elements using thesimplified equations given in Ref. 关21兴:G I 1关 Z 兵 t 共 w w ᐉ * 兲 t 12共 w m w m * 兲 其2 a i 11 ᐉ Z j 兵 t 21共 w ᐉ w ᐉ * 兲 t 22共 w m w m * 兲 其 兴 ,(8)Fig. 7 Virtual crack closure technique for 2Dsolid elements. a兲 Virtual crack closure techniquefor four-noded element 共lower surface forces areomitted for clarity兲 and b兲 virtual crack closuretechnique for eight-noded element 共lower surfaceforces are omitted for clarity兲.

Appl Mech Rev vol 57, no 2, March 2004G II Krueger : Virtual crack closure technique: History, approach, and applications1关 X 兵 t 共 u u ᐉ * 兲 t 12共 u m u m * 兲 其2 a i 11 ᐉ X j 兵 t 21共 u ᐉ u ᐉ * 兲 t 22共 u m u m * 兲 其 兴 ,(9)wheret 11 6 3 ,2t 12 6 20,1t 21 ,2t 22 1.(10)In contrast to regular parabolic elements, Eqs. 共8兲 and 共9兲 forthe quarter-point elements have cross terms involving the115corner and quarter-point forces and the relative displacements at the corner and quarter-point nodes. Equations 共8兲and 共9兲 are also valid if triangular quarter-point elements,obtained by collapsing one side of the rectangular elements,are used at the crack tip as shown in Fig. 8共b兲. Note that fortriangular elements the nodal forces have to be calculatedfrom elements A, B, C, and D around the crack tip. Specialrectangular and collapsed singularity elements with cubicshape functions are also discussed in Ref. 关21兴. A specialsix-noded rectangular element with quarter-point nodes isdescribed in Ref. 关59兴. Due to the fact that these elements arenot readily available in most of the commonly used finiteFig. 8 Singularity elements with quarter-pointnodes at crack tip. a兲 Quadtrilateral elements withquarter-point nodes and b兲 collapsed quarterpoint element.

116Krueger : Virtual crack closure technique: History, approach, and applicationselement codes, equations are not provided. For additionalinformation about singularity elements the interested readeris referred to Refs. 关58兴,关60– 64兴.3.2 Formulas for three-dimensional solids and plateÕshell elementsIn a finite element model made of three-dimensional solidelements 关Fig. 9共a兲兴 or plate or shell type elements 关Fig.9共b兲兴 the delamination of length a is represented as a twodimensional discontinuity by two surfaces. The additionaldimension allows us to calculate the distribution of the energy release rates along the delamination front and makes itpossible to obtain G III , which is identical to zero for twodimensional models. Nodes at the top surface and the bottomAppl Mech Rev vol 57, no 2, March 2004surface have identical coordinates and are not connected witheach other as explained in the preceding section. The delamination front is represented by either a row of single nodes ortwo rows of nodes with identical coordinates, coupledthrough multipoint constraints. The undamaged sectionwhere the delamination is closed and the structure is intact ismodeled using single nodes or two nodes with identical coordinates coupled through multipoint constraints if a delamination propagation analysis is desired. This is discussed indetail in the Appendix, which explains specific modeling issues.3.2.1 Formulas for three-dimensional solidsFor convenience, only a section of the delaminated area thatis modeled with eight-noded three-dimensional solid ele-Fig. 9 Delaminations modeledas two-dimensional discontinuity.a兲 Delamination modeled with bilinear 3D solid elements and b兲delamination modeled with bilinear plate/shell type elements.

Appl Mech Rev vol 57, no 2, March 2004Krueger : Virtual crack closure technique: History, approach, and applicationsments is illustrated in Fig. 10. The mode I, mode II, andmode III components of the strain energy release rate, G I ,G II , and G III , are calculated asG I 1Z 共 w w Lᐉ * 兲 ,2 A Li Lᐉ(11)1X 共 u u Lᐉ * 兲 ,2 A Li Lᐉ(12)1Y 共 v v Lᐉ * 兲 ,2 A Li Lᐉ(13)G II G III with A ab as shown in Fig. 10 关65兴. Here A is thearea virtually closed, a is the length of the elements at thedelamination front, and b is the width of the elements. For117better identification in this and the following figures, columns are identified by capital letters and rows by small letters as illustrated in the top view of the upper surface shownin Fig. 10共b兲. Hence, X Li , Y Li , and Z Li denote the forces atthe delamination front in column L, row i. The correspondingdisplacements behind the delamination at the top face noderow ᐉ are denoted u Lᐉ , v Lᐉ , and w Lᐉ and at the lower facenode row ᐉ* are denoted u Lᐉ * , v Lᐉ * , and w Lᐉ * as shown inFig. 10. All forces and displacements are obtained from thefinite element analysis with respect to the global system. Alocal crack tip coordinate system (x ,y ,z ) that defines thenormal and tangential coordinate directions at the delamination front in the deformed configuration has been added tothe illustration. Its use with respect to geometrically nonlinear analyses will be discussed later.For twenty-noded solid elements, the equations to calculate the strain energy release rate components at the elementcorner nodes 共location Li) as shown in Fig. 11 areFig. 10 Virtual crack closure technique for fournoded plate/shell and eight-noded solid elements.a兲 3D view 共lower surface forces are omitted forclarity兲 and b兲 top view of upper surface 共lowersurface terms are omitted for clarity兲.

118Krueger : Virtual crack closure technique: History, approach, and applicationsG I 1关 1 Z 共 w w Kᐉ * 兲 Z Li 共 w Lᐉ w Lᐉ * 兲2 A L 2 Ki Kᐉ Z L j 共 w Lm w Lm * 兲 21 Z M i 共 w M ᐉ w M ᐉ * 兲兴 ,G II (14)1关 1 X 共 u u Kᐉ * 兲 X Li 共 u Lᐉ u Lᐉ * 兲2 A L 2 Ki Kᐉ X L j 共 u Lm u Lm * 兲 21 X M i 共 u M ᐉ u M ᐉ * 兲兴 ,G III (15)1关 1 Y 共 v v Kᐉ * 兲 Y Li 共 v Lᐉ v Lᐉ * 兲2 A L 2 Ki Kᐉ Y L j 共 v Lm v Lm * 兲 21 Y M i 共 v M ᐉ v M ᐉ * 兲兴 ,(16)where A L ab as shown in Fig. 11 关66兴. Here X Ki , Y Ki ,Appl Mech Rev vol 57, no 2, March 2004and Z Ki denote the forces at the delamination front in columnK, row i. The relative displacements at the correspondingcolumn K are calculated from the displacements behind thedelamination at the lower face node row ᐉ* as u Kᐉ * , v Kᐉ * ,and w Kᐉ * and at the top face node row ᐉ, as u Kᐉ , v Kᐉ andw Kᐉ 关Fig. 11共b兲兴. Similar definitions are applicable in columnM for the forces at node row i and displacements at node rowᐉ and in column L for the forces at node row i and j anddisplacements at node row ᐉ and m, respectively. Only onehalf of the forces at locations Ki and M i contribute to theenergy required to virtually close the area A L . Half of theforces at location Ki contribute to the closure of the adjacentarea A J and half of the forces at location M i contribute tothe closure of the adjacent area A N .The equations to calculate the strain energy release ratecomponents at the midside node 共location M i) as shown inFig. 12 are as follows 关66,67兴:Fig. 11 Virtual crack closure technique for corner nodes in eight-noded plate/shell and twentynoded solid-elements. a兲 3D view 共lower surfaceforces are omitted for clarity兲 and b兲 top view ofupper surface 共lower surface terms are omittedfor clarity兲.

Appl Mech Rev vol 57, no 2, March 2004G I Krueger : Virtual crack closure technique: History, approach, and applications冋111Z 共 w w Lᐉ * 兲 Z L j 共 w Lm w Lm * 兲2 A M 2 Li Lᐉ2册(17)冋111X Li 共 u Lᐉ u Lᐉ * 兲 X L j 共 u Lm u Lm * 兲G II 2 A M 221 X M i 共 u M ᐉ u M ᐉ * 兲 X Ni 共 u Nᐉ u Nᐉ * 兲2册1 X N j 共 u Nm u Nm * 兲 ,2冋111Y 共 v v Lᐉ * 兲 Y L j 共 v Lm v Lm * 兲2 A M 2 Li Lᐉ21 Y M i 共 v M ᐉ v M ᐉ * 兲 Y Ni 共 v Nᐉ v Nᐉ * 兲21 Z M i 共 w M ᐉ w M ᐉ * 兲 Z Ni 共 w Nᐉ w Nᐉ * 兲21 Z N j 共 w Nm w Nm * 兲 ,2G III 119(18)册1 Y N j 共 v Nm v Nm * 兲 ,2(19)where only one half of the forces at locations Li, L j and Ni,N j contribute to the energy required to virtually close thearea A M . Half of the forces at locations Li and L j contribute to the closure of the adjacent area A K and half of theforces at locations Ni and N j contribute to the closure of theadjacent area A 0 .Instead of computing the strain energy release rate components at the corner or midside nodes as described above,G I , G II , and G III may be calculated for an entire element,Fig. 12 Virtual crack closure technique for midside nodes in eight-noded plate/shell and twentynoded solid elements. a兲 3D view 共lower surfaceforces are omitted for clarity兲 and b兲 top view ofupper surface 共lower surface terms are omittedfor clarity兲.

120Krueger : Virtual crack closure technique: History, approach, and applicationswhich may be advantageous in cases where the elements arenot of square or rectangular shape. For example, for the computation of the strain energy release rate components along acircular or elliptical front where elements are trapezoidal theuser may find this approach more suitable. The equations tocalculate the strain energy release rate components for oneelement as shown in Fig. 13 are as follows 关66 – 68兴:G I 1关 Z 共 w w Lᐉ * 兲 Z L j 共 w Lm w Lm * 兲2 A M Li Lᐉ Z M i 共 w M ᐉ w M ᐉ * 兲 Z Ni 共 w Nᐉ w Nᐉ * 兲 Z N j 共 w Nm w Nm * 兲兴 ,G II (20)1关 X 共 u u Lᐉ * 兲 X L j 共 u Lm u Lm * 兲2 A M Li Lᐉ X M i 共 u M ᐉ u M ᐉ * 兲 X Ni 共 u Nᐉ u Nᐉ * 兲 X N j 共 u Nm u Nm * 兲兴 ,(21)G III Appl Mech Rev vol 57, no 2, March 20041关 Y 共 v v Lᐉ * 兲 Y L j 共 v Lm v Lm * 兲2 A M Li Lᐉ Y M i 共 v M ᐉ v M ᐉ * 兲 Y Ni 共 v Nᐉ v Nᐉ * 兲 Y N j 共 v Nm v Nm * 兲兴 ,(22)where the forces at locations Li, L j and Ni, N j are calculated only from elements A and B, which are shaded in Fig.13共b兲. This is unlike the previous equations where four elements contributed to the forces at locations Li, L j and Ni,N j. The force at location M i is also calculated from elements A and B, which is identical to the procedure above.A three-dimensional twenty-noded singular brick elementwith quarter points is shown in Fig. 14. As mentioned abovethe desired 1/冑r singularity of the stress field at the crack tipis achieved by moving the midside node to the quarter position. A prism-shaped singular element is obtained by collaps-Fig. 13 Virtual crack closure technique 共element method兲 for eight-noded plate/shell andtwenty-noded solid elements. a兲 3D view 共lowersurface forces are omitted for clarity兲 and b兲 topview of upper surface 共lower surface terms areomitted for clarity兲.

Appl Mech Rev vol 57, no 2, March 2004Krueger : Virtual crack closure technique: History, approach, and applicationsing a face of the element as shown in Fig. 15. Each set ofcollapsed nodes at the front must either be defined as asingle node or the degrees of freedom must be connectedthrough multipoint constraints as if each set was a singleG I node. The mode I, mode II, and mode III components of thestrain energy release rate, G I , G II , and G III , are calculatedfor a brick or prism singularity element using the simplifiedequations given in Ref. 关68兴:冋 再1t 12Z Li t 11共 w Lm w Lm * 兲 t 12共 w M ᐉ w M ᐉ * 兲 2t 12共 w Lᐉ w Lᐉ * 兲 共 w Nᐉ w Nᐉ * 兲2 A M2再再 Z Ni t 12共 w M ᐉ w M ᐉ * 兲 Z N jt 12共 w Lᐉ w Lᐉ * 兲 2t 12共 w Nᐉ w Nᐉ * 兲 t 11共 w Nm w Nm * 兲2冎 再冎 再121冎冎1t 11共 w Nᐉ w Nᐉ * 兲 共 w Nm w Nm * 兲 Z M i共 w Lm w Lm * 兲 t 21共 w M ᐉ w M ᐉ * 兲 t 22共 w Nᐉ w Nᐉ * 兲22 t 22共 w Lᐉ w Lᐉ * 兲 t 111共 w Nm w Nm * 兲 Z L j 共 w Lm w Lm * 兲 共 w Lᐉ w Lᐉ * 兲22冎册,(23)Fig. 14 Virtual crack closure technique 共element method兲 for twenty-noded quarter point elements. a兲 3D view 共lower surface forces areomitted for clarity兲 and b兲 top view of upper surface 共lower surface terms are omitted for clarity兲.

122Krueger : Virtual crack closure technique: History, approach, and applicationsG II Appl Mech Rev vol 57, no 2, March 2004冋 再1X t 共 u u Lm * 兲 t 12共 u M ᐉ u M ᐉ * 兲2 A M Li 11 Lm 2t 12共 u Lᐉ u Lᐉ * 兲 再再t 12共 u u Nᐉ * 兲2 Nᐉ X Ni t 12共 u M ᐉ u M ᐉ

2.1 Crack closure method using two analysis steps Even though the virtual crack closure technique is the focus of this paper and is generally mentioned in the literature, it appears appropriate to include a related method: the crack closure method or two-step crack closure technique. The ter-minology in the literature is often inexact and this .