An extended finite element method for modeling crack growth with frictional contact

Abstract A new technique for the finite element modeling of crack growth with frictional contact on the crack faces is presented. The eXtended Finite Element Method (X-FEM) is used to discretize the equations, allowing for the modeling of cracks whose geometry are independent of the finite element mesh. This method greatly facilitates the simulation of a growing crack, as no remeshing of the domain is required. The conditions which describe frictional contact are formulated as a non-smooth constitutive law on the interface formed by the crack faces, and the iterative scheme implemented in the LATIN method [Nonlinear Computational Structural Mechanics, Springer, New York, 1998] is applied to resolve the nonlinear boundary value problem. The essential features of the iterative strategy and the X-FEM are reviewed, and the modifications necessary to integrate the constitutive law on the interface are presented. Several benchmark problems are solved to illustrate the robustness of the method and to examine convergence. The method is then applied to simulate crack growth when there is frictional contact on the crack faces, and the results are compared to both analytical and experimental results.


Introduction
Modeling crack growth with frictional contact presents signi®cant challenges for traditional ®nite element methods. In addition to the standard complications presented by remeshing, a possible mismatch in element geometries across the crack faces creates additional burdens for imposing the contact-friction law. In this paper, we present a framework to model frictional contact on arbitrary evolving interfaces and cracks with the ®nite element method. The nonlinear boundary value problem is resolved by adapting the iterative scheme employed by the LATIN method to the problem of interest. The discrete formulation employs the recently developed eXtended Finite Element Method (X-FEM), which allows for the modeling of geometric features which are independent of the ®nite element mesh. The combination of the two methods is shown to provide a straightforward means of circumventing several of the diculties associated with standard ®nite element formulations for these classes of problems. The formulation is used to model crack growth when there is frictional contact on the discontinuity surface, but general constitutive laws are also considered.
Several instances of nonlinear constitutive laws on interfaces arise in classical mechanics. Unilateral contact with dry friction can be formulated as a non-smooth constitutive dissipative law, relating the displacements and tractions on the interface. Dierent numerical methods have been proposed to model 1 frictional contact, including penalty and mixed methods. In the present investigation, the iterative procedure adopted to solve the nonlinear problem on the interface is taken from the LATIN method [15], which separates the formulation into a linear global stage and a nonlinear local stage. The method is quite general, in that it has also been used to solve traditional problems in plasticity and viscoplasticity, as well as domain decomposition [4]. In terms of its application to contact conditions on an interface, it shares similar features with the method of augmented Lagrangians (see [1,22]). The work presented in this paper is distinguished by the manner in which the problem is discretized and the application to frictional crack growth.
The modeling of moving discontinuities with the ®nite element method is cumbersome due to the need to update the mesh`topology' to match the geometry of the discontinuity. The X-FEM [2,18] circumvents these problems by enriching a standard mesh-based approximation with additional discontinuous functions. The form of the enriched approximation follows the partition of unity framework presented in [17]. The geometry of the discontinuity is thus evolved by updating the enrichment scheme; no remeshing of the domain is required. The only interaction between the mesh and the geometry of the discontinuity involves the construction/selection of the enriched basis functions.
The X-FEM shares common features with other recently developed ®nite element techniques to model crack growth without remeshing. These include the incorporation of a discontinuous mode on an element level [23], a moving mesh technique [20], and most recently the element partition method (EPM) [9]. These methods all possess an advantage over boundary element formulations in that the formulations can easily be extended to problems with material nonlinearity.
The present method is applied to model crack growth in a compressive ®eld, a situation which frequently arises in the study of geomaterials. A condition known as axial splitting often occurs in which the initial crack will turn and propagate in the direction of the applied load. The phenomenon has been studied experimentally by Hoek and Bieniawski [13] and Cotterell [5], or more recently by Nemat-Nasser and Horii [19]. For a single crack in a compressive ®eld with frictional contact on the crack faces, Nemat-Nasser and Horii [19] also gives analytical results. Finite element studies of crack growth in rock with continual remeshing can be found in [14], where contact on the crack faces was not considered. This paper is organized as follows. In Section 2, a review of the problem formulation and iterative procedure is given in the context of the LATIN framework. Section 3 presents the discretization scheme using the X-FEM, and discusses the modi®cations necessary to model the interface. In Section 4, the particular constitutive laws considered are reviewed, including unilateral contact with sliding and friction. For each constitutive law studied, convergence results for a benchmark problem are also provided. Numerical simulations of crack growth under compressive loads are presented in Section 5, and closing remarks are given in the last section.

Problem formulation
In this section, we review the nonlinear boundary value problem (BVP) of interest. Consider the domain X & R 2 containing an internal boundary described by the curve C d as shown in Fig. 1. The following discussion is applicable when C d represents a crack or when it separates the body into two distinct pieces, in which case it is better called an interface. Consequently, we refer to the curve C d as both the crack and the interface in the sections which follow. We also distinguish the crack faces C d and C À d as shown in Fig. 2, with the normal to C d denoted by n. We assume quasi-static loading by a body force b and tractions t imposed on the part C t of the boundary of X. The domain is constrained by prescribed displacements u imposed on the part C u of the boundary. The stress and displacement ®elds are denoted by r and u, respectively. We also introduce the displacement and traction on each face of the crack: w , t on C d and w À , t À on C À d . To refer to both w and w À , we use the compact notation w. Similarly, t denotes the pair (t , t À ).
The set of equations describing the problem originates from the kinematics, equilibrium and the constitutive laws. Assuming small strains and displacements, as well as a linear elastic material law, the constitutive law in the body X simply reads where the operator u is the symmetric part of the gradient of u and C is the Hooke tensor. The interfacial constitutive law on the crack faces is expressed in terms of the displacement and tractions on both sides of the crack: The precise form of the G operator depends on the particular constitutive law to be modeled. For unilateral contact with sliding or with friction, the operator is nonlinear. The speci®c equations for each constitutive law studied in this paper are provided in Section 4. We now describe the kinematics and equilibrium equations. The pair u; w is said to be kinematically admissible if u; w P U a , where where the space U is related to the regularity of the kinematic ®elds. The pair u; w is said to be kinematically admissible to zero if u; w P U 0 , where The set of statically admissible stress traction pairs r; t is denoted by T a : T a r; t P T satisfying  where we have used the condensed notation and the space T is related to the regularity of the static ®elds. The kinematic and equilibrium equations characterized by the spaces U a and T a are global in nature but linear. In contrast, the constitutive laws are local in nature, with (2) possibly nonlinear on C d . The LATIN method takes advantage of this distinction by decomposing the equations into a set A related to the linear, possibly global equations and a set I of local, possibly nonlinear equations. Using the compact notation: s u; w; r; t 7 and the space S gathering all such regular s: S fu; w; r; t : u; w P U; r; t P Tg; 8 we may de®ne the sets A and I as A fs P S : u; w P U a ; r; t P T a ; r Cu on Xg; 9 I fs P S : Gw ; w À ; t ; t À 0 on C d g: 10 Since the constitutive law in the body is linear (Hooke's law), we have placed it in the set A. In the case of nonlinear material behavior in the body, for instance plasticity, the constitutive equations should be placed in I.
The problem is to ®nd the element s AI which is located at the intersection of the sets A and I. The iterative strategy employed by the LATIN method to ®nd s AI is described in the following section.

Iterative procedure
The problem of ®nding the intersection of A and I is represented geometrically in Fig. 3. We use the superscript A and I to denote an element s in A or I, respectively. The iterative strategy begins with an initial element s A 0 in A and builds a sequence of approximate solutions s A 1 ; . . . ; s A n ; s A n1 ; . . . until convergence. A given iteration s A n 3 s A n1 involves two steps: s A n 3 s I n and s I n 3 s A n1 as shown in Fig. 3. The iterations stop when the`distance' between s I n and s A n1 (measured with an appropriate norm) is below a speci®ed tolerance.
Two successive approximations are always tied by a given search direction. The direction used when going`down' from the set A to the set I is denoted by E AI whereas the direction when going`up' from I to A is denoted by E IA . These two directions may be ®xed or may evolve with subsequent iterations. The appropriate choice of the search directions in order to achieve (fast) convergence is an important aspect of the LATIN method (see [15]). In the present investigation, the choice of search directions is consistent with the application of the LATIN method to domain decomposition and contact between structural components [4].
The next two sections describe the two steps in the method: · Going down from s A n to s I n in the search direction E AI . We will see that this step is local and may be solved independently for each pair of points facing each other on the crack faces. · Going up from s I n to s A n1 in the search direction E IA . This step involves a global solve.
2.1.1. The local step: update on the crack faces The local update involves determining a new estimate s I n for the tractions and displacements on the interface given quantities obtained from a solution s A n P A. The new approximation is required to satisfy the constitutive law (i.e., s I n P I). Additional equations are provided by the search direction. To move from an element s A n P A to s I n P I, the search direction E AI is associated with a linear operator k 0 . The search equations are then The above can be written more compactly as with the understanding that the relationship holds on both C d and C À d , separately. This compact form will be used in the remainder of this section when convenient.
The process of determining s I n given s A n then involves solving the above search equations in conjunction with the constitutive law on the interface s I n P I A Gw n ; w À n ; t n ; t À n 0 on C d : 13 The resulting closed form equations for t I n ; w I n with the constitutive laws studied in this paper are presented in Section 4.

