Access the full text.
Sign up today, get DeepDyve free for 14 days.
alexandre-externe.martin@edf.fr IMSIA, UMR This paper presents a robust enrichment strategy to model weak and strong EDF/CNRS/CEA/ENSTA ParisTech discontinuities as well as cracks for industrial applications. First, numerical issues 9219, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 encountered with popular extended ﬁnite element approximation spaces are pointed Palaiseau Cedex, France out. Then, the paper gives indications on how to circumvent those issues. The very Full list of author information is originality of the paper relies on questioning the theoretical approximation spaces with available at the end of the article respect to numerical results and to modify accordingly their design. The relationship between the new design and the previous designs is clearly established, in order to highlight the very small implementation cost of the modiﬁcations exposed here. Hence with minimal additional computational cost, gains in accuracy can be signiﬁcant as shown later in the paper. Keywords: Fracture, Quadratic elements, Condition number, X-FEM, GFEM, SGFEM Introduction Strain localization is usually an issue for conventional ﬁnite element approaches due to numerical issues in the softening regime of the stress–strain relation, when the problem becomes ill-posed. Apart from taking care directly of the localization by an adaptation of the mesh to the discontinuity, diﬀerent methods have been used in the literature to circumvent this diﬃculty. Smeared-cracked models were ﬁrst proposed [1,2]withper- turbations of the ﬁelds across the interface. Discontinuous embedded elements appeared almost at the same time with initial papers of Ortiz et al. [3] or Dvorkin et al. [4], with arbitrary orientations through an element but independent from an element to its neigh- bor. A little bit later, special interface elements [5,6] were proposed localized in between conventional elements, which require frequent re-meshing and reﬁned meshes in order to allow for crack propagation in the correct direction. The eXtented Finite Element Method X-FEM [7] and the Generalized Finite Element Method GFEM [8] ﬁnally allow meshes not to respect the crack geometry while providing a continuous transfer of information from one element to the next one about the crack surface localization unlike the gener- alized class of embedded discontinuous ﬁnite element approaches [3,4]. These methods with nodally based enrichments have managed to combine performance and robustness, © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 0123456789().,–: vol Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 2 of 51 considering non-meshed cracks in a ﬁnite element framework. X-FEM and GFEM use the Partition of Unity [9] and enrich the classical basis of shape functions with discontin- uous functions [10]. The discontinuity of the displacement ﬁeld across the crack surface is then introduced by a generalized Heaviside function, and adding asymptotic ﬁelds at the front crack gives good precision in linear elastic fracture mechanics [10,11]. The main advantage of these methods in comparison to mesh-less methods is their easy implemen- tation in a general ﬁnite element software, and their capabilities to be applied to various ﬁelds: large transformations [12]orplasticity[13] in the X-FEM context for example… One can say that X-FEM and GFEM extend the possibilities of the FEM, while keeping all its advantages. A useful amelioration has been proposed by Sukumar et al. [14]withthe introduction of level set functions to represent discontinuities (cracks, voids,…). These approaches are extremely handy in 3D to treat crack propagation [15,16]. Moreover, these methods have been implemented to solve linear elasticity fracture mechanics problems [8,17] with better accuracy than with ﬁnite element methods. How- ever, the extension of those methods to quadratic elements in three dimensions is a big challenge. Meanwhile, industrial numerical tools use extensively quadratic elements and simulations in three dimensions are almost the norm. In the wake of [18–20], this arti- cle attempts to close in the gap between the theoretical methods and issues related to implementation in industrial software. Minor concerns in one dimension with linear elements turn out to be daunting chal- lenges in 3D or with quadratic elements. Those concerns have been described in [17,21– 23] but have raised little interest since recent publications [24,25]. Firstly, they are related to the ill conditioning of strong or weak discontinuous approximations in the general case of non-conforming interfaces and secondly, to the numerical issues related to “geometrical enrichment” techniques, near the crack-tip. The analysis will be developed more particularly in the case of strongly discontinuous approximations because a direct link with conditioning can be clearly established. For weak discontinuities involved in bi-materials for instance, conditioning issues are still present but they are coupled with the quality of the approximation space to represent continuous solutions with discontinuous gradients [14,26–28]. A speciﬁc section will be dedicated to this analysis. In “Strong discontinuity approximation conditioning” and “Singular enrichment space optimality and conditioning” sections, conditioning issues related to strong discontinuity and singular functions are investigated, as well as, strategies available in the literature to solve those issues. In “Approximation spaces in the literature” and “Numerical behav- ior of strongly discontinuous and singular approximations” sections, those strategies are benchmarked with linear and quadratic elements. Further analysis unfolds that quadratic elements emphasize conditioning and accuracy issues almost unseen with linear elements. The results exposed here are general enough, given the wide range of approximation spaces considered in the paper. Strong discontinuity approximation conditioning In mechanical engineering general purpose software with extended ﬁnite elements, strong discontinuities are positioned rather arbitrarily within the bulk of the mesh. A sensitivity study is clearly needed to check out whether strongly discontinuous enrichment strate- Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 3 of 51 gies are dependent on the position of the interface or not (see Fig. 1). In the literature, the answer is far from obvious. Authors suggested that there is indeed a sensitivity of the solution with respect to the position of the interface; nonetheless, they considered this sensitivity as a numerical side-eﬀect and complex numerical solutions have been elab- orated to deal with this issue [21,24]. Although those techniques may work with linear elements, very little is said about their eﬃciency with quadratic elements. As a matter of fact, the asymptotic behavior of strongly discontinuous approximations is not well under- stood. Conditioning and accuracy of solutions with strong discontinuities deteriorate severely when the interface gets close to the vertices of the mesh. For curved interfaces or unstructured meshes, the interface has a strong probability of getting close to one vertex at least. Dealing with this conﬁguration is a requirement to prevent unexpected results when switching meshes or when moving the interface during industrial numerical simulations. Those numerical issues have been studied in the case of X-FEM formulations [24]. With X-FEM, the enrichment functions and the classical shape functions are almost collinear when the interface cuts the mesh near one of its nodes (Fig. 2). Whenever those conﬁg- urations appear, conditioning deteriorates very quickly. In such cases, there are plenty of strategies to limit the condition number large increase. Those strategies are detailed below. Fig. 1 Positions of the interface investigated in the paper. The condition number is optimal in two cases: ﬁrstly, when the interface cuts through the middle of the elements and secondly, when the interface is on the boundary of the elements. Apart from those cases, the condition number soars. This behavior will be investigated later in the paper Fig. 2 Example of conﬁgurations leading to ill conditioning. The speed of increase of the condition number is related to the distance between the interface and the nodes of the mesh. The condition number soars when either the distance N IP goes to zero or when N IP goes to zero 1 2 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 4 of 51 First of all, we can consider the elimination of degrees of freedom, based on straightfor- ward criteria [18,24]. Those criteria are rough estimates of the “Heaviside information” on both sides of the discontinuity. Basically, if the Heaviside information is heavily unbal- anced, we have: – either: measure ∩ Supp(Φ ) measure ∩ Supp(Φ ) (1) + i − i –or: measure ∩ Supp(Φ ) measure ∩ Supp(Φ ) (2) − i + i where, Supp(Φ ) refers to the support of the shape function and and are the i + − domains deﬁned on Fig. 3. Practically, this means that node numbered i is either on the side or on the side − + and does not “see” the information from the opposite side. From now on, the information from the opposite side will be called complementary information. The criterion used here replaces the measurement of volumes [18,24] by the measure- ment of distances between the nodes of the mesh and the interface on cut edges. When these criteria are associated to a relocation of the interface they are usually named ﬁt-to- Fig. 3 Set of enriched nodes for a strong discontinuity. A node is enriched with the jump functions only if his support is cut by the crack’s interface Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 5 of 51 vertex. Along each edge (see example on Fig. 2), the balance of “Heaviside information” is weighted through the comparison of lengths [N IP] and [N IP] with respect to [N N ]. 1 2 1 2 Level-sets (abbreviated “lsn”) are a handy tool to evaluate those lengths, since they do not require computing the coordinates of the intersection points. From the stand of node N ,wehaveif lsn refers to the normal level set: measure( ∩ Supp(Φ ) ) ∝ length([N IP]) ≈ lsn(N ) (3) + 1 1 1 measure( ∩ Supp(Φ ) ) ∝ length([N IP]) ≈ lsn(N ) (4) − 1 2 2 where N and N are the vertices of the linear element pictured (Fig. 2). IP isagiven 1 2 intersection point on the linear element. lsn(x) represents the level-set function at point x. Then, the level-set values of both nodes are compared in order to enrich or not node N . A similar approach consists in changing directly the value of the level-set of the N 1 2 node to zero. In any case, a threshold is needed to decide whether or not the enrichment needs to be modiﬁed. This threshold based on the ratio of characteristic dimensions is −2 −3 generally chosen between 10 and 10 [18,24]. A simple criterion expresses as: • Node N is eliminated or the level set is moved to N by resetting lsn (N )tozeroif, 1 2 2 lsn(N ) −2 < 10 (5) lsn(N ) + lsn(N ) 1 2 • Node N is eliminated or the level set is moved to N by resetting lsn (N )tozeroif, 2 1 1 lsn(N ) −2 < 10 (6) lsn(N ) + lsn(N ) 1 2 Apart from distance or volume weighting criteria associated or not to ﬁt-to-vertex, a more evolved criterion may be used [24]. Here, the additional idea is about weighting the “balance of Heaviside information” inside the stiﬀness matrix, instead of weighting the geometrical information of physical domains. This “stiﬀness” criterion attempts to correlate the condition number of the “stiﬀness” matrix to the numerical threshold of elimination. On second hand, to control the condition number, one can consider an algebraic pre- conditioner [21], dedicated to orthogonalize shape functions and enrichment functions, with a Cholesky factorization. Let’s assume K is the stiﬀness matrix. K is the term of the matrix associated with shape i,x function Φ and the enriched d.o.f. corresponding to the enrichment function F Φ .AsK i x i is bilinear, symmetric and positive deﬁnite there is an inner product ·, · such as: K = Φ ,F Φ (7) i,x i x i K Béchet et al. [21] pre-conditioner orthogonalizes Φ and F Φ (in the sense of ·, · ), i x i K according to the following procedure: K → K = P KP K u˜ = P f Ku = f → (8) u = P u˜ and, K = Φ ,F Φ = P Φ ,P F Φ =0(9) i,x i x i ˜ C i C x i K Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 6 of 51 If the degrees of freedom for each node are adjacent (ranked contiguously), the shape of P is: ⎡ ⎤ P 00 ··· ··· ··· ⎢ ⎥ ⎢ 0 P 0 ··· 0 ···⎥ ⎢ ⎥ ⎢ ⎥ 00 P ··· ··· ··· ⎢ ⎥ ⎢ ⎥ P = (10) ⎢ . ⎥ ··· ··· ··· ··· ··· ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ··· 0 ··· ··· P ··· ⎣ ⎦ ··· ··· ··· ··· ··· The blocks of P express as follows: −1 S if thenodeisenriched P = (11) I if not where the inverse terms are computed from Cholesky’s factorization of local blocks of K : i T K = S S (12) loc i ⎡ ⎤ Φ , Φ Φ ,H Φ Φ ,F Φ i i K i i K i x i K ⎢ ⎥ ⎢ ⎥ K = ⎢ Φ ,H Φ H Φ ,H Φ H Φ ,F Φ ⎥ (13) i i i i i x i K K K loc ⎣ ⎦ Φ ,F Φ H Φ ,F Φ F Φ ,F Φ i x i i x i x i x i K K K Singular enrichment space optimality and conditioning In linear elasticity, the asymptotic displacement solution at the tip of the crack satisﬁes (in the local basis deﬁned, Fig. 4)[10,11]: j j j j j j u(r, θ) = r k u + k u + k u λ = + j (14) 1 I 2 II 3 III j=0 0 1 r θ 1 θ 2 u = cos (κ − cos θ) e + sin (κ − cos θ) e I 2μ 2π 2 2 0 1 r θ 1 θ 2 u = sin (2 + κ + cos θ) e + cos (2 − κ − cos θ) e (15) II 2μ 2π 2 2 0 2 r θ 3 u = sin e III μ 2π 2 where κ = 3 − 4υ (under the plane strain hypothesis) is Kosolov’s constant and υ is Poisson’s ratio. The ﬁrst eigenvalue λ = 0.5 represents the less regular term of this expansion series. The associated functions (∝ r) are responsible for a signiﬁcant loss of accuracy within regular FEM [17]. It is the focus point of the X-FEM enrichment to recover at least an order of convergence in energy norm of min (1.5, m), where m is the interpolation order of the elements (1.5 is related to the regularity of the next eigenvalue λ ). Forlinearelements,anorder of convergence close to m = 1 is expected in the energy norm. For quadratic elements, an order of convergence close to m = 2 is expected in the energy norm, when only the components related to the 0.5 eigenvalue are activated by the boundary conditions. Laborde et al. [17] showed that a ﬁxed enrichment zone is needed, around the crack- tip, to ensure the optimal accuracy of singular approximations (see Fig. 5). This technique uncouples the size of the enrichment zone and the size of the mesh, so that the enrichment zone does not shrink to zero with mesh reﬁnement. However, this enrichment technique has given rise to many concerns almost from the very beginning: Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 7 of 51 Fig. 4 Plane crack local basis at the crack-tip. Deﬁnition of cylindrical coordinates at the crack-tip for plane crack Fig. 5 Sets of enriched nodes around the crack. If the crack ends in the middle of the domain, additional enrichment functions are needed to model the singularity at the vicinity of the crack-tip. a illustration of the sets of nodes for a linear mesh (for which I = I ), b Illustration of the sets of nodes for a quadratic mesh CT CT,1 (for which I = I )and c Illustration of the sets I and δI for a quadratic mesh CT CT,2 CT,1 CT,1 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 8 of 51 • The ﬁrst concern is the deﬁnition of the enrichment area and the related enrichment strategy in blending elements. In the literature, two main strategies emerge: the “cut- oﬀ” strategy of Nicaise et al. [29] and the geometrical enrichment strategy [17]. Firstly, the “cutoﬀ” strategy adds global singular d.o.f. and an additional function to soften the transition between the enriched zone and the non-enriched zone. Secondly, the geometrical enrichment introduces local singular d.o.f. which values are set to zero outside the enrichment zone. Nicaise and et al. [29] shows that both strategies are relevant and do not alter the convergence rates. • The second concern is the swift increase of the condition number. Latest works of Chevaugeon et al. [30] and Guptaa et al. [25] give strong research directions to improve the condition number. A vectorial enrichment reduces drastically the number of d.o.f. needed to describe the singular solution. Guptaa et al. [25] suggests that vectorial enrichment could be further improved to permanently remove conditioning issues. The idea of [25] is to subtract from the enrichment function its linear interpolation on the set of elements where the singular enrichment is deﬁned. Approximation spaces in the literature In the literature, there are numerous approximation spaces to model strong discontinuities and cracks within the general framework of extended ﬁnite element methods and alike. The earliest work of [7] introduced the standard X-FEM approximation space of Table 1 in 1999. Then in 2000, the work of [8] extrapolated GFEM [22] to crack modeling. Although X-FEM and GFEM both aimed to model cracks, GFEM was designed to deal with a wider scope of problems than X-FEM. Initially, X-FEM and GFEM were assumed as diﬀerent approaches. X-FEM introduced scalar functions to model Irwin modes at the crack-tip. On the opposite, GFEM used directly Irwin modes as vector functions to model the singular behavior at the crack-tip. In “Approximation spaces for a cracked domain” section, this diﬀerence is investigated. We will ﬁnd out that X-FEM and GFEM are closely related. In 2004, Laborde et al. [17] paper questioned the numerical accuracy of the X-FEM approximation, particularly with higher order elements. Laborde et al. [17] suggested that X-FEM with one layer of enriched elements at the crack-tip, is not very accurate. When the enrichment zone enlarges, accuracy improves, but conditioning deteriorates. Therefore, it suggested a new approximation space to take care of conditioning issues and with a better behavior for quadratic elements. As this new approximation space is restricted to 2D analysis, this approximation is not ﬁtted for industrial numerical simulations. Finally, to address higher order modeling, X-FEM design evolved into vectorial enrichment [30], close to GFEM. Very recently, in 2013, even the GFEM approximation space has evolved into a new space called SGFEM [25,31]. SGFEM addresses conditioning issues of GFEM with a new deﬁnition of the vectorial functions. Following [32] it is even deﬁned by the fact that the condition number of the associated stiﬀness matrix has to be of the same order with respect to mesh reﬁnement than the one of FEM. However, SGFEM as proposed in [25,31]has some drawbacks that will be discussed in “Approximation spaces for a cracked domain” section. Speciﬁcally, to model strong discontinuities, many approximations have been released in the literature. The earliest work of [7] introduced a Heaviside function to model the Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 9 of 51 jump d.o.f. on nodes in the vicinity of a strongly discontinuous interface (Fig. 3). The next relevant approximation space was the one of Hansbo et al. [33] that introduces discontinuous polynomials within the partition of unity. Numerous spaces were derived from those previous approximation spaces and far too many to be studied thoroughly in this paper [19,20,34,35]. Hence, we will focus, in “Strong discontinuity representation” section, on the ones of Moës et al. [7], Hansbo et al. [33] and Belytchko et al. [36], which is somehow intermediary between those of [7]and [33]. From the mathematical point of view, we will ﬁnd out that the three spaces are very similar. Strong discontinuity representation In this section, we model only a strong discontinuity, with a media fully split by an interface. Hence, we do not need the singular functions in that case. Singular functions will be investigated in the next section. Thus, let us assume the following X-FEM approximation space V , modeling a strong discontinuity, ⎧ ⎫ ⎨ ⎬ h h V = u = a Φ + b H Φ (16) i i j j ⎩ ⎭ i∈I j∈I where Φ are lagrangian polynomials up to order two, in the scope of this article. I denotes the nodes of the ﬁnite element mesh, I the set of enriched nodes and H the Heaviside function. A node belongs to I if, and only if, its support is split by the interface .We recall that H is deﬁned as follows: −1, if lsn(x) < 0 ∀x ∈ ,x ∈ / ,H(x) = (17) +1, if lsn(x) > 0 Let’s consider the basic case where the whole domain is divided into two parts as in Fig. 3: • = (H=+1) + • = . (H=−1) − For which, • u represents the approximation of the solution on the dimain extended by continuity to the domain over the interface • u represents the approximation of the solution on the dimain extended by continuity to the domain over the interface The representation of the kinematic jump is crucial for many interfacial laws, such as cohesive laws and contact laws. The jump at every point (commonly integration point) located on the interface, follows the convention: h h h ∀x ∈ , u (x) = u − u (x) (18) + − In order to simplify the notations in Table 1, we deﬁne for a given node j ∈ I , = H j,1 + − χ and = χ , where χ denotes the characteristic function of the set A.Using j j,2 j A these notations, we have: +2χ =+2 , if j ∈ + j j,1 − H − H x = j j j −2χ − =−2 , if j ∈ j j,2 + Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 10 of 51 Table 1 Strong discontinuity enrichments comparison X-FEM [7] Domain enrich. [33] Shifted enrich. [26] Linear shape functions on enriched nodes h h h u = a Φ + α Φ u = c Φ − 2d Φ i i j,1 j,1 i i j j,2 u = a Φ + b Φ + + j j j j i∈I/I j∈I i∈I j∈I ∩− H H H j∈I j∈I Displacement approximation u = a Φ − b Φ h h j j j j − u = a Φ + α Φ u = c Φ + 2d Φ i i j,2 j,2 i i j j,1 − − j∈I j∈I i∈I/I j∈I i∈I j∈I ∩+ H H H h h h Jump approximation u = 2b Φ u = α Φ − α Φ u =− 2d Φ − 2d Φ j j j,1 j,1 j,2 j,2 j j,2 j j,1 j∈I j∈I j∈I ∩− j∈I ∩+ H H H H X-FEM [7] approximation space is on the left, Hansbo et al. [33] is in the middle, and Fries [26] approximation space is on the right Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 11 of 51 where x is the location of j. The function H − H x is used to deﬁne the formulation j j introduced in [26]. Let us compare side-by-side, the XFEM representation of the displacement jump intro- duced by [7] with the other well-known formulations of [33]and [26]. Even though the jump approximations have diﬀerent expressions (Table 1), we will investigate next whether the approximation spaces are really diﬀerent or not. Areias et al. [37] showed that [7]and [33] involve diﬀerent bases but that both enrich- ments still represent the same approximation space. Each enrichment can be expressed in terms of the other one with a suitable change of variables. With the same analysis as in [37], the diﬀerent enrichments can be expressed in terms of the other ones for the remaining couple of formulations X-FEM/Shifted and Domain/Shifted. Thus, from the mathematical point of view, the approximation spaces involved in the three formulations are equivalent (see Table 2). The same can be said of other formu- lations encountered in the literature [19,34,35]. We would like to stress the fact that this equivalency only holds in the case of strong discontinuity. When singular functions are injected in the approximation, the relationship between the diﬀerent approximation spaces may be diﬀerent as studied in the next section. Approximation spaces for a cracked domain The model problem is a cracked domain (Fig. 6), under a linear elasticity assumption. The material is also assumed to be homogeneous and isotropic. Dirichlet boundary con- ditions are applied on the boundary and Neumann boundary conditions are applied on The space of admissible displacements is: V = v ∈ H (); v = 0 on (19) Table 2 Change of variables when switching enrichment strategies. ∀j ∈ I , the following changes of variables hold Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 12 of 51 Fig. 6 Example of a cracked domain. Deﬁnition of notations used to label the domain and its boundaries The weak form of the equilibrium problem is: Find u ∈ V such that a(u, v) = l(v) ∀v ∈ V, a(u, v) = σ u : ε v d, ( ) ( ) P: (20) ⎪ l(v) = f · vd + g · vd σ (u) = λtr(ε (u))1 + 2με (u) , where u is the displacement, σ is the Cauchy stress, ε is the strain, 1 the identity tensor, f is the body force applied on ω, g is the traction applied on and λ and μ are the Lamé parameters. Now, the discrete approximation spaces from the literature will be investigated. In the fol- lowing sections, we will prove that spaces available in the literature can be gathered into two classes: “straightforward” enrichment and “bubble” enrichment. An illustration of those two classes is provided on Fig. 7. “Straightforward” enrichment means that the sin- gular function is introduced directly into the approximation space without modiﬁcation. “Bubble” enrichment means that the singular function is reshaped into a new function before its introduction into the approximation space. In the following sections, those classes are described: • “Straightforward” singular enrichment class” section shows that X-FEM and GFEM are “straightforward” enrichments. Moreover, GFEM is a subspace of X-FEM. This statement implies that X-FEM is somehow more accurate than GFEM, but both methods are nonetheless very close. • “Bubble” enrichment class” section shows that a “bubble” enrichment as the one of Gupta et al. [25] is not a “straightforward” enrichment. This statement implies that “bubble” enrichments have numerical properties not identical to “straightforward” enrichments that will be investigated in “Singular approximation at a crack tip” sec- tion. Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 13 of 51 Fig. 7 “Straightforward” enrichment class versus “bubble” enrichment class. The diﬀerence in shape highlights the diﬀerence in behavior stated in the previous section “Straightforward” singular enrichment class Let V be the approximation space using vector asymptotic functions to describe the crack kinematics: d d v = a Φ E + b H Φ E i,j i i,j i h i∈I j=1 j i∈I j=1 j V = (21) + c Φ K i,α i i∈I α=1 CT where d is equal to two in the bidimensional case and three in the tridimensional case, E ,E ,E is the usual Cartesian basis, K are the vector asymptotic functions and I CT 1 2 3 α the set of nodes enriched by such functions. In this contribution, I is deﬁned as the set CT of nodes placed at a distance lower than R from the crack-tip (cf. Fig. 5). The vector asymptotic functions K are deﬁned as follows: K = u ⎨ 1 (22) K = u II K = u III 0 0 0 with u ,u u given by (15). I II III h 2,∞ Let V (W ) be the approximation space deﬁned with a cutoﬀ function χ ∈ CUTOFF 2,∞ W [29], which value is set to one if its distance from the crack front is lower than r , is set to zero if its distance from the crack front is greater than r with r > r and varies 1 1 0 linearly in between: Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 14 of 51 d d v = a Φ E + b H Φ E + i,j i i,j i j j h 2,∞ i∈I j=1 i∈I j=1 V (W ) = (23) CUTOFF d c χK α=1 α and let V be the standard X-FEM approximation space [7]: XFEM ⎧ ⎫ d d ⎨ ⎬ v = a Φ E + b H Φ E + i,j i i,j i i∈I j=1 j i∈I j=1 j V = (24) XFEM ⎩ α ⎭ c F Φ E k∈I α=1,4 j=1 j k,α CT where F are deﬁned as: F = r sin ⎪ 2 ⎪ √ ⎨ θ F = r cos √ (25) F = r sin sin θ F = r cos sin θ We remark that: 2 3 r cos cos θ = F − F √ (26) 4 1 r sin cos θ = F − F Consequently, the vector functions K and the scalar functions F are related by the following linear relations: 1 1 2 3 1 1 4 2 √ √ K = (κ − 1) F + F e + (κ + 1) F − F e 2μ 2π 2μ 2π 1 1 4 1 1 2 3 2 √ √ K = (κ + 1) F + F e + (1 − κ) F + F e (27) 2μ 2π 2μ 2π ⎩ 2 1 3 K = F e μ 2π 2,∞ h h h h Lemma 1 V (W ) is a subspace of V and V is a subspace of V , i.e. the CUTOFF h XFEM three approximation spaces are related as follows, 2,∞ h h h V (W ) ⊂ V ⊂ V CUTOFF h XFEM 2,∞ where W is the subspace of interpolated cutoﬀ functions. Proof The proof relies on the work of [29]. 2,∞ h h • V (W ) ⊂ V : CUTOFF 2,∞ h h For v ∈ V (W ), CUTOFF d d d v = a Φ E + b H Φ E + c χ K i,j i i,j i α h j j α i∈I j=1 i∈I j=1 α=1 χ is the interpolation of a cutoﬀ function χ which vanishes outside I , h CT χ = χ (x ) Φ i i i∈I CT Then, d d d v = a Φ E + b H Φ E + c χ (x ) Φ K i,j i i,j i α i i j j α i∈I j=1 i∈I j=1 i∈I α=1 H CT which leads to, h h v ∈ V Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 15 of 51 h h • V ⊂ V : XFEM h h For v ∈ V , d d d v = a Φ E + b H Φ E + c Φ K i,j i i,j i k,α k j j α i∈I j=1 i∈I j=1 k∈I α=1 H CT 4 d l m • Equation (27) shows that K can be expressed as K = k F e , where l,m,α α α l=1 m=1 k is a constant tensor. i,m,α In case of a straight crack, the local basis is constant. Thus, the local basis expresses in cartesian coordinates as, e = μ E , where m,j j=1 μ is a constant tensor. m,j Then, 4 d d K = k μ F E l,m,α m,j α j m=1 l=1 j=1 d d 4 d d c Φ K = c k μ F Φ E k,α k k,α l,m,α m,j k α j α=1 α=1 m=1 l=1 j=1 ⎛ ⎞ 4 d d d ⎝ ⎠ = c k μ F Φ E k,α l,m,α m,j k l=1 j=1 α=1 m=1 With a suitable change of variable, d d c = c k μ k,α l,m,α m,j k,l α=1 m=1 we have ﬁnally, d d 4 d h l v = a Φ E + b H Φ E + c F Φ E i,j i i,j i j j k j k,l i∈I j=1 i∈I j=1 k∈I l=1 j=1 CT which means that: h h v ∈ V XFEM • In case of a curved crack, the result of the lemma holds with a discretization of the 1 2 3 local basis of [30]. Let us consider a discrete local basis at each node e ,e ,e , such k k k as, d d d v = a Φ E + b H Φ E + c Φ K i,j i i,j i k,α k j j α,k i∈I j=1 i∈I j=1 α=1 H k∈I CT 4 d l m K = k F e l,m,α α,k k m=1 l=1 e = μ E m,j,k j=1 The straight crack proof still holds here, with the additional node index “k”. Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 16 of 51 Remark as a practical consequence of Lemma-1, the error bound of X-FEM should be lower than the ones of V and of the cutoﬀ [29]. With X-FEM, the optimization process performs on a larger space than with other formulations and reaches a closer inﬁmum to the exact solution of the problem. Nevertheless, X-FEM introduces more d.o.f. than GFEM and cutoﬀ enrichments which increases its condition number. “Bubble” enrichment class Now, let us consider the SGFEM space for linear elastic fracture mechanics similar to the one proposed in [25,31]: d d 4 v = a Φ E + b H ψ Φ E i,j i i i,j,k i,k h i∈I j=1 j i∈I j=1 k=1 j V = (28) SGFEM d + c Φ K − K i,α i i∈I α=1 α α CT where ψ = 1,x − x ,y − y ,z − z ,and K = K (r , θ )Φ is i i i i,k α { } α k k k k∈ I ∪δI k=1..4 CT CT the usual interpolation. δI is deﬁned here as the transition layer of nodes between crack-tip elements and other CT elements of the mesh for which K unknowns are set to zero. This allows to connect the enriched layer with the remaining of the mesh and to avoid blending issues [26]. h h h h Lemma 2 V is not a subspace of V and V is not a subspace of V ,so XFEM SGFEM SGFEM XFEM that: h h h h V ⊂ V and V ⊂ V (29) XFEM SGFEM SGFEM XFEM h h Proof – V ⊂ V SGFEM XFEM h h h It is obvious that V ⊂ V because V involves higher order discontinuous SGFEM XFEM SGFEM polynomials H ψ Φ , which are out of V . i,k XFEM h h – V ⊂ V XFEM SGFEM 1 h • F E Φ belongs to the test function space V with a = 0,b = k i,j i,j 1 XFEM k∈I CT 0,c = 1and c = 0 elsewhere (for α = 1). k,1 k,α 1 h • However, F E Φ does not belong to V , which can be shown by SGFEM k∈I CT contradiction. 1 h • Let us assume F E Φ belongs to V , SGFEM k∈I CT d d 4 F E Φ = a Φ E + b H ψ Φ E k i,j i i,j,k i,k i 1 j j k∈I i∈I j=1 i∈I j=1 k=1 CT H + c Φ K − K i,α i α α α=1 i∈I CT • Singular functions cannot be represented by discontinuous polynomials basis, so we necessarily have b = 0 i,j,k • Singular functions cannot be represented by continuous polynomials as well, d d so, a Φ E − c Φ K = 0 as Φ K is also a polynomial. i,j i i,α i i j α α i∈I j=1 i∈I α=1 CT • Φ K introduces at least three higher order polynomials per node, α Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 17 of 51 • Which leads to, c Φ K = 0, i,α i α=1 • Then, c = 0, i,α • And, a = 0 because a Φ E = 0. i,j i,j i j=1 • Thus, F E Φ = 0 which is contradictory. k∈I CT Remark As a practical consequence of Lemma-2 the “bubble” class of “Bubble” enrich- ment class” section does not belong to the “straightforward” class of “Straightforward” singular enrichment class” section, in which X-FEM is found to be the largest space. Hence, the deﬁnition of two classes of enrichment makes perfectly sense. Numerical behavior of strongly discontinuous and singular approximations In this section the numerical behavior of the diﬀerent approximation spaces introduced previously is investigated through a couple of benchmarks. First of all, we assess the asymp- totic behavior of strongly discontinuous approximations, through a one dimensional case. Secondly, we consider the convergence study of a two-dimensional problem with a crack and singular approximation. Numerical behavior of strongly discontinuous approximations Numerical analysis of strongly discontinuous approximations with linear elements Let us consider the basic case of the traction of a one-dimensional bar with a ﬁctitious bi- material interface. The arbitrary interface is positioned along the abscissa x (ε) (Fig. 8). The problem even though continuous is treated as a discontinuous one, with gluing interface boundary conditions imposed through a Lagrange multiplier. Actually, if the materials were diﬀerent on each side of the interface, the resulting problem would be equivalent to enforcing a weak discontinuity with a strong discontinuity framework. This artefact is used so as to establish convergence results in energy norm. If it was not the case, the solution obtained would be that of two rigid bodies with prescribed displacements on one of their end. Fig. 8 1D case with X-FEM interface and linear elements. Left, the three elements used in the model. Right, the analytical solution expected. The use of a discontinuous function is possible because discontinuous degrees of freedom bear Neumann conditions that enforce the continuity of the solution (as in a contact problem with Lagrange multipliers at the interface) Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 18 of 51 • The problem is the solution of the following diﬀerential equation: u (x) = 0 ∀x ∈ = [0, 3] (30) • With arbitrary Dirichlet boundary conditions: u(0) = β u(3) = 3α + β (31) • And with continuity conditions at the interface expressed in terms of the continuity of the derivative of the displacement ﬁeld (stress continuity at the interface): du du − + − + = α = α α = α = α (32) dx + dx − x ε x ε ( ) ( ) This continuity condition results from the equivalent Lagrangian form of (30) in which a Lagrange multiplier λ is used to impose a continuous displacement across the interface. (u, λ) is the saddle point of the following functional: + − L (u, λ) = u dx + λ u x ε − u x ε = 0 The expected continuous solution satisfying the boundary conditions above, is: u(x) = αx + β , ∀x ∈ [0, 3] (33) The resulting weak form of the problem is: u ∈ : w ∈ H /w(0) = β,w(3) = 3α + β , λ ∈ Find, ∀v ∈ : w ∈ H /w(0) = 0,w(3) = 0 ∀μ ∈ R u v dx = 0 + + − − → λ − α v x ε − λ − α v x ε = 0 (34) + − μ u x ε − u x ε = 0 The X-FEM discrete linear space is: w/w = a + a + b H + a + b H + a 1 1 2 2 2 2 3 3 3 3 4 4 : (35) a = β a = 3α + β 1 4 Then, the 6 × 6 matrix associated with the discretization of the weak form on the X-FEM space of linear functions (ﬁrst column of Table 1)is: ⎧ ⎫ ⎡ ⎤ 1 −110 0 0 ⎪ a ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ −12 −2 + 2ε −11 − 2ε 0 a ⎢ ⎥ ⎪ 2 ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎨ ⎬ ⎢ ⎥ 1 −2 + 2ε 21 − 2ε −10 b ⎢ ⎥ (36) ⎢ ⎥ ⎪ ⎪ 0 −11 − 2ε 22ε −1 a ⎢ ⎥ ⎪ 3 ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎣ 01 − 2ε −12ε 2 −1⎦ ⎪ b ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ 00 0 −1 −11 a With the domain formulation of the second column of Table 1, the discrete matrix is: ⎧ ⎫ ⎡ ⎤ 1 −10 0 0 0 ⎪ α ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ −12 − ε 0 −1 + ε 00 α ⎢ ⎥ ⎪ ⎪ 2,1 ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎨ ⎬ ⎢ ⎥ 00 ε 0 −ε 0 α 2,2 ⎢ ⎥ (37) ⎢ ⎥ ⎪ ⎪ 0 −1 + ε 01 − ε 00 α ⎢ ⎥ ⎪ 3,1 ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎣ 00 −ε 01 + ε −1⎦ ⎪ α ⎪ 3,2 ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ 00 0 0 −11 α 4 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 19 of 51 With the shifted formulation of the third column of Table 1, the discrete matrix is: ⎡ ⎤ ⎧ ⎫ 1 −10 0 0 0 ⎪ c ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎢ −12 2ε −12 − 2ε 0 ⎥ ⎪ c ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎨ ⎬ ⎢ ⎥ 02ε 4ε −2ε 00 d ⎢ ⎥ (38) ⎢ ⎥ ⎪ ⎪ 0 −1 −2ε 2 −2 + 2ε −1 c ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎣ 02 − 2ε 0 −2 + 2ε 4 − 4ε 0 ⎦ d ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ 00 0 −10 1 c We enforce boundary conditions directly on the discrete problem, as equality constraints through Lagrange multipliers. The position of the interface depends only on ε. Nonethe- less, the solution does not depend on the ε parameter, so that the relative error does not depend on ε, which makes our analysis relevant. Let us consider the following arbitrary numerical values: α = 10/9and β = 10. Then, the solution is computed in standard 64-bit arithmetic, with the linear solver UMFPACK, to reproduce the behavior of direct solvers. The numerical error has to be close to zero, for any given position of the interface. The error is evaluated here in terms of the H -norm: ( ) u − u 1 = (u − u ) + u − u dx (39) h H h Although all three formulations represent the same approximation space, they do not have exactly the same numerical behavior (Fig. 9), at least concerning the error in H -norm. The X-FEM error increases steadily while the “shifted” formulation and the “domain” −16 formulation levels of error stay close to machine precision (about 2.2 · 10 ). Moreover, all three enrichment conditionings increase as the distance between the interface and node N decreases (Fig. 10). This means that all three enrichment strategies are very sensitive to the position of the interface. A pre-conditioner is needed for the three enrichments. In the case of X-FEM, we studied Béchet et al. [21] pre-conditioner. For other formulations we used a diagonal pre-conditioner to scale the d.o.f. The diagonal pre-conditioner allows a scaling of rows and columns through a multipli- cation with a diagonal matrix to the left and to the right: ⎪ K = D KD c c Ku = f → K u = f with u = D u (40) f = D f In linear elasticity, D is usually deﬁned as, 1 max(K ) + min(K ) i,i i,i [D ] = * (41) i,i i,i Thus, on Fig. 11, condition numbers are reduced drastically. However, the error still increases in case of X-FEM. The error is clearly not only related to conditioning. Another explanation should be considered. When looking at the X-FEM discrete stiﬀness matrix, it is noticeable that the unknowns a and b obey to an almost similar equation (as H ≈− ). The diﬀerence between 2 2 2 2 those equations lies in the terms 1 − 2ε and 2 − 2ε. Those terms imply the sum of heterogeneous quantities that leads to a tremendous loss of accuracy on the diﬀerence information 2ε. For instance, let us consider the following string of calculations in “double precision” arithmetic (8-octet storage): 1 − (1 − 6.3847497084E − 12) ≈ 6.348E − 12 (42) Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 20 of 51 Fig. 9 H1-error of strongly discontinuous approximations with linear elements. Shifted and domain enrichments outperform X-FEM formulation when “eps” (ε)goestozero Fig. 10 Conditioning of strongly discontinuous approximations without the use of pre-conditioners. The three approximation spaces yield the same condition number (as they are equivalent) Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 21 of 51 Fig. 11 The use of dedicated pre-conditioners improves the condition number but has little eﬀect on X-FEM accuracy. According to the curves above, conditioning and accuracy are clearly two separated issues: improving the condition number does not automatically solve accuracy issues. Conditioning and accuracy are both symptoms of the faulty behavior of partition of unity with asymptotically small domains The ﬁnal result 6.348E–12, is quite diﬀerent from the exact result 6.3487497084E–12: only four digits remain of the “diﬀerence information” between the equations. Hence ill conditioning indicates also a lesser accuracy within the assembled stiﬀness matrix due to round-oﬀ errors and truncated information. Even highly eﬀective pre-conditioners [21] cannot recover the truncated information. It is not surprising that, besides precondition- ing, authors have used triple precision arithmetic to recover accuracy [24]. Remark Because the error curves are shattered, the phenomenon might be also proba- bilistic and related to random round-oﬀ errors in interaction with the solver algorithm. A comprehensive study with a wide range of solvers is presented (Figs. 12 and 13). This study shows that the type of solver interferes with the level of error, which makes numerical analysis on errors quite sensitive: • For iterative methods (PCG), the error is generally higher because those solvers are very sensitive to conditioning, • For matrix factorization based solvers (UMF), the error is lower than with iterative methods, • For matrix inversion based solvers, good results can be obtained frequently. Numerical analysis of strongly discontinuous approximations with quadratic elements The same one-dimensional case is considered here. The same discretization of the weak form is used with quadratic shape functions. The results are plotted on Fig. 14. For the three approximations, conditioning worsens very quickly, at about three times the rate of linear elements. More worrying is that the error increases. X-FEM enrichment error worsens more quickly than the others. Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 22 of 51 Fig. 12 Overall error analysis of a wide range of solvers for the 1D case. Only X-FEM’s enrichment without pre-conditioner is studied here. The whole trend suggests a global increase of the sensitivity of solvers (with respect to epsilon), when the “epsilon” parameter goes to zero Fig. 13 The use of dedicated pre-conditioner stabilizes the behavior of solvers. But the error still increases with the X-FEM method, as explained in “Numerical analysis of strongly discontinuous approximations with linear elements” section On Fig. 15, the same pre-conditioners than in “Numerical behavior of strongly discon- tinuous approximations” section are considered. The wrong behavior noticed above still applies. This numerical behavior with quadratic elements is well explained by the coupling between enrichment functions, as pictured on Fig. 16 with domain enrichment [33]. When measure( ∩ Supp(Φ ) ) → 0 for a given vertex node, the middle node satis- + S ﬁes likewise measure( ∩ Supp(Φ ) ) → 0, on the vertex patch. For such a conﬁgu- + M ration, the enrichment shape functions χ Φ and χ Φ become almost homothetic + S + M (Fig. 16). Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 23 of 51 Fig. 14 H1-error and conditioning analysis with quadratic elements (Lagrange). The conditioning slopes are about three as shown is “Numerical analysis of strongly discontinuous approximations with quadratic elements” section Fig. 15 The use of linear element pre-conditioners barely improves the accuracy and condition number of quadratic elements (Lagrange) • In case of Lagrange polynomials, we have: Φ (ξ) = (1 − ξ)(1 − 2ξ), ∀ξ ∈ 0, 1 [ ] Φ (ξ) = 4ξ(1 − ξ), ∀ξ ∈ 0, 1 (43) [ ] Both functions can be related with the following expression, Φ (ξ) =−4Φ (ξ) + 4(1 − ξ) (44) S M Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 24 of 51 Fig. 16 Coupling between Lagrange quadratic shape functions. When the crack passes close to the nodes, the discontinuous shape functions of the middle node and one vertex node become almost colinear interface e (x) = (1 − x)χ 1 [] 0,1 (a) 1 2 3 4 5 6 7 interface e (x) = (2 − x) χ 2 [] 2−ε ,2 (b) 1 2 3 4 5 6 7 Fig. 17 Graphs of polynomial functions e (a) and e (b) 1 2 when ξ goes to 1, the diﬀerence information between those functions decreases accord- inglyto(1 − ξ) . This generates a vanishing sub-space (consisting of polynomials only deﬁned on the small support ∩ Supp(Φ ) ), which implies a conditioning slope of + S at least 3 with quadratic shape functions (Fig. 15). In other words, the condition number κ(K)isbounded by: −3 κ(K ) ≥ Cε (45) Proof Let κ(K ) be the condition number of the stiﬀness matrix: max e Ke T e Ke e =1 2 1 κ(K ) = ≥ ∀e ,e / e = 1, e = 1 (46) 1 2 1 2 2 2 T e Ke min e Ke 2 e =1 where is the Euclidian norm 2. In order to estimate the lower bound of the condition number, we consider the polyno- mials e and e illustrated on Fig. 17. 1 2 Those polynomials belong to the discrete space of [33] (likewise [7]and [26]). They can be used to determine a lower bound of the condition number. Functions e and e are 1 2 expressed in the discrete space of [33]as: Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 25 of 51 e =− 2 + ( ) 1 1 2 (47) e = 4χ − + χ − ( ) 2 3 4 In the discrete space of [33], we have: 2 2 2 u = a + α + α (48) i j,1 j,2 i∈I /I j∈I j∈I H H H so that we have e = e = 1. 1 2 2 2 1 T 1 T As e Ke and e Ke represent the elastic energy of these solutions, we have: 1 2 2 1 2 2 1 2 e Ke e dx 1 0 1 = = T 2 2 3 80ε e Ke 2 e dx 2−ε 2 Thus, −3 κ(K ) ≥ Cε The asymptotic behavior of Φ (ξ) [at the ﬁrst order of (1 − ξ)] is: Φ (ξ) ≈−Φ (ξ)/4 M S Then, the related d.o.f. become redundant, which leads to ill conditioning. As enrichment strategies are equivalent, the redundant d.o.f. pollutes also the approximation spaces of the other formulations [7,26]. • In case of Bernstein polynomials, we have: Φ (ξ) = (1 − ξ) , ∀ξ ∈ [0, 1] Φ (ξ) = 2ξ(1 − ξ), ∀ξ ∈ [0, 1] (49) when ξ goes to 1, there is no coupling between Bernstein polynomials as observed with Lagrange polynomials. However, the polynomial Φ cancels out accordingly to (1 − ξ) , which also leads to the same conditioning slope of three (Fig. 18) (as shown in the previous section). It is noticeable that [33] has better accuracy with Bernstein polynomials than with Lagrange polynomials (Figs. 18, 19). Numerical behavior of singular approximations at the crack tip Numerical analysis of crack approximations with linear elements: crack opening in mode 1 Given the capabilities of the software we used, only two approximations were tested here: • X-FEM/GFEM vectorial enrichment [30], • X-FEM scalar enrichment [7]. Other approximations are out of scope because, – Gupta et al. SGFEM kind of enrichment [25,31] (let us recall that accordingly to [32] a genuine SGFEM enrichment is deﬁned by the fact that the condition number of the associated stiﬀness matrix is of the same order—2 [38]—with respect to mesh reﬁnement than the one of FEM) needs far too many additional Heaviside d.o.f. which are diﬃcult to implement and may lead to conditioning issues, Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 26 of 51 Fig. 18 Asymptotic behavior with Bernstein polynomials. Bernstein polynomials allow emphasizing once more that conditioning and accuracy are separated issues: the condition numbers of all three approximation spaces are almost the same, but there is a great diﬀerence in accuracy when “eps” (ε)goestozero – Cut-oﬀ enrichment is not very convenient, as assembling global d.o.f. is out of the scope of the industrial software we used (Code_Aster). Moreover, the method does not extend properly in 3D [29]. Example 1 (Horizontal crack opening in mode I with linear elements). The domain geom- etry is a square deﬁned by = − 0.5, + 0.5 × − 0.5, + 0.5 . The meshing procedure [ ] [ ] subdivides the domain into regular sized cells, which edges are parallel to the crack. The size of the enrichment zone is r = 0.1. The analytical solution in displacement corresponds to the exact mode I limited to the r term in Williams series expansion [11] and is given by the ﬁrst term of (15). Then, Dirichlet boundary conditions are applied on three sides. On the left side, Neumann conditions are preferred. As the left side is cut by the crack, Dirichlet conditions are more diﬃcult to apply given the discontinuity of the displacement ﬁeld (Fig. 20). The condition number is estimated by MUMPS solver, so that the results shown here have to be taken only as rough estimates. On Fig. 21, the condition number of X-FEM almost skyrockets as noticed in [21]. We did not consider the pre-conditioner of [21] here to show how both enrichment strategies, scalar and vectorial, behave without expensive treatment. Even the conditioning slope of the vectorial enrichment is far from optimal. The expected optimal conditioning slope is two [25]. Then, the relative error in energy norm is computed accordingly to the following for- mula: $ $ h h u − u = σ − σ : ε − ε d σ : εd (50) energy Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 27 of 51 Fig. 19 Asymptotic behavior with Bernstein polynomials with dedicated pre-conditioners (same pre-conditioners as with linear elements). As with Lagrange polynomials, dedicated pre-conditioners do not solve accuracy issues Fig. 20 Boundary conditions around the patch-test. Mixed boundary conditions are applied to avoid enforcing displacement on the nodes of the boundary near the crack interface which would add only a technical diﬃculty beyond the scope of the paper On Fig. 22, the rate of convergence in energy norm of the X-FEM scalar enrichment is under-optimal (around 0.91) as noticed in [17]. Optimality is recovered with the vectorial enrichment, the convergence rate being then 0.989. Nonetheless, X-FEM is more accurate than the vectorial enrichment as predicted in Lemma 1. Example 2 (Inclined crack opening in mode I with linear elements). We introduce a more realistic case with a chosen inclination of the crack at 44.9 , in order to test the conditioning with both Heaviside and singular enrichments. The analytical solution is still the one of a mode I problem limited to the r term in Williams series expansion given by the ﬁrst term of (15)(Fig. 23). The inclined crack analytical solution is expressed through a rotation of 44.9 of the horizontal crack solution. The strain and stress tensors are rotated as well of the same angle. Then the same mix of Neumann and Diriclet boundary conditions Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 28 of 51 Fig. 21 Conditioning of crack approximations. The condition number increases swiftly with the geometrical enrichment strategy i.e. the introduction of singular functions on many layers of elements around the crack-tip Fig. 22 Convergence analysis of crack approximation spaces. The single-layer enrichment (topological) convergence curves are added to recall the beneﬁt of the geometrical enrichment strategy are applied (Fig. 24), similarly to the case of the horizontal crack. We consider the same regular meshes as above which are not oriented accordingly to the crack surface here. The pre-conditioner of [21] is not applied here as in the previous case. On Fig. 25, the conditioning of the X-FEM vectorial enrichment becomes very high and reaches a steady value around 10 . This behavior may be explained by the ill con- ditioning of the Heaviside enrichment. As the “Heaviside information” becomes very unbalanced, its ill conditioning overcomes the regular increase of conditioning expected with the vectorial enrichment (see Fig. 21) as explained in “Strong discontinuity approxi- mation conditioning” section. Nonetheless, both enrichments have optimal convergence accuracy (see Fig. 26). Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 29 of 51 Fig. 23 Analytical solution corresponding to a mode 1 displacement ﬁeld for the inclined crack Fig. 24 Boundary conditions for an inclined crack. The boundary conditions are the same with the horizontal crack test Fig. 25 Conditioning of vectorial and scalar enrichments with a crack orientation at 44.9 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 30 of 51 Fig. 26 Convergence analysis for vectorial and scalar enrichments with a crack orientation at 44.9 . Both enrichments have optimal convergence order but error levels are quite high. Results are in accordance with the ones of Fig. 22 and X-FEM is more accurate than GFEM Numerical analysis of crack approximation with quadratic elements: crack opening in mode 1 Only the linear vectorial enrichment approximation [30] is working with a reasonable condition number. We recall that other approximations are not tested because, • X-FEM conditioning skyrockets with quadratic elements, even with a small enrich- ment zone, • Gupta et al. SGFEM kind of enrichment [25,31] needs far too many additional Heav- iside d.o.f., so that it is diﬃcult to implement and may lead to ill conditioning, • Cut-oﬀ with global d.o.f. assembling requires a “macro” element, which is not straight- forward to implement in ﬁnite element software. Moreover the method does not extend properly in 3D [29]. Even if the linear vectorial enrichment converges at an optimal rate (in energy norm) with quadratic elements, the condition number is still very high (see Fig. 27). We did not use the pre-conditioner of [21], in order to test the conditioning of the vectorial enrichment. In “Singular approximation at a crack tip” section, a new enrichment strategy with better intrinsic numerical behavior, will be exposed. Remark a dedicated integration scheme is needed at the tip of the crack, to get an optimal rate of convergence in energy norm with quadratic elements [17]. Here, we used a Gauss- Radau integration rule of order 20 [39]. As the aim of the paper is not about singular integration, we did not try to optimize the procedure. Optimization of integration schemes has been well studied in [30,40]. Improving the design of quadratic approximation spaces to deal with previous numerical issues In this section, basic improvements are suggested to deal with numerical issues stated in previous sections: – For strongly discontinuous approximations, we suggest a correction to the approx- imation spaces for quadratic elements. Note that as denoted by [17] optimal con- Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 31 of 51 Fig. 27 Numerical behavior of the X-FEM vectorial enrichment with quadratic elements vergence rates can only be achieved if the discontinuous approximation space is of the same order than the one of the continuous space, that is to say quadratic. The correction exposed here on the quadratic form is slightly diﬀerent from the ones dis- cussed in “Strong discontinuity approximation conditioning” section (ﬁt-to-vertex and elimination of d.o.f.). – For singular enrichment, a reshape of the approximation space is considered. The new approximation sums up beneﬁts of work on strong approximation spaces with the signiﬁcant improvement associated to the “bubble” space as discussed in the next section. Improvement of strongly discontinuous approximations Partition of unity alteration As discussed in “Numerical analysis of strongly discontinuous approximations with quadratic elements” section, the use of vertex node and middle node shape functions (resp. χ Φ and χ Φ in Fig. 16) leads to an incorrect asymptotic behavior of strongly + S + M discontinuous formulations, both for Bernstein and Lagrange polynomials. −3 When getting rid of the middle node d.o.f. χ Φ for ε around 10 , the condition + M number decreases sharply (Fig. 28). This threshold value is close to estimates of [18,24], but we stress on the fact that only the middle d.o.f. is removed and that the position of the interface is not shifted. This alteration process allows the analysis to proceed beyond −12 ε = 10 as with linear elements. Let us precise the diﬀerences between the removal of the middle node and straightfor- ward criteria from “Strong discontinuity approximation conditioning” section: • The ﬁt-to-vertex and volume criterion remove both the vertex node and the middle node d.o.f. (χ Φ and χ Φ ). Then, with the ﬁt-to-vertex, the interface position + S + M is switched to the closest node. • With the new strategy, only the minimal information is removed from the approxima- tion (a single d.o.f., the nearest to the interface) with a tuned threshold. The interface is not switched. However, the alteration of the partition of unity has a noticeable impact on enrichment [33] with Lagrange polynomials (Fig. 28): the H1-error increases sharply after the alter- ation, instead of decreasing as observed with other formulations. This behavior might be explained by the partition of unity and the nature of the solution represented here: Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 32 of 51 Fig. 28 H1-error with quadratic Lagrange polynomials. The enrichments are alternated through the removal −3 of the middle node beyond a threshold eps of ε = 10 • The exact solution is continuous, •In [7,26] discontinuous functions are added to continuous functions, so that the removal of the discontinuous middle node d.o.f. does not alter the partition of unity of continuous space functions [9] and consequently the representation of continuous functions. This is in contrast with [33], as the removal of the enriched middle node d.o.f. prevents the representation of continuous functions on some parts of the domain where the partition of unity is broken. Extrapolation in 3D Let us consider the patch-test of [24]. It is a 3D block of side 4m deﬁned by = [− 2, + 2]× [− 2, + 2]×[− 2, + 2] and split by an inclined interface. The domain is meshed with regular quadratic hexahedral elements (4 × 4 × 4 elements). A multi-axial loading is applied. Given the elastic behavior, the stress is homogeneous within the block at a value of 10 MPa, for any given position of the interface. The interface equation is parameterized with δ as: x + y + z + δ = 0 (51) The error is analyzed with the energy norm: $ ,$ h h u − u = σ − σ : ε − ε d σ : εd (52) energy As the Hansbo et al. [33] enrichment performs as well as the “shifted” enrichment, only the “shifted” enrichment is tested here. The strategies discussed in 1D are extrapolated in 3D (an example is given Fig. 29): • Preconditioning strategies are unchanged (Béchet et al. pre-conditioner [21]with X-FEM and diagonal scaling with shifted enrichment), Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 33 of 51 Fig. 29 3D test-case with quadratic elements. The strategies to improve the condition number are assessed on a 3D patch-test • The elimination threshold for the middle nodes, is ﬁxed at a ratio of relative distance −3 to the vertex of 10 along the edge (see on the 1D example illustrated by Fig. 16), or, could be expressed with a volume ratio criterion [24]: - . measure( ∩ Supp(Φ ) ), + M −9 min ≤ 10 measure( Supp(Φ ) ) (53) measure( ∩ Supp(Φ ) ) − M On (Fig. 30), the shifted formulation is more accurate than the X-FEM one, as observed in the 1D case. The diﬀerence in accuracy reaches a 10 factor: accuracy issues observed in 1D with X-FEM are heightened in 3D. Furthermore, the volume ratio criterion on the middle node is very satisfactory. The results improve as soon as the elimination process is activated around a node to interface −3 distance (here δ/ 3) of 10 m. Synthesis Although, all three strongly discontinuous formulations considered in the paper represent the same approximation space, they do not have the same numerical performance. From the results above, domain and shifted enrichments are less sensitive to the position of the interface than X-FEM. With quadratic elements, the three formulations need a special care regarding the asymptotic behavior of the shape functions as shown in “Numerical Fig. 30 Behavior of strongly discontinuous enrichments in 3D with the corrections studied previously Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 34 of 51 analysis of strongly discontinuous approximations with quadratic elements” section. The X-FEM jump formulation is far less accurate in 3D. From a practical point of view, the “shifted” formulation seems to oﬀer the best com- promise, because it provides: • A simpler pre-conditioner than X-FEM, which decreases the computational cost, • Less discontinuous d.o.f. than [33], which eases the implementation and improves the accuracy of the formulation in case of elimination of d.o.f. with Lagrange polynomials, • A good accuracy in 3D with quadratic elements. Singular approximation at a crack tip New “bubble” approximation space We introduce a simpler space than the one of Gupta et al. [25,31], to reduce its too many additional Heaviside d.o.f. As a matter of fact Gupta et al. approach requires up to 12 Heaviside d.o.f. per node instead of 3 so that it is diﬃcult to implement and may lead to ill conditioning. In the elements where K is discontinuous, a technique similar to the ghost node interpolation used in [41] is preferred. The principle consists in using two continuous interpolations by prolongation of the one on and the one on to the whole domain + − (an example is given Fig. 31). Hence, two continuous extensions of K (r, θ) are introduced: • K (r, |θ |) overlaps with K (r, θ)on , α α • K (r, − |θ |) overlaps with K (r, θ)on . α α In the elements where K has a singular point, the interpolation of K is smoothened, α α so that subsequent interpolation errors of the singular function K do not pollute the accuracy of the displacement ﬁeld as observed in [25]. Fig. 31 Comparison between double node interpolation and ghost node interpolation for an arbitrary function. The ghost node method uses the same element to interpolate both extrapolated branches. In contrast, the double node reﬁnes the mesh around the discontinuity. Thus, the double node is twice more accurate than the ghost node, but it requires a meshing of the interface (a strategy not relevant in the framework of X-FEM/GFEM) Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 35 of 51 Let n be the order of the interpolation used for the crack tip degrees of freedom. I is CT the corresponding set of nodes enriched by the vectorial asymptotic functions i.e. I CT,1 corresponds to vertex nodes and I contains both vertex nodes and middle nodes. If CT,2 we consider a quadratic interpolation for the crack tip degrees of freedom, we have n = 2 and I = I .Let I be the associated subset of I with singular points i.e. the set of CT,2 CT T CT nodes such that a crack-tip belongs to their support. The nodes of I are not taken into account in the interpolation of K , so that the interpolation of K is smooth when the α α polar coordinate r goes to zero. Let m be the order of the interpolation of K , which can be diﬀerent from n. The new interpolation operator reads: K = K (r , θ ) k k k,m α α k∈ I ∪δI /I /I {{ CT,m CT,m} T } H + K (r , |θ |) χ + + K (r , − |θ |) χ − (54) α k k k,m α k k k,m k∈I ∩{{I ∪δI }/I } H CT,m CT,m T where are the shape functions corresponding to order m interpolation and I is k,m CT,m the subset of I corresponding to order m interpolation, i.e. I is the restriction CT CT,1 of I to vertex nodes and I contains both vertex nodes and middle nodes of I . CT CT,2 CT Since we consider quadratic approximation in this section, we have I = I . The set CT,2 CT I ∪ δI contains all the nodes lying in the support of a node of I . CT,m CT,m CT This interpolation should be interpreted, for an element which is split by the interface (i.e. an element that does not contain a crack-tip), as follows • Assuming the evaluation point (gauss point) is located on , (K ) = K (r , |θ |) (55) α α k k k,m k∈{I ∪δI }/I CT,m CT,m T • Assuming the evaluation point (gauss point) is located on , | | (K ) = K (r , − θ ) (56) k k k,m α α k∈{I ∪δI }/I CT,m CT,m T The new “bubble” approximation space, based on the new interpolation operator, is: ⎧ ⎫ d d ⎪ ⎪ ⎪ ⎪ v = a Φ E + b (H − H(x )) Φ E ⎪ ⎪ i,j i i,j i i j j ⎨ ⎬ i∈I j=1 i∈I j=1 h H V = (57) bub,m,n ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ + c K − K ⎩ i,α i,n ⎭ α α i∈I α=1 CT,n where Φ are the shape functions corresponding to the order of the elements, i.e. p = 2 in this section and n refers to the order of the shape functions used for the vectorial i,n enrichment discretization. Remark • Let us stress on the fact that the ghost node interpolation does not introduce additional nodes. The two branches of discontinuous functions are interpolated over the element through extrapolation of branches (Fig. 31); • The diﬀerence between the original SGFEM and the “new” bubble space V bub,m,n lies in the behavior of the interpolation in the bandwidth of elements where the enrichment function becomes unsmooth or discontinuous. The original SGFEM uses the same interpolation operator whether the enrichment function is smooth or dis- continuous. The new “bubble” space replaces the enrichment function with smooth Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 36 of 51 continuous extensions in the bandwidth of elements where the enrichment func- tion becomes discontinuous. Thanks to the ghost node interpolation, both smooth extensions are interpolated without additional node (Fig. 31). Comparison with previous approximations The mode I problem of “Numerical behavior of singular approximations at the crack tip” section is considered here, with a similar enrichment radius of r = 0.1. In this section, we focus on the convergence analysis carried out with quadratic elements. On Fig. 32, the “bubble” enrichment has a better condition number than the linear vectorial enrichment. On Fig. 33, all enrichment strategies have optimal convergence rates in energy norm. The bubble transformation of the singular function preserves the optimality of the vectorial formulation. h h Nonetheless, bubble spaces V and V areabout ﬁvetimes more bub,m = 2,n = 1 bub, m = 2,n = 2 h h accurate than the bubble space V . The bubble space V has bub, m = 1,n = 1 bub, m = 1,n = 1 the same accuracy than X-FEM vectorial enrichment. Clearly, the bubble spaces with Fig. 32 Conditioning of “bubble” enrichment vs X-FEM vectorial enrichment Fig. 33 convergence analysis of “bubble” enrichment strategies with vectorial enrichment Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 37 of 51 quadratic interpolation of the singular function outperform X-FEM and the linear bub- ble space. As a matter of fact, among the bubble spaces with quadratic interpolation of h h singular function denoted V and V , we recommend to choose bub, m = 2,n = 1 bub, m = 2,n = 2 n = 1 (linear enrichment for the singular part) but m = 2 (quadratic interpolation of the singular function) to have the best accuracy and lowest condition number with an order with respect to mesh reﬁnement around 3.25 for an optimal expected value of 2 [31,32,38] (see Fig. 32). Remark • V is not very interesting because it has a poor numerical behav- bub, m = 1,n = 2 ior. The shift with a lower order of interpolation has not been considered even in the original method [25]. • The good results detailed in this section with quadratic elements, are conﬁrmed with linear elements. In this case the condition number for the bubble space V bub, m = 1,n = 1 evolves with an order with respect to mesh reﬁnement around 2 which corresponds to the optimal expected value (see Fig. 34 related to Figs. 21, 22 and Fig. 35 related to Figs. 25, 26 of “Numerical analysis of crack approximations with linear elements: crack opening in mode 1” section). In addition, an accuracy study of stress intensity factors is also shown with convergence orders with respect to mesh reﬁnement close to the theoretical value of 2 (Fig. 36). Fig. 34 Convergence and conditioning analysis with linear elements for a horizontal crack Fig. 35 Convergence and conditioning analysis with linear elements for inclined crack at an angle of 44.9 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 38 of 51 Fig. 36 Convergence analysis of stress intensity factor for the mode I problem. Stress Intensity Factors (SIF) are estimated with a domain integral method (G-theta method [42,43]). The “bubble” formulation with quadratic elements reaches the limit of numerical accuracy of the estimator used so that a convergence analysis cannot be performed here. For formulations with linear elements, the convergence rates are close to the optimal rate of convergence with mesh reﬁnement, around 2 Extension to weak discontinuities Similarly to the analyses performed for strong discontinuities and singular enrichment, this section will describe several enrichments for weak discontinuities found in the lit- erature. Those are used in case of bi-materials. We will show that getting a result with correct accuracy does not depend only on conditioning issues but also on the quality of the approximation space chosen to represent continuous ﬁelds with discontinuous gradients. These two separate issues are often mixed in the literature as it will be shown later in this section. Enrichment functions for weak discontinuities The general enrichment in case of weak discontinuity takes the following form: u (x) = a Φ (x) + b F (x) Φ (x) (58) i i j j i∈I j∈I with the same notations than in “Strong discontinuity representation” section. The ﬁrst sum in the right hand part of the previous equation corresponds to the standard ﬁnite element approximation while the second one stands for the enriched term. Diﬀerent choices are present in the literature for the function F x and the ones in ( ) [14,27,28] are possible (Table 3). Some of them are related: as a matter of fact, the choice by [28] is a shifted version [26]of[14]and [27] is a bubble transformation of [14]. Table 3 Weak discontinuity enrichments comparison Sukumar [14]Moës[27] Belytschko [28] F(x) lsn Φ x |lsn | Φ x − lsn Φ x lsn Φ x − lsn i i ( ) i i ( ) i i ( ) i i ( ) j i∈I i∈I i∈I i∈I H H H H Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 39 of 51 We analyze in the following the inﬂuence of these choices on the overall behavior of the solution with respect to the satisfaction of simple patch tests convergence and accuracy properties. Remark Sukumar’s and Belytschko’s enrichment functions yield almost the same approx- imation space. It can be shown straightforwardly through a change of variables as in “Strong discontinuity approximation conditioning” section. Let a and b be the classical j j and enriched degrees of freedom for Belytchko’s enrichment and a and b refer to the j j classical and enriched degrees of freedom for Sukumar’s enrichment. We have for j ∈ I a = a − lsn b j j j (59) b = b Thus, in the next sections both approximations may display similar numerical properties. Note also that the zero of the level set term lsn Φ (x) which appearsin[14,27,28] i i i∈I represents a discretization of the interface. In the following we took the same orders of interpolation for the functions F(x) and the level sets, which seemed to be quite natural. Numerical analysis of weakly discontinuous approximation In order to see if conditioning issues are also relevant for weak discontinuities we per- formed a numerical analysis similar to the one of “Numerical behavior of strongly dis- continuous approximations” section, relying on a one dimensional bi-material problem. We will show that even if these conditioning issues can be prevented by a change of for- mulation, convergence issues can still be observed, due to the fact that simple patch tests cannot be represented correctly for linear or quadratic elements. Analytical solution for the one dimensional bi-material problem To capture piecewise linear displacements as in “Numerical analysis of strongly discon- tinuous approximations with linear elements” section, we simulate a one-dimensional bi-material problem of length L in a two dimensional domain since our numerical soft- ware is 2D and 3D and as was done in [14]. To ease calculations, we reduce it to the one-dimensional bi-material problem. • For the one dimensional bi-material bar, the solution veriﬁes the following diﬀerential equation d u = 0, ∀x ∈ = [0,x [∪]x , L] , so that, 0 0 dx α x + β,if x ≤ x u(x) = , where x is the location of the interface, α x + β,if x ≥ x • With arbitrary Dirichlet boundary conditions we choose to impose on the central element nodes for the sake of simplicity, u(0) = 0,u(L) = u • Material properties are given by: E,if x ≤ x 1 0 E (x) = , where x is the location of the interface, i 0 E,if x ≥ x 2 0 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 40 of 51 • Assuming small strains, the expected solution satisfying the boundary conditions above is: E u − E u 2 2 x, if x ≤ x so that α = E L+(E −E )x E L+(E −E )x 1 2 1 0 1 2 1 0 u(x) = (60) E x+(E −E )x E u ⎩ 1 2 1 0 + 1 u, if x ≥ x so that α = E L+(E −E )x E L+(E −E )x 1 2 1 0 1 2 1 0 Numerical analysis of weakly discontinuous approximations with linear elements In the numerical analysis of “Numerical analysis of weakly discontinuous approximations with linear elements” section, to ease the calculations, we consider = [− 1, 2] with the following Dirichlet boundary conditions u (-1) = 0, u (2) = u¯ so that L = 3m.The interface is = {x }. The weak form of the problem becomes: Find: − + u ∈ : w ∈ H () /w(−1) = 0,w(2) = u, ¯ w(x ) = w x 0 0 − + ∀v ∈ : w ∈ H () /w(−1) = 0,w(2) = 0,w(x ) = w x 0 0 0 E u v dx = 0 −1 → (61) − + E u x = E u x 1 2 0 0 The X-FEM discrete linear space is: w/w = a + a + b F + a + b F + a 1 1 2 2 2 2 3 3 3 3 4 4 : (62) a = 0 a = u¯ 1 4 Then, the 6 × 6 matrix associated with the discretization of the weak form on the X-FEM space of linear function with Sukumar et al. formulation [14]is: ⎛ ⎞ ⎛ ⎞ E −E −E (1 − ε)0 0 0 1 1 1 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ −E E (2 − ε) + E ε 2E (1 − ε) −E (1 − ε) − E ε −E ε 0 ⎜ ⎟ 1 1 2 1 1 2 2 ⎜ ⎟ a ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ 3 3 ⎟ 3 3 (2−ε) (1−ε) ⎜ ⎟ ε ε ⎜ ⎟ −E (1 − ε)2E (1 − ε) E + E −E (1 − ε) −E − E 0 b 1 1 1 2 1 1 2 ⎜ 2 ⎟ 3 3 3 3 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 0 −E (1 − ε) − E ε −E (1 − ε) E (1 − ε) + E (ε + 1) 2E ε −E ⎜ a ⎟ 1 2 1 1 2 2 2 3 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 3 3 3 3 ⎜ ⎟ (1−ε) ε (1−ε) (ε+1) ⎜ ⎟ b 0 −E ε −E − E 2E ε E + E −E ε 2 1 2 2 1 2 2 ⎝ ⎠ 3 3 3 3 ⎝ ⎠ 00 0 −E −E ε E 2 2 2 With the formulation of Belytschko et al. [28] it becomes: ⎛ ⎞ ⎛ ⎞ E −E 00 0 0 1 1 ⎜ ⎟ ⎜ ⎟ −E E (2 − ε) + E ε (E − E )(1 − ε)ε −E (1 − ε) − E ε (E − E )(1 − ε)ε 0 1 1 2 1 2 1 2 1 2 ⎜ ⎟ ⎜ ⎟ a ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ E E E E ⎟ 1 2 1 2 ⎜ ⎟ 0(E − E )(1 − ε)ε (u + 2) + v (E − E )(1 − ε)ε u + v 0 1 2 1 1 2 1 2 2 ⎜ ⎟ 3 3 3 3 b ⎜ 2 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ 0 −E (1 − ε) − E ε (E − E )(1 − ε)ε E (1 − ε) + E (ε + 1) (E − E )(1 − ε)ε −E ⎟ 1 2 2 1 1 2 2 1 2 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ E E E E 1 2 1 2 ⎜ ⎟ 0(E − E )(1 − ε)ε u + v (E − E )(1 − ε)ε (u + 1) + (v + 1) 0 1 2 2 2 2 1 1 1 ⎜ ⎟ 3 3 3 3 3 ⎝ ⎠ ⎝ ⎠ 00 0 −E 0 E 2 2 u (ε) =−ε [ε (4ε − 6) + 3] v (ε) = ε [ε (4ε − 6) + 3] 1 1 where: u (ε) = −2ε + 3ε-1 v (ε) = ε [ε (2ε − 6) + 3] 2 2 Finally, with Moës et al. formulation [27], the discrete matrix is: ⎛ ⎞ ⎛ ⎞ E −E 000 0 a 1 1 1 ⎜ ⎟ ⎜ ⎟ 2 2 ⎜ ⎟ ⎜ ⎟ −E E (2 − ε) + E ε 2(E − E )(1 − ε)ε −E (1 − ε) − E ε 2(E − E )(1 − ε) ε 0 a 1 1 2 2 1 1 2 2 1 2 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 2 2 02(E − E )(1 − ε)ε E f + E f 2(E − E )(1 − ε)ε E h + E h 0 b ⎜ 2 1 1 1 2 2 1 2 1 1 2 2 ⎟ ⎜ 2 ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 2 2 ⎜ 0 −E (1 − ε) − E ε 2(E − E )(1 − ε)ε E (1 − ε) + E (ε + 1) 2(E − E )(1 − ε) ε −E ⎟ ⎜ a ⎟ 1 2 1 2 1 2 1 2 2 3 ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ 2 2 ⎜ ⎟ ⎜ ⎟ 02(E − E )(1 − ε) ε E h + E h 2(E − E )(1 − ε) ε E g + E g 0 b 2 1 1 1 2 2 1 2 1 1 2 2 3 ⎝ ⎠ ⎝ ⎠ 00 0 −E 0 E 2 2 4 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 41 of 51 where: 4 16 2 2 3 f = 1 − ε ε ε 4ε − 2 + 1 ; f = 1 − ε ε ; ( ) [ ( ) ] ( ) 1 2 3 3 2 2 h = (1 − ε) ε (−1 + 4ε) 4 16 4 2 2 3 2 2 h = (1 − ε) ε (3 − 4ε); g = (1 − ε) ε ; g = (1 − ε) ε(ε(4ε − 6) + 3) 2 1 2 3 3 3 For all cases, we have taken ε = 1 − x , where x denotes the interface position. 0 0 On considering the 4 × 4 sub-matrix that is free from boundary conditions, the diag- onal pre-conditioner can easily scale the formulation of Moës et al. [27] as entire rows (symmetrically columns) vanishes when ε = 0or ε = 1(seeFig. 37). For the enrichments of Sukumar et al. [14] and Belytschko et al. [28], no conditioning issues can be observed, in case of linear approximation. However, we will see in the next paragraph that these latest two formulations are far from optimal with respect to convergence issues, being unable to represent simple patch tests. Ability of the enrichment to catch piecewise linear and quadratic solution First of all, we would like to test that the proposed enrichments (see Table 4)are able to capture piecewise linear and quadratic solutions, so that the enrichment function for the material interface problem preserves the equivalence between the FEM and X-FEM discretized spaces. For that purpose, we will consider the previous problem of the one- dimensional bi-material bar for the linear case [14]. Then to capture piecewise quadratic displacements we adapt the bi-material volumic load test of [27] using a constant volumic load, instead of a quadratic one. To simplify the calculations, the bar is subdivided into three elements (see Fig. 8)of equal length L, in order to have cut and uncut elements. The central element which carries the interface is such that 0 ≤ x ≤ L. Dirichlet boundary conditions are applied on this central element such that u(0) = 0,u (L) = u. The linear shape functions we use in x x the central element are of the form = 1 − , = , while the quadratic ones are 1 2 L L taken as: Fig. 37 Conditioning of X-FEM formulation with quadratic elements for weak discontinuities representation. With X-FEM’s ridge enrichment functions no conditioning issues are observed through the use of a diagonal pre-conditioner Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 42 of 51 1 2 x x 4x x x 2x = 1 − 2 1 − , = 1 − , =− 1 − 1 2 3 L L L L L L Unknowns that we search to recover the linear and quadratic patch test solutions are the degrees of freedom a ,a ,b ,b and a ,a ,a ,b ,b ,b ,in(58) respectively. 1 2 1 2 1 2 3 1 2 3 For Sukumar and Belytschko formulations, the X-FEM blending elements which are adjacent to those completely cut by the interface have a quadratic interpolation [even with linear elements, due to the interpolation (58) and the expression of F(x)inTable 4] or cubic interpolation [with quadratic elements and the expressions (58)and F(x)] with just one active b degree of freedom (see Fig. 8, for the left and right blending elements). Since the patch test is linear, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence linear patch tests cannot be satisﬁed and are just approximated by these formulations, which will aﬀect the energy norm as noticed by [14]. Note that Moës formulation with the same order of interpolation for F(x) and the shape functions Φ(x) is able to catch the linear patch test solution, and that Moës formulation with the linear interpolation of F(x) is also able to catch the solution for quadratic or linear shape functions. Table 5 gives a summary of the principal results obtained for the linear patch test. Table 6 represents the relative error in energy norm for this bi-material bar under traction for diﬀerent positions of the interface. One can see that the enrichments of Sukumar and Belytschko, even when quadratic (see Table 7), do not yield proper results due to the transition of F(x) in blending elements [26], while the problem is avoided with Moës enrichment. Moreover, as we cannot recover the exact linear solution for this patch test, an order of convergence of 0.5 in energy norm with respect to mesh reﬁnement is obtained for Sukumar and Belytschko enrichments (Fig. 38), which is consistent with the previous work of Moës et al. [44] extended to X-FEM with an enrichment that does not allow to catch the discontinuity. Analytical solution for the volumic load problem Consider a plate of side L with two materials E and E with ν = ν = 0, the interface 1 2 1 2 being localized at x = x in the plate, as in Fig. 39. Plain strain conditions are assumed. The plate is embedded on the left side and a variable volumic force is imposed on the plate 3 3 3 3 along the x direction, given by 3f 3 = f . The case is computed in two dimensions since our numerical software works only in two and three dimensions, using a similar approach to the one in [14]. Table 4 Expressions of F(x) for weakly discontinuous enrichments with linear and quadratic approximations resulting from Table 3 for the one dimensional bi-material bar problem Sukumar [14]Moes[27] Belytschko [28] xx | | | || | | | Linear F(x) x − x x + x − 2 − x − x x − x − x − x 0 0 0 0 J 0 | | | | | | Quadratic F(x) x − x • If L/2 ≥ x x − x − x − x 0 0 0 J 0 x + (L − 6x ) 0 0 | | + 4x − x − x 0 0 • If L/2 ≤ x x + (2x − 3L) 0 0 | | + 4 (L − x ) − x − x 0 0 L Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 43 of 51 ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎧ ⎪ ⎨ ⎪ ⎩ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Table 5 Coeﬃcients used to capture the solution with the ﬁnite element approximation for the linear patch test Sukumar [14]Moës[27] Belytchko [28] Linear φ(x) and The analytical solution a = 0 The analytical solution linear F(x) cannot be represented. a = u¯ cannot be represented. (E −E )¯u 2 1 Incompatible with b = b = Incompatible with 1 2 2(E L+(E −E )x ) 1 2 1 0 blending elements. blending elements. a = 0 u¯ (E −E )x u¯ 2 1 0 a = + 2 2(E L+(E −E )x ) 1 2 1 0 • if L/2 ≥ x a = u¯ (E2−E1)¯u b = b = b = 1 2 3 2(E L+(E −E )x ) 1 2 1 0 Quadratic φ(x) and The analytical solution The analytical solution a = 0 quadratic F(x) cannot be represented. cannot be represented. u¯ (E − E )(L − x )¯u Incompatible with 2 1 0 Incompatible with a = + 2 2(E L + (E − E )x ) blending elements. 1 2 1 0 blending elements. • if L/2 ≤ x a = u¯ (E − E )¯u 2 1 b = b = b = 1 2 3 2(E L + (E − E )x ) 1 2 1 0 a = 0 Quadratic φ(x) and The analytical solution (E − E )x u¯ + E Lu¯ u¯ The analytical solution 2 1 0 1 a = = linear F(x) cannot be represented. 2(E L + (E − E )x ) 2 cannot be represented. 1 2 1 0 (E − E )x u¯ + E Lu¯ Incompatible with 2 1 0 1 Incompatible with a = = u¯ (E L + (E − E )x ) 1 2 1 0 blending elements. blending elements. (E − E )¯u 2 1 b = b = b = 1 2 3 2(E L + (E − E )x ) 1 2 1 0 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 44 of 51 Table 6 Comparison of relative error in energy norm for diﬀerent interface locations between the diﬀerent enrichment functions for linear elements used in the bi-material bar −6 problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa,L=25mand u¯ = 3 · 10 m) Interface location Sukumar Moës Belytschko 12.3 2.38688E−01 3.26389E−15 2.38688E−01 12.4 2.42403E−01 4.30236E−15 2.42403E−01 12.5 2.45626E−01 2.82096E−15 2.45626E−01 12.6 2.48370E−01 5.19794E−15 2.48370E−01 12.7 2.50650E−01 6.25056E−15 2.50650E−01 Table 7 Comparison of relative error in energy norm for diﬀerent interface locations between the diﬀerent enrichment functions for quadratic elements used in the bi-material bar problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa,L=25mand −6 u¯ = 3 · 10 m) Interface location Sukumar Moës Belytschko 12.3 8.37168E−02 2.19329E−14 8.37168E−02 12.4 8.02193E−02 2.72702E−14 8.02193E−02 12.5 7.62151E−02 3.39596E−12 7.62151E−02 12.6 7.17862E−02 1.85735E−14 7.17862E−02 12.7 6.70550E−02 1.43340E−14 6.70550E−02 The exact displacements in the two materials are given by: 1 2 f x u (x) =− − Lx x if x ≤ x 1 0 E 2 4 - . 5 1 2 1 2 2 2 f x 1 1 u (x) =− − Lx + f − Lx − x if x ≥ x 2 0 0 E 2 2 E E 2 1 2 No linear element with linear F(x) is able to capture this quadratic solution directly. For Sukumar and Belytschko formulations with quadratic elements, the xfem blending elements which are adjacent to those completely cut by the interface give a cubic inter- polation (linear F(x)) or order four interpolation (quadratic F(x)) with just one b degree of freedom (see Fig. 8, for the left and right elements). Since the patch test is quadratic, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence Fig. 38 Convergence in energy norm at a rate of 0.5 [44] for Sukumar and Belytschko enrichments. These −6 results are obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and u¯ = 3.10 m Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 45 of 51 Fig. 39 Plate with a material interface and a variable volumic load quadratic patch tests cannot be satisﬁed and are just approximated by these formulations, which will aﬀect the energy norm as noticed by [14]. Note that Moës formulation with quadratic shape functions Φ(x) and linear F(x) is able to catch the quadratic patch test solution, and that Moës formulation with quadratic shape functions Φ(x) and quadratic F(x) is not able to catch it. Table 7 gives a summary of the principal results obtained for the quadratic patch test. Figure 40 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh reﬁnement for the quadratic patch test, using linear shape functions Φ(x) and linear F(x) and quadratic shape functions Φ(x) and quadratic F(x). Using linear elements we obtain the optimal convergence rates of 1. Using quadratic elements we obtain 1.5. In this later case, a convergence rate of 2 would be expected since the analytical solution does not lie in the ﬁnite element space (cf. Table 8). We will not prove formally that the obtained convergence rates are in fact the expected ones, nevertheless we can give some arguments that those rates are indeed the expected ones. In the quadratic patch test case, the ﬁrst and second derivatives of the analytical solu- tion are discontinuous. The ﬁnite element space with enrichment we use with quadratic F(x) and quadratic shape functions Φ(x) does not allow to capture this solution. How- ever it allows to capture the linear patch test solution with a discontinuous ﬁrst deriva- tive. We can deduce from this analysis that its rate of convergence on energy will be greater than 1 but lower than 2. Then the loss of 0.5 in convergence rate observed numerically for quadratic elements, with a value of 1.5 observed numerically instead of an expected 2, is coherent with the analysis for linear elements (0.5 instead of 1 when elements do not catch the discontinuity on the ﬁrst derivative) performed by Moës et al. [44]. From this analysis, we recommend to use Moës formulation with quadratic elements and linear F(x). Convergence analysis for weak discontinuities In this paragraph we will investigate the convergence of our X-FEM formulations dedi- cated to weak discontinuities when the solution is not in the approximation space. Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 46 of 51 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Table 8 coeﬃcients used to capture the solution with the ﬁnite element approximation for the quadratic patch test Sukumar [14]Moës[27] Belytchko [28] Linear φ(x)and The analytical solution cannot be represented. Incompatibility for the central element. linear F(x) Quadratic φ(x)and The analytical solution No solution. The analytical solution quadratic F(x) cannot be represented. Incompatibility for the central element. cannot be represented. Incompatible with blending elements. Incompatible with blending elements. a = 0 1 2 f 5 a = x L − x + 2 0 4E 2 1 2 f 3 5 2 2 L − x L + x 4E 2 2 f f a = x (2L − x ) + (x − L) 3 0 0 0 Quadratic φ(x)and The analytical solution The analytical solution 2E 2E 1 2 1 2 linear F(x) cannot be represented. cannot be represented. f 1 1 b = − (2L − x ) 1 0 Incompatible with blending elements. 4 E E Incompatible with blending elements. 1 2 1 2 f 1 1 b = − (3L − 2x ) 2 0 8 E E 1 2 1 2 f 1 1 b = − (L − x ) 3 0 4 E E 1 2 Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 47 of 51 Fig. 40 Convergence rates on the relative error in terms of the energy norm with respect to mesh reﬁnement for the quadratic patch test, using linear shape functions and a linear enrichment function (linear, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic) Volumic load problem with straight interface The test we propose in the following was already investigated by Moës et al. [18]. Let us consider again the squared plate of Fig. 39 with side L = 2 m separated in its middle by two materials E = 1MPa and E = 10 MPa and with ν = ν = 0. This time the variable 1 2 1 2 3 3 3 3 2 volumic force is f = f imposed on the plate along the x direction. 3 3 The exact displacements in the two materials are: 1 2 f x u (x) = L x − x 3E L 4 6 1 2 1 27 4 2 f x 31f L 1 1 u (x) = L x − + − x 3E L 4 192 E E 2 1 2 We considered an unstructured mesh with 5, 10, 20 and 40 triangular elements on each side. Figure 41 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh reﬁnement for the volumic load problem, using linear shape functions Φ(x) and linear F(x), quadratic shape functions Φ(x) and linear F(x)and quadratic shape functions Φ(x) and quadratic F(x). Using linear elements or quadratic elements and linear enrichment function we obtain the optimal convergence rates. Using quadratic elements and quadratic enrichment function we obtain 1.5. Note that the second derivative of the analytical displacement is discontinuous in this case too. Circular inclusion with imposed displacement To investigate the behavior of weak discontinuity enrichments in presence of curved interfaces, Moës et al. [18] suggested the following test case. A circular inclusion of radius a = 0.4 m is placed in the center of a material . A linear displacement is imposed on the boundary (r = b = 2 m). The same material parameters are chosen: E = 1MPa, ν = 0.25 and E = 10 MPa, ν = 0.3. This problem has an analytical 1 1 2 2 solution provided in [18]. The analytical displacement reads: ⎧ ( ) 2 2 b b 1 − α + , if 0 ≤ r ≤ a 2 2 a a u (r) = , 2 2 b b r − α + , if a ≤ r ≤ b 2 2 a a Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 48 of 51 Fig. 41 Convergence rates on the relative error in terms of the energy for the volumic load problem. The enrichment strategies are: linear shape functions and linear enrichment function (linear, linear), quadratic shape functions and a linear enrichment function (quadratic, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic) Fig. 42 Convergence rates on the relative errors with respect to mesh reﬁnement for the circular inclusion problem. The enrichment strategies are: linear shape functions and linear enrichment function (weak, linear, linear), quadratic shape functions and linear enrichment (weak, quadratic, linear) and quadratic shape functions and quadratic enrichment function (weak, quadratic, quadratic). Two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface have also been investigated: one with quadratic shape functions and a linear approximation of the level set (strong, quadratic, linear) and one with quadratic shape functions and a quadratic approximation of the level set (strong, quadratic, quadratic) with: (λ + μ + μ ) b 1 1 2 α = 2 2 2 2 λ + μ a + λ + μ b − a + μ b ( ) ( ) 2 2 1 1 2 The domain of computation is a square of side L = 2 m. We considered an unstructured mesh with 5, 10, 20, 40, 80 and 160 triangular elements by side. Figure 42 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh reﬁnement for circular inclusion, using linear shape functions Φ(x) and linear F(x), quadratic shape functions Φ(x) and linear F(x) and quadratic shape functions Φ(x)and quadratic F(x). Using linear elements we obtain the optimal convergence rate. Using quadratic elements we obtain 1.5. In the case of quadratic F(x), this behavior arises from the discontinuity of the second derivative of the analytical displacement. The case of linear Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 49 of 51 F(x) is slightly diﬀerent, due to its link with the geometrical approximation of the level sets and more particularly of the interface established at the end of “Enrichment functions for weak discontinuities” section. Ferté et al. [45,46] showed for strong discontinuities that the approximation of the level set impacts the expected theoretical convergence rates, in the case of quadratic elements: 2 is expected if a quadratic approximation of the level set is used (equivalently quadratic F(x)), while 1.5 is expected if a linear representation of the level set is used (equivalently linear F(x)). To illustrate this behavior, we also considered two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface: one with quadratic shape functions and a linear approximation of the level set and one with quadratic shape functions and a quadratic approximation of the level set. The convergence rates obtained for these two models are also shown on Fig. 42. Finally for weak discontinuities, optimal rates of convergence in energy norm are diﬃ- cult to obtain for quadratic elements. If linear approximations of the level sets (equivalently linear F(x)) are used associated to a quadratic approximation of the enrichment, conver- gence rates are limited due to errors on the geometry. If quadratic approximations of the level sets (equivalently quadratic F(x)) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the approximation spaces. It appears that both convergence rates end up being the same, with a loss of 0.5 in energy norm with respect to optimal orders of convergence of 2. Optimal orders of convergence are recovered when using a strong discontinuity approximation space with lagrangian multipliers to represent weak discontinuities. Conclusion This paper outlines numerical issues around popular approximation spaces to model weak discontinuities or strong discontinuities and cracks, particularly with quadratic elements. We analyzed those issues on a few benchmarks. Those benchmarks revealed that popular enrichment strategies do not behave well with asymptotic conﬁgurations (when the dis- continuity gets close to a node of the mesh, typically for a distance below 10% of the length of the cut edge) for strongly discontinuous enrichments and that the singular enrichment could be improved with respect to conditioning and accuracy. Instead of working on complex and expensive strategies to solve those issues for strongly discontinuous enrichments, we preferred to apply a slight reshape of the approximation spaces, as it is in our sense, the most eﬀective approach to deal with those issues for industrial applications. In fact, the new approximation spaces get through the proposed benchmarks, with signiﬁcant improvement of the numerical behavior and accuracy of the formulations studied. Nonetheless, more theoretical work is needed to understand why very unsmooth “bub- ble” functions tend to give better results than “straightforward” singular functions. At least, our results conﬁrmed the numerical properties of “bubble” spaces denoted in the literature [25,31,32]. Moreover, the integration procedure needs to be improved, particu- larly in 3D with quadratic elements, to further the analysis on complex crack surface and crack front geometries. At last, in case of weakly discontinuous enrichments, it was shown that most approxi- mation spaces of the literature did not exhibit conditioning issues. However, for quadratic elements, orders of convergence in energy norm were 0.5 lower than those obtained with Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 50 of 51 an equivalent strongly discontinuous enrichment, so that we recommend using the latest approach even in the case of weak discontinuities, associated to the strategy exposed in the present work to avoid conditioning issues. Authors’ contributions MN conducted the numerical implementation of strong discontinuities approximations exposed in the paper. PM contributed to draft the manuscript and helped assessing the numerical properties of approximations spaces. NM helped strengthening the theoretical background and he provided numerous ideas to shape new approximations spaces. AM and SG participated in the study of weak discontinuities implementation and benchmarks. All authors read and approved the ﬁnal manuscript. Author details IMSIA, UMR EDF/CNRS/CEA/ENSTA ParisTech 9219, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France, Ecole Centrale de Nantes, Laboratoire de Mécanique et Matériaux, 1 Rue de La Noe, 44321 Nantes, France, Rheinisch Westfalische Technische Hochschule Aachen, German Research School for Simulation Sciences GmbH, Schinkelstrasse 2, 52062 Aachen, Nordrhein-Westfalen, Germany. Acknowledgements Not applicable. Competing interests The authors declare that they have no competing interests. Availability of data and materials Not applicable. Consent for publication Not applicable. Ethics approval and consent to participate Not applicable. Funding Not applicable. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aﬃliations. Received: 26 October 2016 Accepted: 18 September 2017 References 1. Rots JG, Nauta P, Kusters GMA, Blaauwendraad J. Smeared crack approach and fracture localization in concrete. Heron. 1985;30(1):1–48. 2. Oliver J. A consistent characteristic length for smeared cracking models. Int J Num Methods Eng. 1989;28(2):461–74. 3. Ortiz M, Leroy Y, Needleman A. A ﬁnite element method for localized failure analysis. Comput Methods Appl Mech Eng. 1987;61(2):189–214. 4. Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Num Methods Eng. 1990;30:541–64. 5. Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids. 1994;42:1397–434. 6. Ortiz M, Pandolﬁ A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Num Methods Eng. 1999;44:1267–82. 7. Moës N, Dolbow J, Belytschko T. A ﬁnite element method for crack growth without remeshing. Int J Num Methods Eng. 1999;46:135–50. 8. Duarte CA, Babuška I, Oden JT. Generalized ﬁnite element methods for three dimensional structural mechanics problems. Comput Struct. 2000;77:215–32. 9. Melenk JM, Babuška I. The partition of unity ﬁnite element method: basic theory and applications. Comput Methods Appl Mech Eng. 1996;139:289–314. 10. Westergaard HMW. Bearing pressures and cracks. J Appl Mech. 1939;6:A49–53. 11. Williams ML. On the stress distribution at the base of a stationary crack. ASME J Appl Mech. 1956;24:109–14. 12. Legrain G, Moës N, Verron E. Stress analysis around crack tips in ﬁnite strain problems using the extended ﬁnite element method. Int J Num Methods Eng. 2005;63(2):290–314. 13. Elguedj T, Gravouil A, Comberscure A. Appropriate extended functions for X-FEM simulation of plastic fracture mechanics. Comput Methods Appl Mech Eng. 2006;195:501–15. 14. Sukumar N, Chopp DL, Moes N, Belytschko T. Modeling holes and inclusions by level sets in the extended ﬁnite element method. Comput Methods Appl Mech Eng. 2001;190:6183–200. 15. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by the extended ﬁnite element and level sets—part II: level set update. Int J Num Methods Eng. 2002;53:2569–86. Ndeﬀo et al. Adv. Model. and Simul. in Eng. Sci. (2017) 4:6 Page 51 of 51 16. Duarte CA, Hamzeh ON, Liszka TJ, Tworzydlo WW. Generalized ﬁnite element method for the simulation of three dimensional crack propagation. Comput Methods Appl Mech Eng. 2001;190:2227–62. 17. Laborde P, Pommier J, Renard Y, Salaün M. High-order extended ﬁnite element method for cracked domains. Int J Num Methods Eng. 2005;64:354–81. 18. Dréau K, Chevaugeon N, Moës N. Studied X-FEM enrichment to handle material interfaces with higher order ﬁnite element. Comput Methods Appl Mech Eng. 2010;199:1922–36. 19. Cheng KW, Fries T-P. Higher-order XFEM for curved strong and weak discontinuities. Int J Num Methods Eng. 2010;82:564–90. 20. Duarte CA, Reno LG, Simone A. A high-order generalized FEM for through-the thickness branched cracks. Int J Num Methods Eng. 2007;72(3):325–51. 21. Bechet E, Minnebo H, Moes N, Burgardt B. Improved implementation and robustness study of the X-FEM method for stress analysis around cracks. Int J Num Methods Eng. 2005;64:1033–56. 22. Strouboulis T, Babuška I, Copps K. The design and analysis of Generalized Finite Element Method. Comput Methods Appl Mech Eng. 2000;181:43–69. 23. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended ﬁnite element method. Int J Num Methods Eng. 2000;48:1741–60. 24. Siavelis M, Guiton MLE, Massin P, Moës N. Large sliding contact along branched discontinuities with X-FEM. Comput Mech. 2013;52:201–19. 25. Guptaa V, Duarte CA, Babuška I, Banerjee U. A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Comput Methods Appl Mech Eng. 2013;266:23–39. 26. Fries T-P. A corrected XFEM approximation without problems in blending elements. Int J Num Methods Eng. 2008;75:503–32. 27. Moes N, Cloirec M, Cartraud P, Remacle J-F. A computational approach to handle complex microstructure geometries. Comput Methods Appl Mech Eng. 2003;192:3163–77. 28. Belytschko T, Black T, Moës N, Sukumar N, Usui S. Structured extended ﬁnite element methods of solids deﬁned by implicit surfaces. Int J Num Methods Eng. 2003;56:609–35. 29. Nicaise S, Renard Y, Chahine E. Optimal convergence analysis for the extended ﬁnite element method. Int J Num Methods Eng. 2011;86:528–48. 30. Chevaugeon N, Moës N, Minnebo H. Improved crack tip enrichment functions and integration for crack modeling using the extended ﬁnite element method. J Multiscale Comput Eng. 2013;11:597–631. 31. Guptaa V, Duarte CA, Babuška I, Banerjee U. Stable GFEM (SGFEM), improved conditioning and accuracy of GFEM/XFEM for three dimensional fracture mechanics. Comput Methods Appl Mech Eng. 2015;289:355–86. 32. Babuška I, Banerjee U. Stable generalized ﬁnite element method (SGFEM). Comput Methods Appl Mech Eng. 2012;201–204:91–111. 33. Hansbo A, Hansbo P. A ﬁnite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Methods Appl Mech Eng. 2004;193:3523–40. 34. Ling D, Yang Q, Cox B. An augmented ﬁnite element method for modeling arbitrary discontinuities in composite materials. Int J Fract. 2009;156:53–73. 35. Simone A, Duarte CA, Van der Giessen E. A generalized ﬁnite element method for polycrystals with discontinuous grain boundaries. Int J Num Methods Eng. 2006;67:1122–45. 36. Belytschko T, Black T. Elastic crack growth with minimal remeshing. Int J Num Methods Eng. 1999;45:601–20. 37. Areias P, Belytschko T. Letter to editor. Comput Methods Appl Mech Eng. 2006;195:1275–6. 38. Du Q, Wang D, Zhu L. On mesh geometry and stiﬀness matrix conditioning for general ﬁnite element spaces. SIAM J Num Anal. 2009;47(2):1421–44. 39. Rathod HT, Hiremath SV. Boundary integration of polynomials over an arbitrary linear tetrahedron in Euclidean three-dimensional space. Comput Methods Appl Mech Eng. 1998;153:81–106. 40. Mousavi SE, Sukumar N. Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons. Comput Mech. 2011;47:535–54. 41. Nistor I, Guiton MLE, Massin P, Moes N, Geniaut S. An X-FEM approach for large sliding contact along discontinuities. Int J Num Methods Eng. 2009;78:1407–35. 42. Destuynder P, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J de Mécanique Théorique Appliquée. 1983;2(1):113–35. 43. Destuynder P, Djaoua M. Sur une interprétation mathématique de l’intégrale de Rice en théorie de la rupture fragile. Math Methods Appl Sci. 1981;3:70–87. 44. Moës N, Oden JT, Zohdi TI. Investigation of the interactions between the numerical and the modeling errors in the homogenized Dirichlet projection method. Comput Methods Appl Mech Eng. 1998;159:79–101. 45. Ferté G, Massin P, Moës N. Convergence analysis of linear or quadratic X-FEM for curved free boundaries. Comput Methods Appl Mech Eng. 2014;278:794–827. 46. Ferté G, Massin P, Moës N. Interface problems with linear or quadratic X-FEM: design of a stable multiplier space and error analysis. Int J Num Methods Eng. 2014;100(11):834–70.
"Advanced Modeling and Simulation in Engineering Sciences" – Springer Journals
Published: Dec 1, 2017
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.