The bi-material interface arc crack has been the focus of interest in the com- posite community, where it is usually referred to as the ber-matrix interface crack. In this work, we investigate the convergence properties of the Virtual Crack Closure Technique (VCCT) when applied to the evaluation of the Mode I, Mode II and total Energy Release Rate of the ber-matrix interface crack in the context of the Finite Element Method (FEM). We rst propose a synthetic vectorial formulation of the VCCT. Thanks to this formulation, we study the convergence properties of the method, both analytically and numerically. It is found that Mode I and Mode II Energy Release Rate (ERR) possess a loga- rithmic dependency with respect to the size of the elements in the crack tip neighborhood, while the total ERR is independent of element size. Keywords: Fiber/matrix interface crack, Bi-material interface arc crack, Linear Elastic Fracture Mechanics (LEFM), Virtual Crack Closure Technique (VCCT), Mode separation, Convergence 1. Introduction Bi-material interfaces represent the basic load transfer mechanism at the heart of Fiber Reinforced Polymer Composite (FRPC) materials. They are present at the macroscale, in the form of adhesive joints; at the mesoscale, as interfaces between layers with di erent orientations; at the microscale, as ber- 5 Preprint submitted to Finite Elements in Analysis and Design July 2, 2019 DRAFTmatrix interfaces. Bi-material interfaces have for long attracted the attention of researchers in Fracture Mechanics [1, 2], due to their hidden complexity. The problem was rst addressed in the 1950's by Williams [3], who derived through a linear elastic asymptotic analysis the stress distribution around an open crack (i.e. with crack faces nowhere in contact for any size of the crack) 10 between two in nite half-planes of dissimilar materials. He found the existence of a strong oscillatory behavior in the stress singularity at the crack tip of the form r1 2sin ("logr) with"=1 2log1 1 +  ; (1) in both Mode I and Mode II. In Eq. 1, is one of the two parameters introduced by Dundurs [4] to characterize bi-material interfaces: 15 =2(11)1(21) 2(1+ 1) +1(2+ 1)(2) where= 34in plane strain and =34 1+in plane stress, is the shear modulus, Poisson's coecient, and indexes 1 ;2 refer to the two bulk materials joined at the interface. De ning aas the length of the crack, it was found that the size of the oscillatory region is in the order of 106a[5]. Given the oscillatory behaviour of the crack tip singularity of Eq. 1, the de nition of 20 Stress Intensity Factor (SIF) lim r!0p 2rdiverges and ceases to be valid [1]. It implies that the Mode mixity problem at the crack tip is ill-posed. It was furthermore observed, always in the context of Linear Elastic Fracture Mechanics (LEFM), that an interpenetration zone exists close to the crack tip [6, 7] with a length in the order of 104a[6]. Following conclusions rstly proposed 25 in [7], the presence of a contact zone in the crack tip neighborhood, of a length to be determined from the solution of the elastic problem, was introduced in [8] and shown to provide a physically consistent solution to the straight bi-material interface crack problem. The curved bi-material interface crack, more often refered to as the ber-matrix 30 interface crack (or debond) due to its relevance in FRPCs, was rst treated by 2 DRAFTEngland [9] and by Perlman and Sih [10], who provided the analytical solution of stress and displacement elds for a circular inclusion with respectively a single debond and an arbitrary number of debonds. Building on their work, Toya [11] particularized the solution and provided the expression of the Energy Release 35 Rate (ERR) at the crack tip. The same problems exposed previously for the open straight bi-material crack were shown to exist also for the open ber-matrix interface crack: the presence of strong oscillations in the crack tip singularity and onset of crack face interpenetration at a critical aw size.1 In order to treat cases more complex than the single partially debonded ber 40 in an in nite matrix of [9, 10, 11], numerical studies followed. In the 1990's, Par s and collaborators [12] developed a Boundary Element Method (BEM) with the use of discontinuous singular elements at the crack tip and the Virtual Crack Closure Integral (VCCI) [13] for the evaluation of the Energy Release Rate (ERR). They validated their results [12] with respect to Toya's analytical 45 solution [11] and analyzed the e ect of BEM interface discretization on the stress eld in the neighborhood of the crack tip [14]. Following Comninou's work on the straight crack [8], they furthermore recognized the importance of contact to retrieve a physical solution avoiding interpenetration [12] and studied the e ect of the contact zone on debond ERR [15]. Their algorithm was then 50 applied to investigate the ber-matrix interface crack under di erent geometrical con gurations and mechanical loadings [16, 17, 18, 19, 20, 21, 22]. Recently the Finite Element Method (FEM) was also applied to the solution of the ber-matrix interface crack problem [23, 24, 25], in conjunction with the Virtual Crack Closure Technique (VCCT) [26, 27] for the evaluation of the 55 ERR at the crack tip. In [23], the authors validated their model with respect to the BEM results of [12], but no analysis of the e ect of the discretization in the crack tip neighborhood comparable to [14] was proposed. Thanks to the interest in evaluating the ERR of interlaminar delamination, di erent studies 1For the ber-matrix interface crack, aw size is measured in terms of the angle  sub- tended by half of the arc-crack, i.e. a= 2. 3 DRAFTexist in the literature on the e ect of mesh discretization on Mode I and Mode 60 II ERR of the bi-material interface crack when evaluated with the VCCT in the context of the FEM [28, 29, 30]. However, no comparable analysis can be found in the literature on mode separation and convergence analysis of the VCCT when applied to the ber-matrix interface crack (circular bi-material interface crack) problem in the context of a linear elastic FEM solution. In the present 65 article, we rst present the FEM formulation of the problem, together with the main geometrical characteristics, material properties, boundary conditions and loading. We then propose a vectorial formulation of the VCCT and express Mode I and Mode II ERR in terms of FEM natural variables. With this tool, we derive an analytical estimate of the ERR convergence and compare it with 70 numerical results. 2. FEM formulation of the ber-matrix interface crack problem In order to investigate the ber-matrix interface crack problem, a 2-dimensional model of a single ber inserted in a rectangular matrix element is considered (see Figure 1). Total element length and height are respectively 2 LandL, whereL 75 is determined by the ber radius Rfand the ber volume fraction Vfby L=Rf 2r Vf: (3) The ber radius Rfis assumed to be equal to 1 m. This choice is not dictated by physical considerations but for simplicity. It is thus useful to re- mark that, in a linear elastic solution as the one considered in the present work, the ERR is proportional to the geometrical dimensions of the model and, 80 consequently, recalculation of the ERR for bers of any size requires a simple multiplication. As shown in Fig. 1, the debond is placed symmetrically with respect to thexaxis and its size is characterized by the angle  (which makes the full debond size equal to 2 and the full crack length equal to Rf2). A region 85  of unknown size appears at the crack tip for large debond sizes (at least 4 DRAFT xz ∆θL L L¯ux=¯εxL ¯ux=−¯εxL ΩfΩm Rf O∆ΦFigure 1: Schematic of the model with its main parameters. 6080), in which the crack faces are in contact with each other and free to slide. Frictionless contact is thus considered between the two crack faces to allow free sliding and avoid interpenetration. Symmetry with respect to the xaxis is applied on the lower boundary while the upper surface is left free. Kinematic 90 coupling on the x-displacement is applied along the left and right sides of the model in the form of a constant x-displacement"xL, which corresponds to transverse strain  "xequal to 1% in the results here presented. Table 1: Summary of the mechanical properties of ber and matrix. Estands for Young's modulus,for shear modulus and for Poisson's ratio. Material E[GPa ][GPa ][] Glass ber 70.0 29.2 0.2 Epoxy 3.5 1.25 0.4 The model problem is solved with the Finite Element Method (FEM) within the Abaqus environment, a commercial FEM software [31]. The model is meshed 95 with second order, 2D, plane strain triangular (CPE6) and rectangular (CPE8) elements. A regular mesh of rectangular elements with almost unitary aspect ratio is used at the crack tip. The angular size of an element in the crack 5 DRAFTtip neighborhood represents the main parameter of the numerical analysis. The crack faces are modeled as element-based surfaces and a small-sliding contact 100 pair interaction with no friction is imposed between them. The Mode I, Mode II and total Energy Release Rates (ERRs) (respectively referred to as GI,GIIand GTOT) are evaluated using the VCCT [27], implemented in a in-house Python routine. A glass ber-epoxy system is considered in the present work, and it is assumed that their response lies always in the linear elastic domain. The elastic 105 properties of glass ber and epoxy are reported in Table 1. 3. Vectorial formulation of the Virtual Crack Closure Technique (VCCT) In order to express the VCCT formulation of the ERR in terms of FEM variables, we need to introduce a few rotation matrices in order to represent the discretized representation (FE mesh) of a crack along a circular interface. The 110 position of the crack tip is characterized by the angular size of the crack (see Sec. 2 and Fig. 1 for reference) and the rotation corresponding to the crack tip reference frame is represented by the matrix Rde ned as R=2 4cos () sin () sin () cos ()3 5: (4) Nodes belonging to the elements sharing the crack tip are involved in the VCCT estimation of the ERR and it is assumed that, given a suciently ne 115 discretization, they are aligned with the crack propagation direction de ned at the crack tip. However, irrespectively of how small the elements in the crack tip neighborhood are, a misalignment always exists with respect to the assumed crack propagation direction (in the crack tip reference frame). This is measured by the matrices P(p), de ned as 120 P(p) =2 4cos 1 +1p m  sin 1 +1p m  sin 1 +1p m  cos 1 +1p m 3 5 (5) 6 DRAFTandQ (q), equal to Q (q) =2 4cosq1 m sinq1 m sinq1 m cosq1 m3 5; (6) respectively for the free and bonded nodes involved in the VCCT estimation. In Eqs. 5 and 6, is the angular size of an element in the crack tip neighborhood (see Sec. 2 and Fig. 1), mis the order of the element shape functions and p;qare indices referring to the nodes belonging respectively to free and bonded elements 125 sharing the crack tip. Introducing the permutation matrix P=2 40 1 1 03 5; (7) it is possible to express the derivatives of rotation matrices R,PandQ  with respect to their argument: @R @=DR;@P @= 1 +1p m DP;@Q  @=q1 mDQ :(8) By means of Eqs. 5 and 6, we can express the crack tip forces Fxy=2 4Fx Fy3 5 and crack displacements uxy=2 4ux uy3 5in the crack tip reference frame (where the 130 tangential direction correspond to the direction of crack propagation) while taking into account the misalignment to the nite discretization as Fr=Q RFxyur=PRuxy (9) whereFr=2 4Fr F3 5andur=2 4ur u3 5. The crack tip forces can be expressed as a function of the crack opening displacement as 135 Fxy=Kxyuxy+eFxy; (10) 7 DRAFTwhereKxyis in general a full matrix of the form Kxy=2 4KxxKxy KyxKyy3 5andeFxy represents the e ect of the rest of the FE solution through the remaining nodes of the elements attached to the crack tip. As such, the term eFxycan be expressed as a linear combination of the solution vector uNof nodal displacements of the formeKNuN. Equation 10 thus become 140 Fxy=Kxyuxy+eKNuN: (11) An exemplifying derivation of the relationships expressed in Equations 10 and 11 can be found in Appendix A. It is worthwhile to observe that another author [32] proposed a relationship of the form Fxy=Kxyuxy. However, in [32], this relationship is assumed a priori and manipulated to propose a revised ver- sion of the VCCT, based on the assumption that the matrix Kxyshould be 145 diagonal to provide physically-consistent fracture mode partitioning. On the other hand, in the present work we derive the relationships of Eqs. 10 and 11 from the formulation of the Finite Element Method. According to our deriva- tion, it seems correct that the matrix Kxyshould not in general be diagonal in order to take into account Poisson's e ect. In fact, a positive crack opening 150 displacement would cause a transverse displacement in the neighborhood of the crack tip. Given that material properties are di erent on the two sides of a bi-material interface, a net shear would be applied to the crack tip which would correspond to a net contribution to the crack tip force related to crack shear dis- placement. The analytical derivations presented in Appendix A con rm these 155 physical considerations. Based upon the work of Raju [33], we introduce the matrix Tpqto represent the weights needed in the VCCT to account for the use of singular elements. As already done previously, indices pandqrefer to nodes placed respectively on the free (crack face) and bonded side of the crack tip. Nodes are enumerated 160 so that the crack tip has always index 1, i.e. the higher the index the further the node is from the crack tip. Matrix Tpqhas always a size of dd, where 8 DRAFTd= 2 for a 2 Dproblem and d= 3 for a 3 Dproblem. An element Tpq(i;j) withi;j= 1;:::;d represents the weight to be assigned to the product of com- ponentiof the displacement extracted at node pwith component jof the force 165 extracted at node q. The expression of Tpqfor quadrilateral elements with or without singularity is reported in Appendix B. Notice that, given mis the order of the element shape functions, the element side has m+ 1 nodes and this represents the upper limit of indices pandq. By using matrix Tpq, it is possible to express the total ERR Gevaluated with 170 the VCCT as GTOT =1 2Rfm+1X p=1m+1X q=1Tr uT r;pTT pqFr;q : (12) Introducing the vector G=2 4GI GII3 5of fracture mode ERRs, Mode I and Mode II ERR evaluated with the VCCT can be expressed as G=1 2Rfm+1X p=1m+1X q=1Diag Fr;quT r;pTT pq ; (13) whereDiag () is the function that extracts the main diagonal of the input matrix as a column vector. Substituting Equations 9 and 11 in Equations 12 175 and 13, we can express the Mode I, Mode II and total Energy Release Rate as a function of the crack displacements and the FE solution (more details in Appendix A) as GTOT =1 2Rfm+1X p=1m+1X q=1Tr Q RKxy;quxy;quT xy;pRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Tr Q ReFxy;quT xy;pRT PT TT pq(14) 9 DRAFTand G=2 4GI GII3 5=1 2Rfm+1X p=1m+1X q=1Diag Q RKxy;quxy;quT xy;pRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Diag Q ReKN;quNuT xy;pRT PT TT pq (15) 4. Rotational invariance of GTOT 180 Recalling Equation 14 and observing that matrix Tpqis always equal to the identity matrix pre-multiplied by a suitable real constant (see Eq. B.1 in Appendix B), the total Energy Release Rate can be rewritten as GTOT =1 2Rfm+1X p=1m+1X q=1Tr Q R Kxy;quxy;q+eFxy;q uT xy;pTT pqRT PT  = =1 2Rfm+1X p=1m+1X q=1Tr Q RFxy;quT xy;pTT pqRT PT  ; (16) whereFxyanduxyare the vectors of respectively the crack tip forces and crack displacements in the global ( xy) reference frame. Given that Q ,P185 andRall represent a linear transformation (a rigid rotation in particular), the invariance of the trace to linear transformations ensures that GTOT =1 2Rfm+1X p=1m+1X q=1Tr Q RFxy;quT xy;pTT pqRT PT  = =1 2Rfm+1X p=1m+1X q=1Tr Fxy;quT xy;pTT pq :(17) AsGTOT was de ned according to Equation 12 and given that Tr(AB) = 10 DRAFTTr(BA), it holds that GTOT =1 2Rfm+1X p=1m+1X q=1uT r;pTT pqFr;q=1 2Rfm+1X p=1m+1X q=1Tr Fr;quT r;pTT pq = =1 2Rfm+1X p=1m+1X q=1Tr Fxy;quT xy;pTT pq =1 2Rfm+1X p=1m+1X q=1uT xy;pTT pqFxy;q (18) which shows that the total Energy Release Rate is invariant to rigid rotations 190 and can be calculated equivalently with forces and displacements expressed in the local crack tip reference frame or the global reference frame. The analytical result is con rmed by the numerical solution of the ber-matrix interface crack with di erent element orders and model ber volume fractions, as shown in Figure 2. 195 The result of Equation 18 has also physical implications: { given that stress and displacement elds at the crack tip are the same, two cracks with di erent crack paths are energetically equivalent with respect to the total Energy Release Rate; { given that laws of the type GTOTGcgovern crack propagation, if Gc 200 do not depend on mode ratio, crack orientation will not a ect its growth. 11 DRAFT 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 J-integral VCCT, rotated VCCT, unrotated(a)Vf= 0:1%, 1storder elements, = 0:05. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 J-integral VCCT, rotated VCCT, unrotated(b)Vf= 0:1%, 2ndorder elements, = 0:05. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.911.11.21.3 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 J-integral VCCT, rotated VCCT, unrotated (c)Vf= 40%, 1storder elements, = 0:05. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.911.11.21.3 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 J-integral VCCT, rotated VCCT, unrotated(d)Vf= 40%, 2ndorder elements, = 0:05. Figure 2: Numerical invariance of the total Energy Release Rate: GTOT computed with the VCCT with rotated forces and displacements (label rotated ), with the VCCT with forces and displacements in the global reference frame (label unrotated ) and with J-integral method (label J-integral ). 5. Convergence analysis 5.1. Analytical considerations Substituting Equations 8 in the derivative of Equation 13, we can investigate the dependency of Mode I and Mode II ERR with respect to the size of an 205 12 DRAFTelement in the crack tip neighborhood through @G @=1 2Rf2m+1X p=1m+1X q=1Diag Q RKxyuxyuT xyRT PT TT pq 1 2Rf2m+1X p=1m+1X q=1Diag Q ReKNuNuT xyRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Diag Q RKxyuxyuT xyRT PT DTTT pq +1 2Rfm+1X p=1m+1X q=1Diag Q ReKNuNuT xyRT PT DTTT pq + +1 2Rfm+1X p=1m+1X q=1Diag DQ RKxyuxyuT xyRT PT TT pq +1 2Rfm+1X p=1m+1X q=1Diag DQ ReKNuNuT xyRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Diag Q RKxy@uxy @uT xyRT PT TT pq +1 2Rfm+1X p=1m+1X q=1Diag Q ReKN@uN @uT xyRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Diag Q RKxyuxy@uT xy @RT PT TT pq! +1 2Rfm+1X p=1m+1X q=1Diag Q ReKNuN@uT xy @RT PT TT pq! ; (19) which, after refactoring, provides @G @=1 G+1 2Rfm+1X p=1m+1X q=1Diag Q R Kxyuxy+eKNuN uT xyRT PT DTTT pq + +1 2Rfm+1X p=1m+1X q=1Diag DQ R Kxyuxy+eKNuN uT xyRT PT TT pq + +1 Rfm+1X p=1m+1X q=1Diag Q RKxy@uxy @uT xyRT PT TT pq +1 2Rfm+1X p=1m+1X q=1Diag Q ReKN@uN @uT xyRT PT TT pq + +1 2Rfm+1X p=1m+1X q=1Diag Q ReKNuN@uT xy @RT PT TT pq! : (20) Following the asymptotic analysis of [3, 1], in the case of an open crack the displacement in the crack tip neighborhood will have a functional form of the type 210 u()p (sin;cos) (log) with=1 2log1 1 +  (21) and is Dundurs' parameter introduced in Section 1. Application of Equa- tion 21 to the terms on the right hand side of Eq. 20 provides: uxy;uNu()p (sin;cos) (log)!0! 0; (22) uxyuT xy;uNuT xyu2() sin2;cos2;sincos (log)!0! 0; (23) 13 DRAFT@uxy @uT xy;@uN @uT xy1 2 sin2;cos2;sincos (log)+ sin2;cos2;sincos (log)!0!finite ; (24) G1 uxyuT xy1 u2() sin2;cos2;sincos (log)!0!finite: (25) In Equations 22-25, the multiplication by a trigonometric function of the type  sin;cos;sin2;cos2;sincos prevents the divergence of the asymptote. Recalling Eqs. 5 and 6, in the limit of !0 the rotation matrices become equal to the 215 identity matrix: P;Q !0!2 41 0 0 13 5: (26) Applying the results of Equations 22-26 to Eq. 20, it can be shown that the derivative of Gcan be split in a factor that goes to 0 in the limit of !0 and in a factor independent of : lim !0@G @1   *0 F() +C! : (27) Thus, asymptotically, the Mode I and Mode II Energy Release Rate be- 220 have like the logarithm of the angular size of the elements in the crack tip neighborhood: lim !0@G @1 R d! lim !0GAlog() +B: (28) 5.2. Numerical results Evaluations of the Mode I, Mode II and total Energy Release Rate using the VCCT applied to the FE solution of the ber-matrix interface crack in the 225 single ber model of Sec. 2 are reported respectively in Fig. 3, Fig. 4 and Fig. 5. Results for Mode I ERR in Fig. 3 show clearly the transition from the open crack regime, where Mode I ERR is di erent from zero, to the closed crack 14 DRAFT 0 10 20 30 40 50 60 70 80 9005·10−20.10.150.20.250.30.35 ∆θ/bracketleftBigg◦/bracketrightBiggGI J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦(a)Vf= 0:1%, 1storder elements. 0 10 20 30 40 50 60 70 80 9005·10−20.10.150.20.250.30.35 ∆θ/bracketleftBigg◦/bracketrightBiggGI J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (b)Vf= 0:1%, 2ndorder elements. 0 10 20 30 40 50 60 70 80 9000.10.20.30.40.50.60.70.80.91 ∆θ/bracketleftBigg◦/bracketrightBiggGI J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (c)Vf= 40%, 1storder elements. 0 10 20 30 40 50 60 70 80 9000.10.20.30.40.50.60.70.80.91 ∆θ/bracketleftBigg◦/bracketrightBiggGI J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (d)Vf= 40%, 2ndorder elements. Figure 3: E ect of the size of an element at the crack tip on Mode I ERR. regime of the debond, where GI= 0. Looking at Fig. 3, the crack is open for 60and it is closed , i.e. a contact zone is present, for  70. As 230 expected from the analysis of the previous section, and given that Mode I ERR is di erent from zero only in the open crack regime, a signi cant dependence on the element size can be observed in Fig. 3 when using both 1stand 2nd order elements and with both an e ectively in nite ( Vf= 0:1%) and nite size (Vf= 40%) matrix. At rst sight, it is immediate to see from Fig. 3 that a 235 decrease in leads to a decrease in GI. However, two further e ects can be observed due to the re nement of the mesh at the crack tip, i.e. the decrease of the element size . First, the occurrence of the peak GIis shifted to lower angles for very low volume fractions: it occurs at  = 30with= 1:0;0:5and at 15 DRAFT= 20with0:25for both 1stand 2ndorder elements and Vf= 0:1%. 240 Second, the apperance of the contact zone, i.e. the switch to the closed crack regime, is anticipated to smaller debonds: it occurs at  = 70with0:2 and at = 60with <0:2for both 1stand 2ndorder elements and both Vf= 0:1% andVf= 40%. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGII J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (a)Vf= 0:1%, 1storder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGII J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (b)Vf= 0:1%, 2ndorder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGII J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (c)Vf= 40%, 1storder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGII J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (d)Vf= 40%, 2ndorder elements. Figure 4: E ect of the size of an element at the crack tip on Mode II ERR. Observing Figure 4, it possible to notice the existence of two distinct regimes 245 in the behavior of GIIwith respect to the element size . For 60GII depends on the value of , while 70it is e ectively independent of the element size at the crack tip for both 1stand 2ndorder elements and both an e ectively in nite ( Vf= 0:1%) and nite size ( Vf= 40%) matrix. Comparing the value of  at which the change from the -dependency regime to the - 250 16 DRAFTindependency regime occurs for GIIwith Mode I ERR in Fig. 3, it is possible to observe that the -dependency regime change of Mode II ERR coincides with the onset of the contact zone, i.e. the transition from open crack regime to the closed crack regime. The result con rms the analytical considerations of the previous section: for an open crack both Mode I and Mode II ERR depend on 255 the element size at the crack tip. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (a)Vf= 0:1%, 1storder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.9 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (b)Vf= 0:1%, 2ndorder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.911.11.21.3 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (c)Vf= 40%, 1storder elements. 010 30 50 70 90 110 130 15000.10.20.30.40.50.60.70.80.911.11.21.3 ∆θ/bracketleftBigg◦/bracketrightBiggGTOT J/m2 1.0◦ 0.5◦ 0.25◦ 0.2◦ 0.125◦ 0.1◦ 0.0625◦ 0.05◦ (d)Vf= 40%, 2ndorder elements. Figure 5: E ect of the size of an element at the crack tip on total ERR. Further observation of Figure 4 reveals that, in the open crack regime, de- creasing the element size causes an increase of Mode II ERR. Similarly to Mode I ERR, a shift of the peak GIIcan also observed for Vf= 0:1%: the maximum value of GIIoccurs at  = 70for>0:25for 1storder elements 260 and for > 0:5for 2ndorder elements, while it is shifted to  = 60for 17 DRAFT0:25for 1storder elements and for 0:5for 2ndorder elements. Table 2: Summary of linear regression results and main statistical tests for Mode I ERR Vf[%] Order  []AJ m2 BJ m2 r[]r2[]p(A) []p(B) [] 0.1 1 10.0 0.0064 0.2113 0.9933 0.9866 7.48E-07 3.49E-14 20.0 0.0183 0.3331 0.9996 0.9992 1.44E-10 2.40E-16 30.0 0.0280 0.3392 1.0000 1.0000 2.25E-16 4.26E-21 40.0 0.0304 0.2524 0.9997 0.9995 4.38E-11 7.94E-15 50.0 0.0235 0.1278 0.9985 0.9970 8.61E-09 2.01E-11 60.0 0.0094 0.0284 0.9854 0.9709 7.75E-06 6.14E-07 0.1 2 10.0 0.0069 0.2103 0.9962 0.9924 1.36E-07 1.03E-14 20.0 0.0187 0.3277 0.9997 0.9994 7.85E-11 1.62E-16 30.0 0.0280 0.3296 1.0000 1.0000 3.28E-16 7.29E-21 40.0 0.0298 0.2408 0.9997 0.9995 4.82E-11 1.04E-14 50.0 0.0225 0.1177 0.9984 0.9967 1.10E-08 3.27E-11 60.0 0.0081 0.0228 0.9811 0.9626 1.66E-05 2.17E-06 40 1 10.0 0.0311 0.9196 0.9963 0.9927 1.03E-07 9.33E-15 20.0 0.0501 0.8882 1.0000 0.9999 1.21E-13 2.33E-19 30.0 0.0510 0.6374 0.9998 0.9996 1.66E-11 2.58E-16 40.0 0.0419 0.3760 0.9988 0.9976 4.56E-09 5.25E-13 50.0 0.0279 0.1713 0.9980 0.9961 2.22E-08 2.52E-11 60.0 0.0108 0.0391 0.9901 0.9804 3.44E-06 9.46E-08 40 2 10.0 0.0336 0.9148 0.9988 0.9977 3.45E-09 5.09E-16 20.0 0.0504 0.8719 1.0000 1.0000 3.70E-14 8.26E-20 30.0 0.0506 0.6191 0.9999 0.9997 7.63E-12 1.35E-16 40.0 0.0414 0.3608 0.9994 0.9989 4.95E-10 6.80E-14 50.0 0.0269 0.1593 0.9982 0.9964 1.66E-08 2.31E-11 60.0 0.0097 0.0329 0.9890 0.9781 4.96E-06 1.99E-07 Analysis of the total ERR in Figure 5 leads to an observation that was 18 DRAFTnot predicted by the considerations of the previous section: GTOT is e ectively independent of the element size in both the open and the closed crack regimes, 265 at least for reasonably small elements ( 1:0). Given that GII=GTOT for theclosed crack, it explains the independency of GIIfromafter the onset of the contact zone. Analysis of Fig. 3, Fig. 4 and Fig. 5 has shown the dependency of Mode I and Mode II ERR on the element size . Following the derivations of the previous 270 section, we model the dependency of GIandGIIwith respect to as G()=A() ln+B(); (29) whereA() andB() are parameters dependent on  estimated through linear regression (with x= ln) of the numerical results. 00.10.20.30.40.50.60.70.80.911.105·10−20.10.150.20.250.30.35 0.0064 lnδ+0.21130.0183 lnδ+0.3331 0.028lnδ+0.3392 0.0304 lnδ+0.2524 0.0235 lnδ+0.1278 0.0094 lnδ+0.0284 δ/bracketleftBigg◦/bracketrightBiggGI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (a) 1storder elements. 10−110005·10−20.10.150.20.250.30.35 0.0064 lnδ+0.21130.0183 lnδ+0.3331 0.028lnδ+0.3392 0.0304 lnδ+0.2524 0.0235 lnδ+0.1278 0.0094 lnδ+0.0284 lnδ[−]GI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (b) 1storder elements. 00.10.20.30.40.50.60.70.80.911.105·10−20.10.150.20.250.30.35 0.0069 lnδ+0.21030.0187 lnδ+0.3277 0.028lnδ+0.3296 0.0298 lnδ+0.2408 0.0225 lnδ+0.1177 0.0081 lnδ+0.0228 δ/bracketleftBigg◦/bracketrightBiggGI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (c) 2ndorder elements. 10−110005·10−20.10.150.20.250.30.35 0.0069 lnδ+0.21030.0187 lnδ+0.3277 0.028lnδ+0.3296 0.0298 lnδ+0.2408 0.0225 lnδ+0.1177 0.0081 lnδ+0.0228 lnδ[−]GI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (d) 2ndorder elements. Figure 6: Logarithmic dependence on of Mode I ERR: interpolation of numerical results for Vf= 0:1%. 19 DRAFTTable 3: Summary of linear regression results and main statistical tests for Mode II ERR Vf[%] Order  []AJ m2 BJ m2 r[]r2[]p(A) []p(B) [] 0.1 1.0 10.0 -0.0076 0.0228 -0.9996 0.9991 2.09E-10 1.64E-11 20.0 -0.0194 0.1211 -1.0000 1.0000 1.99E-15 2.02E-18 30.0 -0.0290 0.3007 -0.9999 0.9998 4.12E-12 1.97E-16 40.0 -0.0311 0.5270 -0.9995 0.9989 4.13E-10 1.05E-15 50.0 -0.0240 0.7375 -0.9979 0.9958 2.32E-08 1.66E-15 60.0 -0.0095 0.8685 -0.9835 0.9672 1.12E-05 1.22E-15 0.1 2.0 10.0 -0.0078 0.0249 -0.9996 0.9992 1.91E-10 1.06E-11 20.0 -0.0196 0.1272 -1.0000 1.0000 3.48E-15 2.78E-18 30.0 -0.0288 0.3108 -0.9999 0.9998 1.45E-12 5.47E-17 40.0 -0.0305 0.5387 -0.9995 0.9990 3.32E-10 6.55E-16 50.0 -0.0229 0.7478 -0.9979 0.9959 2.17E-08 1.09E-15 60.0 -0.0082 0.8744 -0.9806 0.9615 1.81E-05 8.26E-16 40.0 1.0 10.0 -0.0344 0.1055 -0.9997 0.9995 3.82E-11 2.73E-12 20.0 -0.0500 0.2977 -1.0000 0.9999 4.22E-14 5.66E-17 30.0 -0.0505 0.4866 -0.9999 0.9997 6.44E-12 4.82E-16 40.0 -0.0420 0.6454 -0.9996 0.9991 2.12E-10 9.66E-16 50.0 -0.0275 0.7386 -0.9985 0.9971 9.01E-09 1.44E-15 60.0 -0.0099 0.7402 -0.9926 0.9853 1.41E-06 5.13E-16 40.0 2.0 10.0 -0.0353 0.1145 -0.9998 0.9995 2.92E-11 1.50E-12 20.0 -0.0504 0.3130 -1.0000 0.9999 4.00E-14 4.17E-17 30.0 -0.0502 0.5039 -0.9999 0.9998 2.87E-12 1.69E-16 40.0 -0.0410 0.6615 -0.9996 0.9992 2.02E-10 6.89E-16 50.0 -0.0263 0.7502 -0.9987 0.9973 6.87E-09 7.76E-16 60.0 -0.0090 0.7458 -0.9921 0.9842 1.79E-06 3.37E-16 As shown in Fig. 6, Fig. 7, Fig. 8 and Fig. 9 both in linear and logarithmic scales of, the result is remarkable: both the correlation coecient rand the 275 20 DRAFT 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.91 0.0311 lnδ+0.9196 0.0501 lnδ+0.8882 0.0510 lnδ+0.6374 0.0419 lnδ+0.3760 0.0279 lnδ+0.1713 0.0108 lnδ+0.0391 δ/bracketleftBigg◦/bracketrightBiggGI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦(a) 1storder elements. 10−110000.10.20.30.40.50.60.70.80.91 0.0311 lnδ+0.9196 0.0501 lnδ+0.8882 0.0510 lnδ+0.6374 0.0419 lnδ+0.3760 0.0279 lnδ+0.1713 0.0108 lnδ+0.0391 lnδ[−]GI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (b) 1storder elements. 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.91 0.0336 lnδ+0.9148 0.0504 lnδ+0.8719 0.0506 lnδ+0.6191 0.0414 lnδ+0.3608 0.0269 lnδ+0.1593 0.0097 lnδ+0.0329 δ/bracketleftBigg◦/bracketrightBiggGI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (c) 2ndorder elements. 10−110000.10.20.30.40.50.60.70.80.91 0.0336 lnδ+0.9148 0.0504 lnδ+0.8719 0.0506 lnδ+0.6191 0.0414 lnδ+0.3608 0.0269 lnδ+0.1593 0.0097 lnδ+0.0329 lnδ[−]GI J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (d) 2ndorder elements. Figure 7: Logarithmic dependence on of Mode I ERR: interpolation of numerical results for Vf= 40%. r2ratio (of explained to total variance) are always greater than 0 :95 and the p-values of the coecients AandBare at least <1E6 and often <1E11 (see Table 2 for GIand Table 3 for GII). The results of the linear regression con rm the analytical derivations of the previous section, which showed the logarithmic behavior of Mode I and Mode II ERR. Similar conclusions were 280 reached in [28, 29] for a straight bi-material crack with respect to the parameter a=a; however, no functional expression of G()was proposed. 21 DRAFT 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.91 −0.0076 lnδ+0.0228−0.0194 lnδ+0.1211−0.0290 lnδ+0.3007−0.0311 lnδ+0.5270−0.0240 lnδ+0.7375−0.0095 lnδ+0.8685 δ/bracketleftBigg◦/bracketrightBiggGII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦(a) 1storder elements. 10−110000.10.20.30.40.50.60.70.80.9 −0.0076 lnδ+0.0228−0.0194 lnδ+0.1211−0.0290 lnδ+0.3007−0.0311 lnδ+0.5270−0.0240 lnδ+0.7375−0.0095 lnδ+0.8685 lnδ[−]GII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (b) 1storder elements. 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.91 −0.0078 lnδ+0.0249−0.0195 lnδ+0.1272−0.0288 lnδ+0.3108−0.0305 lnδ+0.5387−0.0229 lnδ+0.7478−0.0082 lnδ+0.8744 δ/bracketleftBigg◦/bracketrightBiggGII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (c) 2ndorder elements. 10−110000.10.20.30.40.50.60.70.80.9 −0.0078 lnδ+0.0249−0.0195 lnδ+0.1272−0.0288 lnδ+0.3108−0.0305 lnδ+0.5387−0.0229 lnδ+0.7478−0.0082 lnδ+0.8744 lnδ[−]GII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (d) 2ndorder elements. Figure 8: Logarithmic dependence on of Mode II ERR: interpolation of numerical results forVf= 0:1%. 22 DRAFT 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.9 −0.0344 lnδ+0.1055−0.0500 lnδ+0.2977−0.0505 lnδ+0.4866−0.0420 lnδ+0.6454−0.0275 lnδ+0.7386 −0.0099 lnδ+0.7402 δ/bracketleftBigg◦/bracketrightBiggGII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦(a) 1storder elements. 10−110000.10.20.30.40.50.60.70.80.9 −0.0344 lnδ+0.1055−0.0500 lnδ+0.2977−0.0505 lnδ+0.4866−0.0420 lnδ+0.6454−0.0275 lnδ+0.7386 −0.0099 lnδ+0.7402 lnδ[−]GII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (b) 1storder elements. 00.10.20.30.40.50.60.70.80.911.100.10.20.30.40.50.60.70.80.9 −0.0353 lnδ+0.1145−0.0504 lnδ+0.3130−0.0502 lnδ+0.5039−0.0410 lnδ+0.6615−0.0263 lnδ+0.7502 −0.0090 lnδ+0.7458 δ/bracketleftBigg◦/bracketrightBiggGII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (c) 2ndorder elements. 10−110000.10.20.30.40.50.60.70.80.9 −0.0353 lnδ+0.1145−0.0504 lnδ+0.3130−0.0502 lnδ+0.5039−0.0410 lnδ+0.6615−0.0263 lnδ+0.7502 −0.0090 lnδ+0.7458 lnδ[−]GII J/m2 10.0◦ 20.0◦ 30.0◦ 40.0◦ 50.0◦ 60.0◦ (d) 2ndorder elements. Figure 9: Logarithmic dependence on of Mode II ERR: interpolation of numerical results forVf= 40%. 23 DRAFT6. Conclusions & Outlook The application of the Virtual Crack Closure Technique to the calculation of Mode I, Mode II and total Energy Release Rate was analyzed in the context 285 of the Finite Element solution of the bi-material circular arc crack, or ber- matrix interface crack. A synthetic vectorial formulation of the VCCT has been proposed and its usefulness exempli ed in the analysis of the mesh dependency. By both analytical considerations and numerical simulations, it has been shown that: 290 the total ERR is invariant to rotations of the reference frame (and more in general to linear transformations), which implies that rotation of crack tip forces and displacement is actually not required in the use of the VCCT for the calculation of GTOT; the total ERR does not depend on the size of the elements at the crack 295 tip, at least for reasonably small elements ( 1:0) ; as a consequence, Mode II ERR for the closed interface crack does not depend on, asGII=GTOT after the onset of the contact zone; for the open interface crack, Mode I and Mode II ERR depend on the element size through a logarithmic law of the type A() ln+B(); 300 the sign of the logarithm is always positive for GI, i.e. it decreases when decreases, and negative for GII, i.e. it increases when decreases. The conclusion is signi cant: as the behavior of Mode I and Mode II is loga- rithmic with respect to mesh size, there exists no asymptotic limit and thus no convergence of the values. A convergence analysis based on the reduction of the 305 error between successive iterations would not provide a reliable assessment of the accuracy of the FE solution of Mode I and Mode II Energy Release Rate of the ber-matrix interface crack. A validation is thus required with respect to data obtained through a di erent method, be it analytical, numerical or experimen- tal. Moreover, it has been shown that: rst, the same behavior appears when 310 24 DRAFTusing 1stas well as 2ndorder elements; second, no improvement is expected with the use of singular elements, as the logarithmic dependency of GIandGIIis governed by the de nition of ERR itself together with the asymptotic behavior of the displacement eld at the crack tip. These two conclusions run contrary to the suggestions provided in the manuals of many commercial FEM packages, 315 such as Abaqus [31] which suggests that (Section 11.4.2 of the Abaqus Analysis User's Guide ): \Sharp cracks (where the crack faces lie on top of one another in the undeformed con guration) are usually modeled using small-strain assump- tions. Focused meshes, [...], should normally be used for small-strain fracture mechanics evaluations. However, for a sharp crack the strain eld becomes sin- 320 gular at the crack tip. [...] In most cases the singularity at the crack tip should be considered in small-strain analysis (when geometric nonlinearities are ignored). Including the singularity often improves the accuracy of the J-integral, the stress intensity factors, and the stress and strain calculations because the stresses and strains in the region close to the crack tip are more accurate. ". We have shown 325 that, in the context of the ber/matrix interface crack, the convergence of the Energy Release Rate is determined by the asymptotic behavior of the elastic solution and only marginally by the choice of element order and type, thus contradicting the statements in [31]. Acknowledgements 330 Luca Di Stasio thanks Prof. Janis Varna for the useful discussions and sug- gestions. Luca Di Stasio gratefully acknowledges the support of the European School of Materials (EUSMAT) through the DocMASE Doctoral Programme and the European Commission through the Erasmus Mundus Programme.