The global step
We wish to move from s I n P I to s A n1 P A. This is accomplished through a global solve in conjunction with additional equations provided by the search direction E IA which is also associated with the k 0 operator on the interface. With a known element s I n , the search equations are then This equation is used in conjunction with the weak form of the equilibrium equation (5), by solving for t A n1 and making the substitution in the surface integrals on C d . The problem to be solved is then: Find u n1 ; w n1 P U a which satis®es and we have used the compact notation introduced in Eq. (6). After solving (15) and (16), the local tractions on the interface at step n 1 are given by

Convergence and error criterion
The search directions for the global and local steps, Eqs. (12) and (14), are conjugated, i.e., expressed up to a sign by the same linear operator k 0 . This ensures the convergence of the iterative strategy provided k 0 is symmetric positive de®nite and some conditions on the material operator G described in [15] are satis®ed. In this paper, we shall consider k 0 as an isotropic operator where I d is the identity matrix. The parameter k is strictly positive and has the dimension of a stiffness over a length. The extent to which this parameter k affects the rate of convergence to the solution s AI is examined in Section 4. The error indicator of convergence is computed after each global solve as the distance between the current and previous approximations: with the norm When the indicator g is below some tolerance, the iteration process is stopped. Note that if g 0, s A n1 s I n is the exact solution. Finally, concerning the initialization, the solution s A 0 is built by solving the global step, Eqs. (15) (17), with t I n w I n 0 and setting s A 0 u n1 ; w n1 ; r n1 ; t n1 .

