An eﬀective XFEM with equivalent eigenstrain for stress intensity factors of homogeneous plates

Based on the concept of equivalent eigenstrain, a low-order accurate eXtended Finite Element Method (XFEM) is presented. The aim is the determination of the stress intensity factors of cracked homogeneous specimens. The proposed approach diﬀers from the conspicuous amount of existing contributions on this topic. The ﬁndings of the present paper highlight aspects so far neglected in the literature, such as the mechanical meaning of additional ﬁelds and equations speciﬁc to the XFEM approximation of the displacement ﬁeld. Moreover, based on the plane strain examples simulated in the present study, the proposed XFEM is generally computationally more robust and accurate than existing comparable XFEMs, while keeping a minimal implementation eﬀort.


1
The present study proposes an eXtended Finite Element Method (XFEM) 2 where the crack tips are modeled as Eshelby's elastic singularities. The target is 3 to devise an effective model with minimal computational effort. With respect to The solid V displayed in Fig. 1 is loaded with tractions p i on the region ∂V p of the surface ∂V . The body is subjected to vanishing displacements u i on the external surface portion ∂V d . In the absence of the crack, for an isotropic material with Young's modulus E and Poisson coefficient ν, the continuum problem is formulated as where δ ij is the Kronecker delta, E 1 is the first strain invariant and λ = νE (1+ν)(1−2ν) 118 and µ = E 2(1+ν) are the Lamé constants.
I The region Ω is removed. Thus, it undergoes the strain −e * ij . Let σ * ij be the 125 stress associated with −e * ij by means of Hooke's law. 126 II The surface tractions σ * ij n j are applied on ∂Ω to bring Ω back to its primary 127 shape. The inclusion is rewelded in the matrix. The tractions σ * ij n j have 128 originated a layer of body forces b * i = −σ * ij,j spread over ∂Ω.

129
III an opposite tractions field σ * ij n j is applied at ∂Ω so to recover the initial 130 state.

131
Following Eshelby's seminal work [37] and a later, corrected, version [42], both 132 based on Love's classic treatise [43], the displacement at r due to a point-force 135 where 136 U ij = 1 4πµ 137 Therefore, the displacement in the inclusion in stage III due to σ * ij n j is [37] 138 139 where σ * ij is the stress associated with e * ij by means of Hooke's law 141 Remarkably, Eq. (7) suggests that the compatible strain field ij = 1 2 (u i,j + u j,i ) 142 in the inclusion is, in turn, a function of the eigenstrain e * ij .

143
While the strain in both matrix and inclusion is ij = 1 2 (u i,j + u j,i ), the 144 stress in the matrix writes as 146 The stress in the inclusion is cast as 148 with σ * ij given by Eq. (8). The crack tip at which the strain and the stress are singular is regarded as 151 a point singularity, and the crack tip is replaced by a region Ω surrounding the 152 crack tip as shown in Fig. 2. 153 In Ω, the crack is regarded as a disturbance of the strain field. The aim is 154 to correctly capture this disturbance. The first key assumption relies on the Therefore, the displacement field depends on two distinct vector fields of 168 components v i and a i .

169
By compatibility, for infinitesimal displacements and strains, the total strain

172
where the strain term is bounded, while the strain term is singular, because f has singular derivative. The strain field in the whole solid is 180 ij and e * ij being defined by Eqs. (13) and (14), respectively.

181
However, the stress in the inclusion turns out being in Ω. Furthermore, being u i ,and thus ij , function of e * ij through Eqs. (7) and 184 Eq. (8), the stress inside the inclusion is a function of e * ij in its turn, to be cast

187
Somehow reminiscent of Eshelby's tensor [37], the constitutive tensor S ijhk will 188 be specified in the forthcoming Sec. 4.2. Thus, in the inclusion, 190 Equality (18) plays a crucial role in the subsequent developments.

191
The energy momentum. Eshelby's force on an elastic singularity inside a closed 192 surface Γ is the integral holds also in the present case. This can be shown by computing The focus of this section is on the problem of the equilibrium of the solid V of Fig. 1 governed by  Imposition δW = 0 for any admissible variation of the primal fields implies the Euler-Lagrange conditions (3), it can be drawn that the stress σ 22 across the plane of the crack ahead of the tip and the relative displacement of the crack faces just behind the tip are  To determinate S ijhk (17), a heuristic strategy based on Eqs. (26) is adopted.

228
In particular, it is useful for the subsequent developments to compute the vari-229 ation of ∆u 2 along x 1 Attention is restricted to mode I and to the case where function f in the 232 approximated displacement (11) boils down to

239
At distances x 1 sufficiently small from the crack tip,

241
The variation of Eq. (31) is where Eq. (32) has been replaced.  Table 1: Comparison between the standard formulation and the present formulation in Sec. 5.2 is solved by assuming the following expression of S ijhk 254 α being the inclination of the crack with respect to the loading direction.

255
In the applications described in Sec. 5, parameter γ is set equal to 0 when 256 pure mode I opening is expected, while it is non-vanishing for mixed shearing-257 opening mode. Following the classic approach adopted for linear elastic bodies, the total 261 work writes

263
where the stress is σ ij = C ijhk ( hk + e * hk ). Imposition of the stationarity equations for any virtual variation δu i and δa i leads to the system of equations It can be observed that the mechanical meaning of Eq. (37b) is not evident. 282 The vector field f a is approximated by means of in the present paper, the first-order crack tip functions [1] 287 Typically, the approximated domain is split into three sets, namely the set of the 294 finite elements whose nodes are not enriched, the set of the partially enriched 295 elements, called transition elements, and the set of the elements whose nodes The discrete form of the displacement field (41) is rewritten as

311
The compatible strain field is Matrices B M and B F are associated with the gradient of M and with the tensor product ∇F ⊗ Na, respectively. Moreover,F contains the crack tip enrichment functions but does not coincide with F for dimensional consistency reasons. More details on the structure of the aforementioned matrices are in the appendix. Finally, the discrete form of the stresses is are introduced, where S is given by Eq. (35).
for any virtual δv and δa that vanish at the boundaries, being σ and s given by Eqs. (46). The solving equations are obtained by computing the discrete form of the Euler-Lagrange conditions for W Therefore, the stiffness matrix is where S follows from Eq. (34). Let the approximating space be the same as that introduced in the previous section. For the standard XFEM approach, the discrete form of the virtual work function is for any virtual δv and δa, where σ = C(Bδv +F B M δa + B F δa).

323
The stationarity equations of W (50) are The associated stiffness matrix is  Let Eq. (51b) be rewritten as The latter equation is quite restrictive to be satisfied as it has to be satisfied not The "intruder" in Eq (54)  present study is the first to adopt the latter standpoint in fracture mechanics.

361
The main numerical consequence of the adopted approach is the elimination reasoning. One of the crucial relationships of the present XFEM is Eq.(22c), whose discrete one-dimensional expression reduces to • for the displacement discontinuity Heaviside enrichment H,

378
where t is a unit length to be introduced for dimensional consistency;

379
• for the material discontinuity enrichment sign(x), will be different, because matrix S of the eigenstrain approach appears only in 402 the proposed XFEM.

403
The case of the material discontinuity. As a collateral remark, one can observe 404 that, for the present 1D element placed at the right of the singularity distant 405 x 1 from node 1, Eq. (58) writes  constitutive matrix E appears.

548
In particular, the stiffness matrix of the present method can indeed be ob-  where a x ij and a y ij are the components along x and y of vector a ij associated 635 with node i and enrichment function f j . Matrix F is the 2 × (2 N N k ) matrix