Discretization
This section presents the application of discontinuous enrichment in the X-FEM to represent the crack C d , as well as the approximation for the quantities u, w and t. In addition, we present the integration schemes developed to construct element stiness matrices local to the crack and the surface integrals on C d .

The X-FEM approximation
Consider a typical ®nite element mesh and a crack on a rectangular domain as shown in Fig. 4, where the crack geometry C d is independent of the mesh. On the crack, we wish to satisfy the constitutive law (2). With the exception of perfect contact, the constitutive laws discussed in Section 4 allow for discontinuous displacements on C d . The X-FEM enriches a standard approximation locally about the crack with discontinuous functions. The approximation to the displacement ®eld ux takes the form where the notations used are summarized below (see [18] for more details) · I is the set of all nodes in the mesh. · u i is the classical (vectorial) degrees of freedom at node i. · / i is the shape function associated with node i. Each shape function / i has a compact support x i given by the union of the element domains connected to node i. · J & I is the subset of nodes that are enriched for the crack discontinuity and b j are the corresponding additional degrees of freedom. The nodes in J are such that their support (we mean the support of the nodal shape function) intersects the crack but do not cover any crack tips. For example, the circled nodes in Fig. 4 belong to J. · K 1 & I and K 2 & I are the subset of nodes that are enriched for the ®rst and second crack tips, respectively. The corresponding additional degrees of freedom are c l1 k and c l2 k , l 1; . . . ; 4 as the near-tip regions are enriched with four dierent crack functions. The nodes in K 1 (K 2 ) are such that their support contain the ®rst (second) crack tip. Fig. 4 represents the nodes in K 1 and K 2 by squares and diamonds, respectively. Finally, note that if the crack was to cut the domain into two pieces (i.e., no tips) the sets K 1 and K 2 would be empty. · The function H x is discontinuous across C d whereas the functions F 1 l r; h and F 2 l r; h, l 1; . . . ; 4 are the sets of near-tip functions for the tip 1 and 2, respectively. They are described more precisely below. The function H x is a`generalized Heaviside' function, in which the discontinuity is aligned with the surface C d . The surface C d is considered to be a curve parametrized by the curvilinear coordinate s, as in Fig. 5. Given a point x in the domain, we denote by x Ã the closest point on C d to x. At x Ã , we construct the tangential and normal vector to the curve, e s and e n , with the orientation of e n taken such that e s Â e n e z , where the vector e z points out of the page. The function H x is then given by the sign of the scalar product x À x Ã Á e n , i.e., À1 for x À x Ã Á e n < 0: In the case of a kinked crack as shown in Fig. 5(b), where no unique normal but a cone of normals is de®ned at x Ã , H x 1 if the vector x À x Ã belongs to the cone of normals at x Ã and )1 otherwise. The set of near-tip functions F l r; h are a set of functions which span the exact asymptotic crack-tip ®elds for a linear elastic material: where r; h are the local polar coordinates at the crack tip (see [12]). Note that the ®rst function in (24), r p sin h 2 , is discontinuous across the crack faces whereas the last three functions are continuous. It bears emphasis that the near-tip discontinuity can be represented with other sets of functions, or even a single function which is discontinuous across the crack tip geometry (see [8]). In the context of linear elastic fracture mechanics, the incorporation of the exact near-tip asymptotic functions signi®cantly improves the accuracy of the formulation.
Multiple cracks can be treated in the above framework, by incorporating additional discontinuous and near-tip enrichment. A simulation of multiple cracks growing in a compressive ®eld is given in Section 5. Some special care must be taken for intersecting and branching cracks to represent all possible displacement modes at the junction. The modes can be represented by enriching with a junction function [6] or with products of discontinuous functions H x [7].
The approximation for the displacement ®eld is then substituted into the linear system (15), yielding a system of linear algebraic equations, which can be written as where K C is the stiness matrix coming from the ®rst term in (15), K I is the stiness coming from the second term in (15), d is the vector of nodal unknowns, and f is the vector of nodal forces. The vector of nodal unknowns gathers the degrees of freedom fu i ; b j ; c l1 k ; c l2 k g. We note that in this method the term K C K I is a constant matrix which only needs to be assembled and factored at the ®rst iteration.

Numerical integration
Since the geometry of the crack is independent of the ®nite element mesh, it is necessary to modify the quadrature routines for the domain and line integrals in (15). In this section, we describe the integration schemes implemented in the present investigation. We also provide the approximations for the traction and displacement ®elds t; w on the crack in relation to the construction of the line integrals on C d .
For elements cut by the crack and enriched with the jump function H x, it is useful to modify the standard element quadrature strategy to account for the discontinuity. As the crack is allowed to be arbitrarily oriented in an element, standard Gauss quadrature may not be adequate. If the integral of the function H x is indistinguishable from that of a constant function, the system of equations may be singular.
The discrete weak form is normally constructed with a loop over all elements, as the domain is approximated by where m is the number of elements, and X e is the element subdomain. For elements cut by a crack, we de®ne the element subdomain to be a union of a set of subpolygons whose boundaries align with the crack geometry: X e mes e1 X es ; 27 where m es denotes the number of subpolygons for the element. The subtriangles shown in Fig. 6 work well for integration across the discontinuity C d . These subtriangles are also used to select nodes enriched with the discontinuous function H x, by calculating the percentage of a node's support above and below the crack. Simpler schemes, such as the trapezoids used by Fish [11] and the formulation suggested in [7] may be adequate. It is emphasized that the subpolygons are only necessary for quadrature; no additional degrees of freedom are associated with their construction. In the construction of the matrix equations, the element loop is replaced by a loop over the subpolygons for those elements cut by the crack.
To construct the integrals on the crack surface, it is necessary to discretize C d . In traditional ®nite element discretizations [3], Gauss integration points are placed on the element faces which align with the interface. Here, since the crack and mesh geometry are independent, we ®rst divide C d into one-dimensional segments. The segments are determined according to the interaction of the interface geometry with the mesh as illustrated in Fig. 7.
The curve C d is composed of a set of one-dimensional segments. For each segment, we determine all of the intersections with those elements cut by the interface, which results in a set of one-dimensional subelements along C d . These subelements are similar to the two-dimensional subpolygons used for the volume integrals, in that no new displacement degrees of freedom are associated with their generation. We note that in the case when the interface aligns with the element edges, this strategy reduces to the traditional approach.
In order to numerically integrate the terms in (15) on C d , Gauss quadrature points are used along each of the one-dimensional subelements. The ®elds w and t are then interpolated discretely at the Gauss points on the crack, leading to a list of pairs t À g ; w À g on C À d and t g ; w g on C d , where n 1; . . . ; n g is the number of Gauss points on each side of the crack. This discretization is shown in Fig. 8. The iterative process is then a means of determining the coecients (u i , b j , c k , t g , w g ) in conjunction with the equilibrium equations and constitutive law on the interface. In the case of a perfect bond of the crack faces, for example, we expect the nodal coecients b j and c 1 k to be driven to zero by the iterative procedure such that the displacement ®eld is continuous across C d . We also expect t g Àt À g and w g w À g at the discrete points on the interface. These relationships for other constitutive laws are provided in the following section.

Benchmark problems
The two problems reported in this section are a plate with an interface and a crack shown in Figs. 9 and 10, respectively. In the examples which follow, the plate is discretized with a uniform grid of 20 Â 20 quadrilateral elements, and the rigid body modes are constrained. The interface is modeled in the mesh using discontinuous enrichment whereas the crack needs both discontinuous and near-tip enrichment. Unless otherwise noted, the material parameters are: E 100 MPa and m 0:3.
The benchmark problems serve two purposes: · They examine the performance of the LATIN method for this class of discretizations, i.e. those with discontinuous and near-tip enrichment. · They study the eect of the parameter k in (19) on the rate of convergence. Three dierent interfacial constitutive laws on the crack faces are considered: perfect contact and unilateral contact with or without friction.

Perfect contact
A perfect bond between C d and C À d corresponds to the following constraints: w w À ; t Àt À 28  for which the local update equations are given in Appendix B. This law on the crack faces is not very interesting since the structure behaves as if there was no crack. However, this case serves to illustrate several features of the iterative procedure.
Consider the interface problem shown in Fig. 9 with two angles, h 0°and h 20°, and k set to E=W , 2E=W and 10E=W . We note that when h 20°, the interface is not aligned with the mesh. The results shown in Fig. 11 indicate that whether or not the interface aligns with the mesh, there is little change in the convergence for a given k. Less than 1% error is obtained in 20 iterations for all 3 values of k and k 2E=W seems to be an effective choice to achieve fast convergence. These results show trends similar to those reported in the application of the LATIN method to domain decomposition [4].
A similar study is performed with a crack (Fig. 10) of half-length a W =5. Results are shown in Fig. 12 for k 0 E=W , and k 0 E=a, with the angles chosen as before. For this problem, the optimal choice of k  appears to be k E=a, as opposed to k E=W , though the difference is slight. In all calculations the error drops below 1% between 10 and 20 iterations.
We note that for this benchmark problem, the exact solution does not contain any discontinuous or singular ®elds. The iterative procedure therefore results in the enriched degrees of freedom for the generalized Heaviside and near-tip functions (24) to be driven towards zero. After the ®nal iteration, these degrees of freedom will not be exactly zero and some form of the singularity was observed in the stress ®elds near the crack tips. The eect of these singular stress ®elds is suciently small, however, as to result in computed stress intensity factors less than 1% of those obtained in the following section where unilateral contact is modeled on the crack faces.

Unilateral contact
We now consider unilateral contact with or without friction on the interface. It bears emphasis that the only change necessary from modeling perfect contact is in the local update equations.
The constitutive law is easily described by decomposing the traction and displacements on the crack faces in terms of their normal and tangential components. The normal components are related by a unilateral constraint: 29a t Á n 6 0; t À Á n P 0; t Á n Àt À Á n; 29b where n is the unit normal vector to C d . The tangential components must be in equilibrium: where P T is the tangential projection operator: and Â denotes the vector cross-product. The additional equations depend on whether or not there is friction. When the friction is idealized using a Coulomb law, the maximum frictional force supported on the interface is where l is the coecient of friction and the equations to be satis®ed are kP T tk 6 g; 33a where k P 0 is a positive constant.
In the case of unilateral contact without friction, Eq. (33a) are replaced by The local update equations for unilateral contact derived in [3] are recalled in Appendix B.
To examine the rate of convergence with unilateral contact, we consider the problem of Fig. 10 when a W =5 and h 20°. Fig. 13 shows the results for k E=W , k 2E=W , and k E=a. The results are similar to those of perfect contact for the same geometry, and indicate that k E=a provides the best rate of convergence. Contour plots of the displacement in the crack coordinate system for this choice of k are shown in Fig. 14 at iteration 20. It is clear that the displacement tangential to the crack is discontinuous, while the normal displacement is continuous across C d .
A relatively simple test to assess the performance of contact with friction is considered next. We wish to determine the coecient of friction l which prevents the top part from sliding, given the angle h in Fig. 9. This is a classical problem in mechanics, with the solution l tanh: 35 We consider the case when h 11:31°, for which l must be taken to be just above 0:2 to prevent slip on the interface. Fig. 15 shows the convergence results when l 0:202 for various values of k.
The interesting feature of the iterative method is that for the initial iterations, sliding is predicted kG I n k P g. It is only when s n becomes suciently close to s AI that the correct`stick' solution is achieved. This occurs at dierent stages depending on the choice of k, as indicated in the ®gure. For this particular benchmark problem, the method nevertheless requires several more iterations to obtain the correct converged solution. We note that the true solution lies very close to the`slip' threshold, and the large number of iterations required by the method is characteristic of similar contact algorithms that initially predict slip for near-slip' solutions [16].

Comparison with experiments
The study of crack growth in a compressive ®eld has a long history. See for example the early work by Hoek and Bieniawski [13] and Cotterell [5], or the more recent papers of Ingraea and Heuze [14] and Nemat-Nasser and Horii [19]. Crack growth in a compressive ®eld is observed in the fracture of geomaterials such as rocks under large geotectonic states of stress. A condition known as axial splitting often occurs in which the initial crack in an overall compressive ®eld will turn and propagate in the direction of the applied load. Nemat-Nasser and Horii [19] presented experimental results in which thin slits were cut in glass and resin plates. They found that cracks oriented at an angle to the principle compressive load almost invariably propagated in the direction of the load. If no lateral loads were applied, the fracture tended to be stable and the propagating cracks arrest if the magnitude of the applied load is not suciently increased.
The numerical simulation of fracture in a compressive ®eld was presented by Ingraea and Heuze [14]. In that study, the contact between the crack faces was not considered as the cracks emanated from a notch cut in the specimen. Finite elements with remeshing were then used to model the crack propagation. In the present investigation, we explicitly model frictional contact on the crack faces and compare the numerical results to those presented in [19]. 14 The geometry under consideration is shown in Fig. 10. The dimensions of the plate and initial crack are taken from Nemat-Nasser and Horii [19], with W 4 in., and a 0:4 in. A mesh of 13,508 quadrilateral elements is used, with re®nement near the crack. The material properties are taken to be that of PMMA: E 0:45 Â 10 6 psi and m 0:35.
We begin by considering a single crack with orientation h 0:7p. Nemat-Nasser and Horii [19] provides both analytical and experimental results for this situation. We consider the case of frictionless contact on the crack faces, and simulate quasi-static crack growth with Da 0:06 in for a total of 10 steps. The direction of crack growth is determined from the maximum hoop stress criterion [10]. An important consideration in determining the stress intensity factors is to account for the tractions on the crack faces (see Appendix A). The trajectory of the top crack tip is compared to the results of Nemat-Nasser and Horii [19] in Fig. 16. The crack is observed to grow in the direction of the applied load, and shows good areement with the experiments. We note that the numerical solution oscillates about the experimental observation, a phenomenon attributable to changes of sign in K II as noted in [14].
The accuracy of the predicted fracture paths in these studies obviously depends on the choice of step size Da and the procedure used to determine path evolution. In the present investigation, we have employed the simple scheme of using a ®xed value throughout the simulation and an explicit approach for path evolution, as in [18]. With this approach, crack growth increments of approximately one tenth of the physical crack length and twice the average element diameter have provided the best comparisons with experiments. An alternative technique for determining crack growth increment and trajectory was presented in [24]. Fig. 17 plots the equivalent mode I stress intensity factor as a function of crack length, normalized by an appropriate measure. The equivalent mode I stress intensity factor is a measure of the mixed-mode stress ®eld, and is given by where h h is the propagation angle, determined from the maximum hoop stress criterion. The ®gure shows K eq I decreases with increasing crack length, implying stable crack growth.  As a last example, we simulate crack growth for a collinear array of four cracks oriented at h p=4 for a total of six steps. The ®nal con®guration and displacement contours are shown in Fig. 18. The cracks grow in the direction of the applied load, and do not intersect one another. These results are consistent with the experimental observations for the same con®guration presented by Nemat-Nasser and Horii [19].

Summary and concluding remarks
The objective of this work was to develop a method to enforce nonlinear constitutive laws on arbitrary interfaces. This was eected through a combination of enriched X-FEM approximations and the iterative strategy adopted from the LATIN method. With a single mesh, several dierent interface geometries were  considered by changing the enrichment scheme. Once the geometry of the interface was set, several dierent constitutive laws were modeled by only changing the local update scheme. It bears emphasis that this technique is applicable to problems outside of fracture mechanics. Future work includes the extension of the method to time dependent problems with evolving interfaces.

E
for plane stress; E 1 À m 2 for plane strain: Expanding and rearranging terms gives where I 1;2 is called the interaction integral for states 1 and 2 where W 1;2 is the interaction strain energy Writing Eq. (A.1) for the combined states gives after rearranging terms Equating (A.4) with (A.7) leads to the following relationship: Making the judicious choice of state 2 as the pure Mode I asymptotic ®elds with K 2 I 1, K 2 II 0 gives the Mode I stress intensity factor for state 1 in terms of the interaction integral The Mode II stress intensity factor can be determined in a similar fashion. The contour integral (A.5) is not in a form best suited for ®nite element calculations. We therefore recast the integral into an equivalent domain form by multiplying the integrand by a suciently smooth weighting function qx which takes a value of unity on an open set containing the crack tip and vanishes on an outer prescribed contour C 0 . Then for each contour C as in Fig. A.1, assuming the crack faces are straight in the region A bounded by the contour C 0 , the interaction integral may be written as where the contour C C C C À C 0 andm is the unit outward normal. In deriving the above, we have used the relations m j Àn j on C and m j n j on C 0 , C and C À , and m 1 0 on C and C À . Now using the divergence theorem on the closed integral over the contour C and passing to the limit as the contour C is shrunk to the crack tip, we arrive at the following equation for the interaction integral in domain form: In deriving the above, we have used the knowledge that the expression r 2 i2 m 2 ; i 1; 2 vanishes on the crack faces as the auxiliary ®elds satisfy traction-free crack faces. We also note the importance of including the integral over C and C À in the case of contact between the crack faces. While the normal tractions which arise on the crack faces during contact do not contribute to the J-integral (the normal forces do no work during crack extension), they do contribute to the interaction integral. The actual normal forces do work in conjunction with the auxiliary displacements.
For the numerical evaluation of the above integral, the domain A is set from the collection of elements about the crack tip. In this paper, we ®rst determine the characteristic length of an element touched by the Fig. A.1. Conventions at crack tip. Domain A is enclosed by C, C , C À , and C 0 . Unit normal m j n j on C , C À , and C 0 and m j Àn j on C. crack tip and designate this quantity as h local . For two-dimensional analysis, this quantity is calculated as the square root of the element area. The domain A is then set of all elements which have a node within a ball of radius r d about the crack tip. Fig. A.2 shows a typical set of elements for the domain A with the domain radius r d taken to be twice the length h local . Fig. A.3 shows the contour plot of the weight function q for these elements. The weight function q is taken to have a value of unity for all nodes within the ball r d , and zero on the outer contour. The function is then easily interpolated within the elements using the nodal shape functions.

Appendix B. Local update equations
In this appendix, we provide the local update equations for the constitutive laws considered on the interface. These are obtained by solving the constitutive law in conjunction with the search equations, restated here as t I n À t A n k 0 w I n À w A n VM P C d ; B:1a t IÀ n À t AÀ n k 0 w IÀ n À w AÀ n VM P C À d : B:1b We begin with the simple case of a perfect bond on C d , and then consider unilateral contact with friction.

B.1. Perfect contact
The local update equations for perfect contact are

B.2. Unilateral contact
The indicator C I is introduced 2C I w AÀ À w A Á n À 1 k t AÀ À t A Á n: B:3 If C I > 0 there is separation and the update is: w IÀ w AÀ À k À1 0 t AÀ ; B:4b whereas in the case of contact C I 6 0, the update is t I Á n Àt IÀ Á n kC I ; B:5 Á n: B:6 Concerning the update of the tangential components, we must distinguish between the case with and without friction. If the model does not include friction, the update for the tangential components is simply If the model includes friction, the G I indicator is introduced G I 1 2 kP T w AÀ È À w A À P T t AÀ À t A É B:7 and the sliping threshold g is computed by g jt I Á nj (or g jt IÀ Á nj, since Eq. (B.5) holds) When kG I k < g, there is a`stick' condition on the crack and the update is P T t I ÀP T t IÀ G I ; B:8a P T w I P T w IÀ w A k À1 0 P T t I À t A : B:8b When kG I k P g, there is a`slip' condition on the crack and the update is P T t I ÀP T t IÀ g G I kG I k ; P T w IÀ P T w AÀ 1 k P T t IÀ À t AÀ :