Qi-Zhi Zhu
email: qizhi.zhu@gmail.com

On enrichment functions in the extended finite element method

Keywords: Enrichment functions, extended finite element method (XFEM), discontinuities, domain decomposition, level set method, blending elements, imperfect interfaces

This paper presents mathematical derivation of enrichment functions in the extended finite element method (XFEM) for numerical modelling of strong and weak discontinuities. The proposed approach consists in combining the level set method with characteristic functions as well as domain decomposition and reproduction technique. We start with the simple case of a triangular linear element cut by one interface across which displacement field suffers a jump. The main steps towards the derivation of enrichment functions are as follows: (1) extension of the subfields separated by the interface to the whole element domain and definition of complementary nodal variables, (2) construction of characteristic functions for describing the geometry and physical field, (3) determination of the sets of basic nodal variables, (4) domain decompositions according to Step 3 and then reproduction of the physical field in terms of characteristic functions and nodal variables, (5) comparison of the piecewise interpolations formulated at Steps 3 and 4 with the standard XFEM form, which yields enrichment functions. In this process, the physical meanings of both the basic and complementary nodal variables are clarified, which helps to impose Dirichlet boundary conditions. Enrichment functions for weak discontinuities are constructed from deeper insights into the structure of the functions for strong discontinuities. Relationships between the two classes of functions are naturally established. Improvements upon basic enrichment functions for weak discontinuities are performed so as to achieve satisfactory convergence and accuracy. From numerical viewpoints, a simple and efficient treatment on the issue of blending elements is also proposed with implementation details. For validation purposes, applications of the derived functions to heterogeneous problems with imperfect interfaces are presented.

INTRODUCTION

Proper description of discontinuities (cracks, shear bands, slip lines, etc.) and efficient modelling of their evolution in materials and structures are still the main challenges in computational mechanics and engineering science. In that context, initiated by Belytschko and Black [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] and Moës et al. [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the extended finite element method (XFEM) has been intensively developed during the last ten years (see eg, the reviews [START_REF] Karihaloo | Modelling of stationary and growing cracks in fe framework withou remeshing: a sate-of-the-art review[END_REF] [4] [START_REF] Belytschko | A review of extended/generalized finite element methods for material modelling[END_REF] and references therein) and extensively applied to various discontinuous problems, such as propagation of cracks ( [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF] [START_REF] Rabczuk | On three-dimensional modelling of crack growth using partition of unity method[END_REF], among others), evolution of dislocations and grain boundaries (eg, [START_REF] Simone | A generalized finite element method for polycrystals with discontinuous grain boundaries[END_REF] [9] [START_REF] Menk | Crack growth calculations in solider joints based on microstructural phenomena with X-FEM[END_REF] ZHU QI-ZHI [START_REF] Mousavi | Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons[END_REF]), size-dependent effects by imperfect interfaces [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF] [13] [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF], etc.. Details on numerical implementation within the XFEM context can be found in [START_REF] Sukumar | Modelling quasi-static crack growth with the extended finite element method Part I: computer implementation[END_REF] [START_REF] Bordas | An extended finite element library[END_REF]. Theoretically, the basic idea behind the XFEM is to superpose one or several enrichment terms to the conventional continuous interpolation. Enrichment functions involved in the complementary term(s) are mandatory for describing a prescribed type of discontinuities. Both theoretical basis and numerical experiences have largely evidenced the predominant role played by enrichment functions in XFEM-framed modelling and simulation of discontinuous problems in the sense that they influence numerical accuracy and computational efficiency (convergence) directly and strongly. The current study attempts to derive mathematically and systematically enrichment functions used in the XFEM. To this end, a consistent two-step procedure is developed mainly based on the level set method, domain decomposition and reproduction technique as well as characteristic functions.

In literatures, discontinuities encountered in heterogeneous materials and structures are commonly divided into two classes: strong discontinuities where the physical field itself suffers a jump, and weak discontinuities where the gradient field suffers a jump. Enrichment functions in the XFEM are accordingly classified into two categories: functions belonging to the first one are discontinuous and dedicated to describing strong discontinuities. It is found in the literature that the Heaviside function first gained widespread use ( [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] [17] [START_REF] Iarve | Mesh independent modeling of cracks by using higher order shape functions[END_REF] [START_REF] Hettich | Interface material failure modeled by the extended finite-element method and level sets[END_REF]). From practical point of view, Fries [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF] commented that the sign step function can provide the same approximation. Mohammadi [START_REF] Mohammadi | Extended finite element method for fracture analysis of structures[END_REF] addressed in his book a review on the two types of discontinuous enrichment functions. Required for modelling weak discontinuities, functions in the second category are continuous but their derivatives are discontinuous (eg, [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] [START_REF] Farsad | An extended finite element/level set method to study surface effects on the mechanical behavior and properties of nanomaterials[END_REF]). Particularly, based on the Nitche's method, Hansbo and Hansbo [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] constructed two one-sided independent approximations for modelling both weak and strong discontinuities, which was later proved to be a linear combination of the XFEM enrichment basis of the Heaviside type [START_REF] Pedro | A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF]. Recently, some numerically determined enrichment functions have been proposed. For example, Menk and Bordas [START_REF] Menk | Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals[END_REF] realized numerical determination of enrichment functions for describing strain singularities, which enhanced significantly the convergence rates in numerical experiments. Mousavi et al. [START_REF] Mousavi1 | Higher-order extended finite elements with harmonic enrichment functions for complex crack problems[END_REF] constructed numerically harmonic enrichment functions for handling complex crack problems, which is based on the solution of the relevant Laplace equation.

The above propositions are, to a large extent, capable of describing either strong or weak discontinuities. However, to use enrichment functions correctly and efficiently in the XFEM, there is still a lack of systematic studies in view of the following questions remaining to be answered:

• How to justify or derive mathematically the existing functions, do there exist other enrichment functions, and which is better? • Is there some relationship between the two classes of enrichment functions?

• Given an enrichment function, what are the physical meaning of basic nodal variables associated with the continuous interpolation and that of complementary nodal variables accompanying with enrichment functions?

The third question usually arises with the introduction of enrichment functions, that is due to the fact that two variables will be associated with an enriched node, as will be discussed in Section 3.4. Also keep in mind that basic and complementary nodal variables involved respectively in the classical and enrichment interpolations are closely related to imposing Dirichlet boundary conditions, especially when discontinuities intersect with external boundaries [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF] [START_REF] Yvonnet | A general and efficient computational procedure for modelling the Kapitza thermal resistance based on XFEM[END_REF]. Till now, no direct attentions have yet been paid to this numerical issue.

In this paper, similarly to the first steps followed in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] [27] [START_REF] Rabczuk | A new crack tip element for the phantom-node method with arbitrary cohesive cracks[END_REF], we start with geometrical description of a triangular linear element where a strong discontinuity is present. Decompositions of both geometrical domain and displacement field in this cut element are then carried out (see Section 2). After this starting-point, the Hansbo and Hansbo's work [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] and the phantom-node-based method [START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF] [START_REF] Rabczuk | A new crack tip element for the phantom-node method with arbitrary cohesive cracks[END_REF] focused on description of the displacement field by making direct use of the Heaviside function. Finally, the former yielded two onesided independent interpolations while the latter formulated the displacement field as the sum of two element displacement fields, each containing two parts. In the present work, by extending linearly the subfields generated and separated by the interface to the whole element domain, we chose to define in Section 3 complementary nodal variables and construct some characteristic functions, aiming at a complete description of the discontinuous problem under investigation. A consistent two-step derivation procedure is then elaborated. In that process, four sets of basic nodal variables are found. For each set, the displacement field in the piecewise linear but discontinuous element is reproduced as the sum of two or three piecewise linear interpolations. One enrichment function turns up when comparing these interpolations with the standard XFEM form [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. Moreover, two general functions, one being node-dependent and the other node-independent, are obtained based on the four specific functions. Enrichment functions suitable for describing weak discontinuities are formulated in Section 4, which mainly relies on deeper understanding in terms of how enrichment functions for strong discontinuities are constructed. To achieve satisfactory computational efficiency and numerical accuracy, a hyperbolic tangent function serving as smoothing factor will be originally used instead of the level set function itself. The well-known numerical issue of blending elements will be addressed in Section 5. For validation purposes, we will perform in Section 6 numerical applications of the derived functions to heterogeneous problems with imperfect interfaces. Global convergence rates and distributions of stress and displacement will be presented.

Throughout the paper, we will use interfaces to refer to geometrical entities across which discontinuities take place. Cracks and perfect bonding between bulk phases then correspond to two extreme cases: the interface moduli are null for the former (cracks) but infinite for the latter.

PROBLEM SETTINGS

Let us consider a two-dimensional (2D) heterogeneous problem depicted in Figure 1. The geometrical domain Ω is separated by an interface S I into two parts Ω (1) and Ω (2) . Suppose that across S I the physical field defined in Ω suffers a strong or weak discontinuity. For numerical purposes, the domain Ω is discretized by a regular mesh of three-node triangular elements, into which the interface is then introduced. The mesh and the interface thus intersect each other, which may generate four types of elements, as shown in Figure 1.

Further assume that the interface S I can be defined mathematically by the zero level-set of a function φ : R 2 → R S

I (x) = {x ∈ R 2 | φ(x) = 0}. (1)
The value of φ at the i th node is determined as φ(x i ), x i being the position vector of the node.

In what follows, φ(x) and φ(x i ) will be symbolized respectively by φ x and φ i for brevity. It is convenient for later use to divide mesh nodes into two subsets: the subset of non enriched nodes and the subset of enriched nodes. The latter one, denoted by I * , includes the nodes of the elements whose supports intersect with the interface. The relevant numerical problem on blending elements will be addressed in Section 5. Nowadays, it is popular to deal with the above discontinuous problem within the context of the extended finite element method (XFEM). The basic idea behind the XFEM is to superpose one enrichment term to the classical continuous interpolation. To be definite, consider a piecewise linear displacement field u(x) defined in a cut triangular element Ω ijk (see Fig. 2a). According to [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the field u(x) is formulated in an additive form:

u(x) = ne i=1 N i (x)a i + me r=1 M r (x)ψ r (x)b r . (2)
Above, vectors a i and b i are, respectively, the basic and complementary nodal variables affiliated to an enriched node. Note that the number of nodal variables for an enriched node

4 ZHU QI-ZHI Ð (1)
Ð (1) S I S I cut elements blending elements (type 1)

blending elements (type 2) standard elements

Ð (2)
Ð (2) enriched nodes on the interface enriched nodes near the interface

S I S I Figure 1.
Illustration of different kinds of triangular linear elements which are generated when introducing the interface S I into the 2D regular mesh is doubled with respect to a non enriched one. N i (x) are shape functions, just like those used in the standard finite element method, which usually have to satisfy the requirement of the partition of unity, namely ne i=1 N i = 1. M r (x) represent shape functions related to the variables b r and most often take N r (x). n e is the number of all nodes of an element while m e ≤ n e is the number of enriched nodes of the same element. Thus, we can distinguish three situations:

• m e = n e for elements which are cut by the interface into two parts. All the nodes in such an element are enriched; • 0 < m e < n e for blending elements, only some of whose nodes are enriched; and

• m e = 0 for standard elements without enriched nodes.

Finally, ψ r (x) was coined in the literature as enrichment function, whose mathematical derivation is the main goal of the present work. The enrichment function is introduced to describe a prescribed type of discontinuities. An appropriate enrichment function with clear physical meanings plays a central role in XFEM-framed numerical modellings and simulations of heterogeneous problems with discontinuities. In that context, both numerical accuracy and computational convergence may depend on ψ (x) directly and significantly.

k k i i j j u i u i u j u j u k u k i i j j Ð + e Ð + e Ð ¡ e Ð ¡ e k k ® ® a) a) b) b) u + (x) u + (x) u ¡ (x) u ¡ (x) u (x) u (x) u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k u ¡ i u ¡ i u + j u + j c) c)

ENRICHMENT FUNCTIONS FOR STRONG DISCONTINUITIES

Our analyses will focus on linear elasticity problem. Nevertheless, application of the results to other discontinuous problems (e.g. thermal conductivity, piezoelectricity) will be direct. In literatures, a jump in displacement across an interface is commonly called strong discontinuity.

In contrast, a jump in strain across the interface, which generally leads to a jump in traction vector, is called weak discontinuity. This section addresses mathematical derivation of enrichment functions for strongly discontinuous problems. It is known beforehand that such functions must be discontinuous at the interface.

Domain decomposition and nodal variables

Referring to Figure 2, the triangular linear element Ω e is divided by the interface segment S αβ into two subdomains: Ω + e (x) where φ x > 0, and Ω - e (x) where φ x < 0 (Fig. 2a). Correspondingly, the displacement field u(x) is separated by S αβ into two subfields, i.e., u + (x) with x ∈ Ω + e and u -(x) with x ∈ Ω - e (see Fig. 2b). Now, we extend linearly both u + (x) and u -(x) towards the whole domain Ω e , as shown in Fig. 2c. The nodal values in the extended field u + (x) are denoted by u + i , u + j and u + k , and those in u -(x) by u - i , u - j and u - k . For later use, we also define the variables δu i , δu j and δu k each by a difference

δu i = u + i -u - i , δu j = u + j -u - j , δu k = u + k -u - k . (3)
These nine nodal variables will be used to reproduce the piecewise linear but discontinuous displacement field u(x) through shape functions and characteristic functions. In that process we will derive our enrichment functions ψ(x).

Characteristic functions

We first elaborate characteristic functions for describing the geometry and physical field in the cut element Ω e . Like the displacement, the field defined by the value of the function φ (x) can also be approximated by a linear interpolation, denoted by φ (x) or φx , such that φ (x) = ne i=1 N i (x)φ(x i ). The subdomains Ω + e (x) and Ω - e (x) then correspond to φx > 0 and φx < 0, respectively, while the interface, initially described rigorously by the level set function φ (x), is now approximately positioned by the function φx = 0. Therefore, the subdomains and the interface, to which point x may belong, can be fully identified by the signed value of the function φx as follows

sign( φx ) =      1, x ∈ Ω + e ; 0, x ∈ S I ; -1, x ∈ Ω - e . (4)
To complete the description, functions permitting to distinguish the subfields u + (x) and u -(x) are still required. For consistency's sake, the following two sign( φx )-based step functions are proposed:

X + (x) = 1 2 1 + sign( φx ) , (5)
X -(x) = 1 2 1 -sign( φx ) . (6)
As expected, functions ( 5) and ( 6) allow the interface to be stripped from bulk phases, because at the interface

X + (x ∈ S I ) = X + (x / ∈ S I ), X -(x ∈ S I ) = X -(x / ∈ S I ) but X + (x ∈ S I ) = X -(x ∈ S I ).
Moreover, for any point x, the above functions provide the partitions of the sign function and the unity,

X + ( φx ) -X -( φx ) = sign( φx ), (7)
X + ( φx ) + X -( φx ) = 1. ( 8
)
The last feature allows a two-part decomposition of physical field defined in the domain. For example, we have for the displacement field u(x):

u(x) = X + ( φx )u(x) + X -( φx )u(x) = u + (x) + u -(x). ( 9
)
ZHU QI-ZHI

In particular, when point x is on the interface, one has X + ( φx ) = X -( φx ) = 1/2. It is then possible to measure the displacement field at S I by the mean of u( φx = 0 + ) and u( φx = 0 -), such that

u( φx = 0) = 1 2 u( φx = 0 + ) + u( φx = 0 -) . (10)
Further, by introducing the generalized version of the Heaviside step function

H α (x) =    1, x > 0, α, x = 0, 0, x < 0, (11)
Equations ( 5) and ( 6) are respectively reformulated as

X + (x) = H α ( φx ), (12)
X -(x) = 1 -H α ( φx ). ( 13
)
The above two functions are illustrated in Figure 3. Remark that functions ( 5) and ( 6) are refound by taking α = 1 2 and that the basic functions used by Hansbo and Hansbo [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] are obtained by setting α = 1. Also note that the former are of Sign type while the latter are of Heaviside type. Therefore, the general forms ( 12) and ( 13) can be applied to derive in a same framework enrichment functions of both the Sign-and Heaviside types. To simplify the statement, we choose to present in detail our derivations and discussions for the Sign-type enrichment functions, but only provide in Appendix the results of the Heaviside type.

X + (x) X + (x) X ¡ (x) X ¡ (x) x x x x ® ® 1 ¡ ® 1 ¡ ® 0 1 0 1 Figure 3.
Illustration of the general characteristic functions X + (x) and X -(x). Solid circles define the function values at the interface. The values α = 1/2 and α = 1 yield the Sign-type functions ( 5) and ( 6) and the Heaviside functions, respectively.

Main steps towards the derivation of enrichment functions

Mathematical derivations are mainly based on domain decomposition and field reproduction. Note that in the standard form (2) the displacement field u(x) has been formulated as the sum of a continuous term and a discontinuous term. Thus, it is natural to reproduce u(x) by the following two steps:

(i) The first step consists in choosing a triplet of basic element nodal variables a i among

u + i , u + j , u + k , u - i , u - j and u - k .
The purpose is to construct the continuous term in Equation [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. To this end, two principles have to be followed: (a) according to the interpolation principle, two variables defined at a same node (with the same subscript) must not appear simultaneously. Under such restriction, the eight possible combinations corresponding to the cut element shown in Figure 2 are:

{u - i , u + j , u - k }, {u + i , u - j , u + k }, {u - i , u - j , u - k }, {u + i , u + j , u + k }, {u - i , u + j , u + k }, {u - i , u - j , u + k }, {u + i , u - j , u - k }, {u + i , u + j , u - k }; (
b) displacement variables should be placed in a same set when their associated nodes situate on a side of S I . For the example depicted in Figure 2, u - i and u - k , both locating on the same subdomain Ω - e and generated by extending linearly the subfield u -(x), should appear in pairs and so do u + i and u + k . Under such restriction, only the first four sets in the above list preserve, i.e.

{u - i , u + j , u - k }, {u + i , u - j , u + k }, {u - i , u - j , u - k }, {u + i , u + j , u + k }. (14)
It is emphasized that the above operations are equivalent to the following more general procedure: according to the sign of the values φ r , the nodes of the element Ω e are geometrically separated by the interface into two groups, i.e. (i, k) and (j), and correspondingly, the nodal displacements {u + i , u + j , u + k } in the extended field u + (x) are decomposed into two subsets, and so do {u - i , u - j , u - k }. We thus have four subsets of nodal variables:

{u + i , u + k }, {u + j }, {u - i , u - k } and {u - j }.
According to the principle of superposition, these subsets can and only can provide four combinations given in [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF]. It is worthy indicating that the above procedure is also valid for 2D and 3D meshes of higher order elements. (ii) The continuous displacement interpolation already formulated at step (i) will serve as a reference field. Thus, the second step is to establish the complementary discontinuous term so as that the sum of the two parts allows reproducing the initial piecewise linear displacement field which suffers a jump across the interface.

Derivation of basic enrichment functions

After the foregoing preliminary work, we are now in a position to proceed to derive enrichment functions. For each set of basic nodal variables previously determined, domain decomposition and reproduction will be performed upon the initial displacement field. This process justifies the use of {δu i , δu j , δu k } as the set of complementary variables and allows deriving enrichment function ψ r (x) for each case.

3.4.1. Decomposition I (D-I) with a = u - i , u + j , u - k
In the case in question, the discontinuous displacement field is divided into three parts, as shown in Figure 4. The first one F I-1 is continuous over the whole domain and the resulting linear interpolation is function of the nodal displacements {u - i , u + j , u - k }. Both the second one F I-2 and the third one F I-3 are discontinuous across S αβ . It is easily proved that the field u(x) in Ω e is reproduced by superposition of these three fields such that

u(x) = {N i , N j , N k } u - i , u + j , u - k T + X + (x) {N i , N j , N k } {δu i , 0, δu k } T -X -(x) {N i , N j , N k } {0, δu j , 0} T . ( 15
)
Note that δu i and δu k take effect only on subdomain Ω + e while δu j on Ω - e . Further, making use of ( 5) and [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF] allows to reformulate the terms on the right-hand side as follows

{N i , N j , N k } u - i , u + j , u - k T = 3 i=1 N i (x)a i , (16)
X + (x) {N i , N j , N k } {δu i , 0, δu k } T = X + (x) 3 r=1 N r (x)X + (x r )b r , (17)
X -(x) {N i , N j , N k } {0, δu j , 0} T = X -(x) 3 r=1 N r (x)X -(x r )b r , (18)
where a = u - i , u + j , u - k and b = {δu i , δu j , δu k } are the set of basic nodal variables and the set of complementary nodal variables, respectively. By substituting ( 16)-( 18), [START_REF] Sukumar | Modelling quasi-static crack growth with the extended finite element method Part I: computer implementation[END_REF] takes the form

u(x) = ne i=1 N i (x)a i (x) + ne r=1 N r (x) [X + (x)X + (x r ) -X -(x)X -(x r )] b r . (19)
Now, comparing [START_REF] Hettich | Interface material failure modeled by the extended finite-element method and level sets[END_REF] with the general form (2), it follows that

ψ 1 (x, x r ) = X + (x)X + (x r ) -X -(x)X -(x r ) = 1 2 sign( φx ) + sign(φ r ) . ( 20
)
ZHU QI-ZHI Equation ( 20) is the enrichment function corresponding to the set u - i , u + j , u - k and the domain decomposition shown in Figure 4.

u ¡ k u ¡ k k k i i j j u ¡ i u ¡ i u + j u + j ±u k ±u k k k i i j j ±u i ±u i k k i i j j F 0 F 0 F I¡1 F I¡1 F I¡2 F I¡2 F I¡3 F I¡3 u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k u ¡ i u ¡ i u + j u + j ¡±u j ¡±u j Ð + e Ð + e Ð ¡ e Ð ¡ e
+ i , u - j , u + k }. Accordingly, the displacement field u(x) is given by u(x) = {N i , N j , N k } u + i , u - j , u + k T + X + (x) {N i , N j , N k } {0, δu j , 0} T -X -(x) {N i , N j , N k } {δu i , 0, δu k } T . ( 21
)
Remark that the (δu i , δu k )-related nodes and the subdomain X -(x) locate on the opposite sides with respect to S I , and so do δu j and X + (x). By taking a = u + i , u - j , u + k and b = {δu i , δu j , δu k } and using the same procedure as that for decomposition (D-I), the second enrichment function is found:

ψ 2 (x, x r ) = X + (x)X -(x r ) -X -(x)X + (x r ) = 1 2 sign( φx ) -sign(φ r ) . (22)

Decomposition III (D-III)

Referring to Figure 6, the continuous part is obtained by extending the subfield u + (x) to the whole domain Ω e and the discontinuous interpolation involving complementary variables is employed only to subdomain Ω - e . In this case, the linear displacement field u(x) in Ω e is reproduced as

u(x) = {N i , N j , N k } u + i , u + j , u + k T -X -(x) {N i , N j , N k } {δu i , δu j , δu k } T , (23)
and the following enrichment function is derived

ψ 3 (x) = -X -(x) = 1 2 sign( φx ) -1 . (24)
Accordingly, the sets of basic and complementary nodal variables are

u + i , u + j , u + k and {δu i , δu j , δu k }, respectively. k k i i j j F II¡1 F II¡1 F 0 F 0 u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k u ¡ i u ¡ i u + j u + j u + k u + k u + i u + i u ¡ j u ¡ j ¡±u k ¡±u k ¡±u i ¡±u i ±u j ±u j k k i i j j F II¡2 F II¡2 k k i i j j F II¡3 F II¡3 Ð + e Ð + e Ð ¡ e Ð ¡ e
u(x) = {N i , N j , N k } u - i , u - j , u - k T + X + (x) {N i , N j , N k } {δu i , δu j , δu k } T , (25)
which leads to the fourth enrichment function,

ψ 4 (x) = X + (x) = 1 2 sign( φx ) + 1 . ( 26
)
The sets of basic and complementary variables associated with ( 26) are u - i , u - j , u - k and {δu i , δu j , δu k }, respectively. Verification for the case of nodes on the interface Although it's not always the case, we may meet particular configurations of cut elements where one or several nodes locate on the interface. Take as an example the element shown in Fig. 8a where, by assigning, the node x i locates on S I . From geometrical point of view, this node belongs not only to subdomain Ω + e (x) but also to subdomain Ω - e (x). Therefore, there exist two (D-I ) domain decompositions,

k k i i j j u + j u + j F III¡1 F III¡1 F 0 F 0 u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k Ð + e Ð + e Ð ¡ e Ð ¡ e u ¡ i u ¡ i u + j u + j u + k u + k u + i u + i ¡±u k ¡±u k ¡±u j ¡±u j i i j j k k F III¡2 F III¡2 ¡±u i ¡±u i
ZHU QI-ZHI u ¡ k u ¡ k k k i i j j u ¡ j u ¡ j u ¡ i u ¡ i k k i i j j F IV ¡1 F IV ¡1 F IV ¡2 F IV ¡2 ±u i ±u i F 0 F 0 u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k u ¡ i u ¡ i u + j u + j ±u k ±u k ±u j ±u j Ð + e Ð + e Ð ¡ e Ð ¡ e Figure 7
. Two-part domain decomposition of the displacement field with a = u - i , u - j , u - k as illustrated in Fig. 8b and8c. Accordingly, the displacement field can be reproduced by

u 1 (x) = {N i , N j , N k } u - i , u + j , u - k T + X + (x) {N i , N j , N k } {δu i , 0, δu k } T -X -(x) {N i , N j , N k } {0, δu j , 0} T , ( 27
)
for the domain decomposition shown in Fig. 8b, or by

u 2 (x) = {N i , N j , N k } u + i , u + j , u - k T + X + (x) {N i , N j , N k } {0, 0, δu k } T -X -(x) {N i , N j , N k } {δu i , δu j , 0} T . ( 28
)
for the domain decomposition shown in Fig. 8c. According to the relation [START_REF] Menk | Crack growth calculations in solider joints based on microstructural phenomena with X-FEM[END_REF], u(x) is defined as the mean of u 1 (x) and u 2 (x)

u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j Ð + e Ð + e Ð ¡ e Ð ¡ e u + k u + k u + j u + j u ¡ i u ¡ i S®S®¯u ¡ k u ¡ k k k i i j j u + j u + j u ¡ i u ¡ i k k i i j j ±u k ±u k ±u i ±u i ¡±u j ¡±u j k k i i j j u ¡ k u ¡ k k k i i j j u + j u + j k k i i j j ±u k ±u k ¡±u i ¡±u i ¡±u j ¡±u j k k i i j j u + i u + i (a) (a) (b) (b) (c) (c)
u(x) = 1 2 [u 1 (x) + u 2 (x)] . (29)
Then insertion of (??) and ( 28) into [START_REF] Olivier | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals[END_REF] gives

u(x) = {N i , N j , N k } (u + i + u - i )/2, u + j , u - k T + X + (x) {N i , N j , N k } {δu i /2, 0, δu k } T -X -{N i , N j , N k } {δu i /2, δu j , 0} T . ( 30
)
Note that at node x i , one has

X + (x i ) = X -(x i ) = 1/2. Hence, by setting a = (u + i + u - i )/2, u + j , u - k
as the set of basic nodal variables and with the set of complementary variables b = {δu i , δu j , δu k } unchanged, Equations ( 17) and ( 18) hold for [START_REF] Merle | Solving thermal and phase change problems with the extended finite element method[END_REF] and the same enrichment function ( 20) can be obtained.

The above analyses are completely applicable to the (D-II ) domain decomposition, and the same conclusion can be drawn. However, no additional checking is needed for functions [START_REF] Pedro | A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF] and ( 26) because of their node-independent nature.

Formulation of general enrichment functions

Based on the four basic enrichment functions we obtained, two general functions, one being node-dependent and the other node-independent, can be constructed. As a particular case, a fifth specific enrichment function will be reported.

It follows by substituting the node-dependent functions ψ 1 (x, x r ) and ψ 2 (x, x r ) into the general form ( 2)

u(x) = ne i=1 N i (x)a 1 i + me r=1 N r (x)ψ 1 (x, x r )b r = ne i=1 N i (x)a 2 i + me r=1 N r (x)ψ 2 (x, x r )b r , (31)
where

a 1 = u - i , u + j , u - k and a 2 = u + i , u - j , u + k
are the sets of basic nodal variables associated with ψ 1 (x, x r ) and ψ 2 (x, x r ), respectively. We then construct, through one parameter λ ranging from 0 to 1, a linear combination of the terms after the first and second equal sign in [START_REF] Ji | On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method[END_REF], such that

u(x) = ne i=1 N i (x) (1 -λ)a 1 i + λa 2 i + me r=1 N r (x) (1 -λ)ψ 1 (x, x r ) + λψ 2 (x, x r ) b r , (32)
from which the general node-dependent enrichment function is derived

ψ(x, x r ) = (1 -λ)ψ 1 (x, x r ) + λψ 2 (x, x r ) = 1 2 sign( φx ) + (1 -2λ)sign(φ r ) . (33)
Correspondingly, the set of basic nodal variables is also a linear combination:

a i (x) = (1 -λ)a 1 i + λa 2 i .
It can be easily proved that Equation (33) reduces to ψ 1 (x, x r ) for λ = 0 and to ψ 2 (x, x r ) for λ = 1. In particular, λ = 1 2 leads to the fifth specific function

ψ 5 (x) = 1 2 sign( φx ), (34)
which is accompanied by the set of basic nodal variables

a = 1 2 a 1 + a 2 = u + i + u - i 2 , u + j + u - j 2 , u + k + u - k 2 , (35)
and the set of complementary nodal variables

b = {δu i , δu j , δu k } = u + i -u - i , u + j -u - j , u + k -u - k . ( 36
)
12

ZHU QI-ZHI
Remark that although the basis functions ψ 1 (x, x r ) and ψ 2 (x, x r ) are node-dependent, ψ 5 (x) is node-independent.

In the same way, when applying the two node-independent function ψ 3 (x) and ψ 4 (x), we obtain the following interpolation

u(x) = ne i=1 N i (x) (1 -λ)a 3 i + λa 4 i + me r=1 N r (x) (1 -λ)ψ 3 (x) + λψ 4 (x) b r , (37)
with

a 3 = u + i , u + j , u + k and a 4 = u - i , u - j , u - k , which
gives the following general nodeindependent enrichment function

ψ(x) = (1 -λ)ψ 3 (x) + λψ 4 (x) = 1 2 sign( φx ) -(1 -2λ) . (38)
Again, λ = 1 2 yields function [START_REF] Tarancón | Enhanced blending elements for XFEM applied to linear elastic fracture mechanics[END_REF] with the same sets of basic nodal variables a and complementary nodal variables b as those already given in [START_REF] Ventura | Fast integration and weight function blending in the extended finite element method[END_REF] and [START_REF] Bordas | Enriched finite element short course: class notes[END_REF], respectively.

Validation and extension

The foregoing mathematical derivations have focused on a 2D linear but discontinuous field which is defined on a triangular linear element cut by one interface. It will be our purpose in this part to extend applications of the derived functions.

Elements cut by multiple interfaces

In the previous parts, only one unique interface (discontinuity) is considered. Now we deal with the cases where an element intersects with multiple interfaces. For simplicity, the interfaces intersecting with the element are assumed to be geometrically separated from each other. More precisely, there exist no intersections between the interfaces and no interface is included by another. Figure 9 illustrates a threecutting configuration and the associated domain decomposition (Fig. 9b). Suppose that the three interfaces can be fully described by the level sets φ i

x , φ j x and φ k x , respectively. In this example, the domain Ω e consists of four parts, i.e. Ω + e0 , Ω - ei , Ω - ej and Ω - ek and so does the displacement field, marked as u + 0 (x), u - i (x), u - j (x) and u - k (x). As performed in Section 3, the displacement subfields are then extended to the whole element domain. In the case shown in Figure 9, the nodal variables include u + 0r (r = i, j, k) representing the displacement vectors at nodes x r which are determined in the extended field u + 0 (x), and u - sr denoting the node x rrelated displacement vector in the extended field u - s (x). It is also noted that for the cutting pattern in question the number of variables for each node becomes four. By making direct use of the enrichment functions already obtained in the one-cutting case, the piecewise linear but discontinuous displacement field in Ω e is expressed by:

i j k i j k Ð + e Á j Á i Á k + + + ¡ ¡ ¡ Ð + e0 Ð ¡ ei Ð ¡ ek Ð ¡ ej
u(x) = ns r=1 X r + (x) ne i=1 N i (x)u + 0i + ns r=1 X r -(x) ne i=1 N i (x)a r i + ne i=1 N i (x)ψ r i (x) b r
or also formulated equivalently as the sum of continuous interpolations and discontinuous interpolations

u(x) = ns r=1 X r + (x) ne i=1 N i (x)u + 0i + ns r=1 X r -(x) ne i=1 N i (x)a r i + ns r=1 X r -(x) ne i=1 N i (x)ψ r i (x) b r i , (40)
with the interface-dependent characteristic functions

X r + (x) = 1 2 1 + sign( φr x ) , X r -(x) = 1 2 1 -sign( φr x ) , (41)
where n s is the number of interfaces which intersect with the edges of an element and here we take n s = 3. The enrichment functions are now interface-dependent and also probably nodedependent. In the one-cutting case, their expressions have been given in ( 20), ( 22), ( 24), ( 26) and [START_REF] Tarancón | Enhanced blending elements for XFEM applied to linear elastic fracture mechanics[END_REF]. Both the nodal variables a r i (x) and b r i (x) are now interface-related while the former also depends on the choice of enrichment functions.

We are particularly interested in the enrichment function

ψ r 3 (x) = 1 2 sign( φr x ) -1 , (42)
which corresponds to the triplet of basic nodal variables a = u + 0i , u + 0j , u + 0k . With the help of the property

ns r=1 X r + (x) + ns r=1 X r -(x) = 1, (43)
the displacement formulation ( 40) is then simplified as the sum of two terms

u(x) = ne i=1 N i (x)a i + ns r=1 X r -(x) ne i=1 N i (x)ψ r 3 (x) b r i (44) with b r i = u + 0i -u - ri and a i = u + 0i .
It is stressed that such two-term compact form (44) is not obtainable for either of enrichment functions [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF], [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF], [START_REF] Yvonnet | A general and efficient computational procedure for modelling the Kapitza thermal resistance based on XFEM[END_REF] or [START_REF] Tarancón | Enhanced blending elements for XFEM applied to linear elastic fracture mechanics[END_REF]. In this sense, function ψ r 3 (x) would be of particular interest for modelling strongly discontinuous problems with complex geometries.

Remark: A triple junction can be regarded as the asymptotic case of a three-cutting configuration, like that depicted in Figure 9. For polycristals with discontinuous grain boundaries, Simone and his co-workers [START_REF] Simone | A generalized finite element method for polycrystals with discontinuous grain boundaries[END_REF] and Menk and Bordas [START_REF] Menk | Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals[END_REF] made use of discrete 0-1-valued characteristic function (named as indicator function in [START_REF] Menk | Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals[END_REF]) to identify grains which constitute the structure. The enrichment function ψ r 3 we derived plays the same role (Inside the zone of a grain, the associated function takes the value -1, and 0 otherwise) and therefore we can adapt the interpolations (44) for handling discontinuous problems with triple junctions.

3.6.2. Three-dimensional linear element cut by one interface Our concern here is with the most popular mesh of four-node tetrahedral linear elements. Figure 10 shows four possible configurations of a tetrahedron cut by one interface into two parts (Ω + e and Ω - e ). Intersection points between the interface and the element supports (edges) can be easily detected (using the level-set method).

The displacement vectors at nodes i, j, k, l are denoted as u i , u j , u k and u l , respectively. The displacement field u(x) is separated by the interface into two parts: u + (x) with x ∈ Ω + e and u -(x) with x ∈ Ω - e , both linear in their respective subdomains. Similarly to the 2D case, we extend linearly u + (x) and u -(x) to the whole element domain. By following the ZHU QI-ZHI procedure elaborated in Section 3.3, for each cutting configuration shown in Figure 10, we can generate four sets of basic nodal variables. Take the case in Figure 10b as an example, the four sets of basic nodal variables are {u

- i , u + j , u + k , u + l }, {u + i , u - j , u - k , u - l }, {u + i , u + j , u + k , u + l } and {u - i , u - j , u - k , u - l }.
To check the applicability of the derived functions for three dimensional analyses, let us turn back to the 2D formulations. By taking advantage of the property (8), the continuous term ne i=1 N i (x)a i (x), which is formally common to (15), ( 21), ( 23) and [START_REF] Mousavi1 | Higher-order extended finite elements with harmonic enrichment functions for complex crack problems[END_REF], is cast into the form

ne i=1 N i (x)a i (x) = ne i=1 X + (x)N i (x)a i + ne i=1 X -(x)N i (x)a i . (45)
Then, the formulae ( 15), ( 21), ( 23) and ( 25) are all rearranged as

u(x) = ne i=1 X + (x)N i (x)u + i + ne r=1 X -(x)N i (x)u - i . ( 46
)
Obviously, the displacement field (46) in a triangular linear element reads as the sum of two independent one-sided linear interpolations, which are associated with the subdomains Ω + e (x) and Ω - e (x), respectively. Extension of (46) to the case of 3D linear elements is natural and direct. After this extension, in order to derive enrichment functions, only an inverse operation is needed to make occurrence of a continuous interpolation and one or two discontinuous interpolations (referring to Section 3.4). Finally, the same general and specific enrichment functions as those presented in the 2D linear case can be obtained. In other words, the enrichment functions we derived for the 2D case are completely applicable to 3D piecewise linear meshes. 3.6.3. Higher order elements In order to verify the applicability of the derived functions to higher order elements, let us summarize the two-step derivation procedure already elaborated for linear elements: firstly, according to the value of the level set φ x , both the geometry and physical field are decomposed into two parts and correspondingly element nodes are divided into two subsets, one related to subfield Ω + e where φ x > 0 and the other to Ω - e with φ x < 0; secondly, for each node of a cut element, two variables are generated during linear extension of the subfields. The element nodes, geometrically separated by the interface, can and only can provide four combinations of basic nodal variables. For each combination, one enrichment function has been derived. Note that in this process, except the extension of the two subdomains, all other steps are independent on the order of a mesh which is reflected by shape functions. Therefore, when no multiple cuttings by a same interface take place, extension of the subfields of a higher order to the whole element domain is also performable and the above derivation procedure is still valid.

Basic features of the obtained enrichment functions

The enrichment functions we derived show the following basic features: • ψ 1 (x, x r ) and ψ 2 (x, x r ) are node-dependent while ψ 3 (x), ψ 4 (x) and ψ 5 (x) are nodeindependent. The node-related term sign(φ r ) in both ψ 1 (x, x r ) and ψ 2 (x, x r ) allows describing relative position between integration point x and node x r . The enrichment functions ψ 1 (x, x r ) and ψ 2 (x, x r ) are illustrated in Figures 11 and12, respectively, for which use has been made of the cut element Ω ijk shown in Figure 2. • At node x r , ψ 1 (x r , x r ) = sign(φ r ), which is generally non zero and results in the vector of nodal displacement: u(x r ) = a(x r ) + sign(φ r )b(x r ). Therefore, for blending elements, complementary integrations related to the enrichment term will be indispensable (see Figure 14a). This feature is shared by function ψ 5 (x). • In contrast with ψ 1 (x, x r ), one has ψ 2 (x r , x r ) = 0, leading to u(x r ) = a(x r ). Thanks to this feature, for blending elements of type 2, numerical integrations upon the enrichment term will be removed automatically (see Figure 14b). • Both functions ψ 3 (x) and ψ 4 (x) provide one-sided enrichment. Precisely, with respect to the interface, ψ 3 (x) only enriches the inner subdomain Ω - e (x) where sign( φx ) ≤ 0 (see Figure 14c), while ψ 4 (x) only enriches the outer subdomain Ω + e (x) where sign( φx ) ≥ 0 (see Figure 14d). Accordingly, the issue on blending elements is partially removed for these two cases. We will call inner enrichment for ψ 3 (x) and outer enrichment for ψ 4 (x). The node-independent enrichment functions ψ 3 (x), ψ 4 (x) and ψ 5 (x) are depicted in Figure 13. It is noted that functions ψ i (x), i = 1, 2, 3, 4 are 0-1 valued in bulk phases but ψ 5 (x) is not the case.

k k i i j j Ð (+) e Ð (+) e Ð (¡) e Ð (¡) e k k i i j j ® ® ¯® ® à 1 j = 0 à 1 j = 0 à 1 i = à 1 k = 1 à 1 i = à 1 k = 1 à 1 i = à 1 k = 0 à 1 i = à 1 k = 0 à 1 j = ¡1 à 1 j = ¡1
k k i i j j ® ® k k i i j j ® ® Ã 2 i = Ã 2 k = 0 Ã 2 i = Ã 2 k = 0 Ã 2 i = Ã 2 k = ¡1 Ã 2 i = Ã 2 k = ¡1 Ã 2 j = 1 Ã 2 j = 1 Ã 2 j = 0 Ã 2 j = 0
• For the specific functions, their complementary integration zones over the regular mesh shown in Figure 1 are illustrated and compared in Fig. 14. • For the general and specific enrichment functions we obtained, the set of basic nodal variables a r differs from each other while the set of complementary nodal unknowns b r is identical, which results in one unique form of the displacement jump vector

u (x) = u( φx = 0 + ) -u( φx = 0 -) = ne r=1 N r (x)b r . ( 47
)
Obviously, there is no enrichment function involved in (47), just like in the standard finite element formulation. Moreover, at node x r , u (x r ) = b r , implying that the Figure 14. Complementary integration (shaded) zones over the mesh shown in Fig. 1, which include a) all blending elements for function [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF], b) no blending elements for function [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF], c) inner enrichment with partial blending elements for function [START_REF] Pedro | A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF], and d) outer enrichment with partial blending elements for function [START_REF] Yvonnet | A general and efficient computational procedure for modelling the Kapitza thermal resistance based on XFEM[END_REF].

ZHU QI-ZHI k k i i j j ® ® Ã 3 = 0 Ã 3 = 0 k k i i j j ® ® Ã4 = 1 Ã 4 = 1 Ã 3 = ¡1 Ã 3 = ¡1 Ã 4 = 0 Ã 4 = 0 k k i i j j ® ® Ã5 = ¡ 1 2 Ã 5 = ¡ 1 2 Ã 5 = 1 2 Ã 5 = 1 2
complementary variables b r can be physically interpreted as the vector of displacement jump at enriched nodes x r . Based on these observations, the standard formula (2) can then be reformulated by replacing the variables b r by u r , such that

u(x) = ne i=1 N i (x)a i + me r=1 M r (x)ψ r (x) u r . (48)
Equation ( 48) coincides with the form proposed in [START_REF] Olivier | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals[END_REF] where ψ(x) has taken the classic Heaviside function and the first continuous term on the right-hand side of (48) was called the regular part of the displacement field.

Comparisons between enrichment functions

In this part, we address comparisons between the obtained functions and establish their relationships with the existing ones.

• The Heaviside-type enrichment functions in Appendix include the popular ones reported in the literature and also provide some new choices. Originally, for each function, the physical meanings of the basic and complementary nodal variables are revealed. This work clearly confirms both the classic Heaviside step function itself and its shifted form can be used for XFEM-framed modelling of strong discontinuities.

• The function ψ 1 (x, x r ) takes effect only when point x and node x r takes a same sign, implying that the enrichment term is the sum of two one-sided independent parts, which is just the basic idea behind the work of Hansbo and Hansbo [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]. However, differently from the latter, our formulation within the XFEM context also contains a continuous term (see Equation ( 19)). Moreover, according to the comment by Areias Pedro and Belytschko [START_REF] Pedro | A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF], Hansbo and Hansbo's formulations [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] are equivalent to use two Heaviside-type step functions. The present work shows that from numerical viewpoint, these two functions can be combined into one single function ψ 1 (x, x r ). The Hansbo-type interpolations in terms of the characteristic functions X + (x) and X -(x) have been given in (46). The corresponding domain decomposition is depicted in Figure 15.

k k i i j j u + j u + j u ¡ k u ¡ k k k i i j j u ¡ i u ¡ i u ¡ k u ¡ k k k i i j j u + i u + i u ¡ j u ¡ j u + k u + k Ð (+) e Ð (+) e Ð (¡) e Ð (¡) e u ¡ i u ¡ i u + j u + j u + k u + k u + i u + i u ¡ j u ¡ j Figure 15.
A two-sided domain decomposition according to [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] • The function ψ 3 (x) provides an anti-Heaviside enrichment. It is also worth recalling that for multiple cutting configurations this inner enrichment function provides a more compact form and therefore shows more interests for numerical implementation. • The function ψ 4 (x) corresponds to the modified Heaviside function

ψ 4 (x) = 1 2 sign( φx ) + 1 = H 1 2 (x) (49)
where the subscript 1 2 defines the value at the origin, i.e. H(x = 0) = 1 2 , which is exactly the Heaviside function used in MATLAB .

• The functions ψ 3 (x), ψ 4 (x) and ψ 5 (x) are linked to ψ 2 (x, x r ) through their respective shifted forms, i.e.

ψ 3 (x) -ψ 3 (x r ) = ψ 4 (x) -ψ 4 (x r ) = ψ 5 (x) -ψ 5 (x r ) = ψ 2 (x,x r ) . (50)
In addition, the mean of ψ 3 (x) and ψ 4 (x) leads to the function ψ 5 (x) and so does that of ψ 1 (x, x r ) and ψ 2 (x, x r ). • The derivation of the function ψ 5 (x) not only justifies the use of the sign function itself in describing strong discontinuities but also enriches the existing work by clarifying the physical meanings of the basic and complementary nodal unknowns. Note that ψ 2 (x, x r ) is the shifted form of ψ 5 (x) and that these two functions have often been used in XFEM formulations. • Finally, the factor 1 2 , which appears in all the obtained enrichment functions, allows avoiding the occurrence of possible multipliers in the formulae u and b, which now take the terse forms

u (x) = ne r=1 N r (x)b r , b r = u + r -u - r . (51)
From numerical point of view, it is, of course, possible to shift this factor of 1 2 to the term b (see e.g. [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF]), which will add an factor of 2 to the interpolation of u , however.

ZHU QI-ZHI

In summary, the widely used XFEM enrichment functions for strong discontinuities can be derived in the current framework. To the author's acknowledge, the general enrichment functions [START_REF] Chessa | On the construction of blending elements for local partition of unity enriched finite elements[END_REF] and [START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modelling of thin interphases as imperfect interfaces[END_REF] as well as the specific forms (20), ( 24) and ( 26) have never been reported.

Imposing boundary conditions

In the XFEM context, Dirichlet continuous boundary conditions for standard elements are imposed directly, just like that in standard finite element analyses. However, some difficulties may arise when discontinuities intersect with external mesh boundaries [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF] [START_REF] Yvonnet | A general and efficient computational procedure for modelling the Kapitza thermal resistance based on XFEM[END_REF]. Several measures have been taken to handle this numerical issue [START_REF] Merle | Solving thermal and phase change problems with the extended finite element method[END_REF] [31] [START_REF] Béchet | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF]. For instance, Lagrange multipliers have been employed to relax the strong imposition of such conditions. The present work helps to impose external boundary conditions directly and efficiently.

Note that for each enrichment function, the associated physical meanings of both the basic and complementary variables have been clarified. The one-to-one relationships between nodal variables and enrichment functions are provided in Table I. To impose Dirichlet boundary conditions for nodes enriched and located on external boundaries, we only need to set the values of the basic and complementary variables a and b according to the choice of enrichment function.

Table I. Imposing Dirichlet boundary conditions for an enriched node

x i (ū i = (u + i + u - i )/2) Case ψ 1 (x, x r ) ψ 2 (x, x r ) ψ 3 (x) ψ 4 (x) ψ 5 (x) all 5 eqns. x i ∈ I * φ i > 0 a i = u - i a i = u + i a i = u + i a i = u - i a i = ūi b i = u + i -u - i ,x i / ∈ S I φ i < 0 a i = u + i a i = u - i a i = u + i a i = u - i a i = ūi b i = u + i -u - i x i ∈ S I φ i = 0 a i = ūi a i = ūi a i = u + i a i = u - i a i = ūi b i = u + i -u - i

ENRICHMENT FUNCTIONS FOR WEAK DISCONTINUITIES

The previous section has addressed the derivation of enrichment functions for strong discontinuities. In practical applications especially for nanocomposites and nanostructures, there also exists another type of heterogeneous problems where weak discontinuities may be present, for instance, problems with coherent imperfect interfaces, in which the displacement vector remains continuous across an interface while the traction vector suffers a jump across the same interface. This section is dedicated to establishing enrichment functions for describing and modelling such weak discontinuities. The way in which we achieve the objective is somewhat subtle and mainly based on deeper understanding in terms of how the enrichment functions we obtained for strong discontinuities are constructed in a unified way.

Deeper insights into the derived functions

Let us put together the basic enrichment functions ψ i which have been derived in Section 3 for strong discontinuities

ψ 1 ( φx , φ r ) = 1 2 sign( φx ) + sign(φ r ) , (52a)
ψ 2 ( φx , φ r ) = 1 2 sign( φx ) -sign(φ r ) , (52b)
ψ 3 ( φx ) = 1 2 sign( φx ) -1 , (52c)
ψ 4 ( φx ) = 1 2 sign( φx ) + 1 . ( 52d
)
Recall that no enrichment functions but only the jump of displacement vector intervenes in numerical integrations related to interfaces [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF] [26]. Hence, in numerical implementations, the discontinuous enrichment functions ψ i are involved only in weak formulations for bulk phases.

In practice, integration points in a cut element are usually determined by two successive steps: firstly, according to its cutting configuration, the domain Ω e is divided into several subelements ∆Ω i e such that ∪∆Ω i e = Ω e and ∆Ω i e ∩ ∆Ω j e = ⊘ for i = j; secondly, for each subelement, integration points are found most often by using a Gaussian quadrature. Integration points thus determined are all inside Ω e , implying no integration points locate on interfaces or element edges. In summary, we have sign( φx ) = 0 from practical viewpoint, leading to sign( φx )sign( φx ) = 1. By applying this relation, functions (52a)-(52d) are rewritten in a unified form

ψ i ( φx , φ r ) = 1 2 sign( φx )P i ( φx , φ r ), i = 1, 2, 3, 4 (53)
with

P 1 ( φx , φ r ) = 1 + sign( φx )sign(φ r ); (54a) P 2 ( φx , φ r ) = 1 -sign( φx )sign(φ r ); (54b) P 3 ( φx ) = 1 -sign( φx ); (54c) P 4 ( φx ) = 1 + sign( φx ). ( 54d
)
Obviously, the above four functions ψ i differ from each other only by the terms P i ( φx , φ r ) which in fact characterize their respective complementary integration zones as have been remarked in Section 3.7. The combination of the first two terms, 1 2 sign( φx ), is exactly the fifth specific enrichment function ψ 5 which has been found in Section 3.5. The use of the factor 1 2 is argued in Section 3.8. For completeness, we can impose P 5 = 1 for ψ 5 . It is finally known from ψ 5 that the second term sign( φx ) in (53) plays the crucial role of describing strong discontinuities.

Enrichment functions for weakly discontinuous problems

We now proceed to construct an appropriate function characteristic to describe weak discontinuities, just like sign( φx ) for strong discontinuities. By its nature, such a function must be continuous at an interface while its derivative must be discontinuous across the same interface so as to trigger requisite weak discontinuities. To be consistent, candidates in terms of the function φx are usually attractable for combined XFEM/LSM modellings. Among those fulfilling the above requirements, the piecewise linear absolute function φx may be the simplest one. By replacing sign( φx ), the second term on the right-hand side of the functions in (53), by φx , we formulate enrichment functions for weak discontinuities as

ϕ i ( φx , φ r ) = 1 2 φx P i ( φx , φ r ), i = 1, 2, 3, 4, 5 (55)
It follows with the help of the property φx = φx sign( φx ) :

ϕ i ( φx , φ r ) = 1 2 φx sign( φx )P i ( φx , φ r ), i = 1, 2, 3, 4, 5 (56)
which are equivalent to the formulae

ϕ i ( φx , φ r ) = φx ψ i ( φx , φ r ), i = 1, 2, 3, 4, 5 (57)
Function (61) also allows us to avoid the aforementioned disadvantage residing in function (60) but it may cause the following numerical problem. In Figure 17, the element Ω ijk is divided by the interface S I into two subdomains Ω αiβ and Ω αβjk . The fictive perfect interface segment S kβ ′ is parallel to S I and passes the node x k . According to function (61), the value of ϕ(x, x k ) = φxφk will change its sign across S kβ ′ . It is, however, quite difficult to justify this change in sign from geometrical or physical viewpoints. In fact, it may be a disadvantage when taking into account the partition of unity. A simple study shows that functions (57) and ( 60) can ensure the uniformity of the function sign in the subdomains.

S I S I i i j j k k zone with ' k < 0 zone with ' k < 0 ® ® ¯¯0 ¯0 zone with ' k > 0 zone with ' k > 0

Improvements upon the functions

Functions (57) are actually constructed by multiplying (53) with φx . The continuous factor φx plays the role of smoothing the step function sign( φx ). Also note that functions (53) for strong discontinuities are bounded by -1 and 1 while (57) are not. To be consistent and inspired by [START_REF] Bordas | Enriched finite element short course: class notes[END_REF], the continuous factor φx is generalized and replaced by its hyperbolic tangent function, such that ϕ i ( φx , φ r ) = tanh(ρ φx )ψ i ( φx , φ r ), i = 1, 2, 3, 4, 5.

(62)

Above, the coefficient ρ is introduced to control the kinematics of the value of the function tanh(ρ φx ) towards the asymptotic value -1 or 1. When ρ = 1, tanh( φx ) adopts the following series of Taylor expansion tanh( φx ) = φx -φ3

x 3 + 2 φ5 x 15 - 17 φ7 x 315 + • • • (63)
which implies that the primary proposition φx is exactly the first-order approximation of tanh( φx ).

The coefficient ρ can be constant or evolve with mesh refinement. For the latter, we can use a characteristic value φ c to (partially) normalize the function φx in such a way that

ϕ i ( φx , φ r ) = tanh φx φ c ψ i ( φx , φ r ), (64)
which is equivalent to take in (62) ρ = 1/φ c . Comparisons with different values of ρ will be performed in Section 6.1.

BLENDING ELEMENTS

In the previous sections, the derived enrichment functions have been verified for elements whose supports (edges) are cut by interfaces and whose nodes are all enriched. When not all ZHU QI-ZHI nodes of an element are enriched, this element is termed as blending element. Such elements blend fully enriched elements (with all their nodes enriched) with standard elements (without nodes enriched) [START_REF] Chessa | On the construction of blending elements for local partition of unity enriched finite elements[END_REF]. Because the number of enriched nodes in a blending element m e is less than that of all nodes of an element n e , the requirement on the partition of the unity is no longer satisfied. Recently, Fries [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF] addressed some comments on this numerical topic and proposed a solution by reproducing the enrichment functions for all element types. However, the Fries' technique increases the order of enrichment term from linear to quadratic and needs to enrich all nodes in blending elements. Later, Tarancón et al. [START_REF] Tarancón | Enhanced blending elements for XFEM applied to linear elastic fracture mechanics[END_REF] proposed to enhance blending elements by increasing the polynomial order of approximation, which has provided greater accuracy in LEFM analyses. Ventura et al. [START_REF] Ventura | Fast integration and weight function blending in the extended finite element method[END_REF] proposed a method for alleviating the errors in the blending elements. The idea behind this method is to pre-multiply the enrichment function by a smooth weight function with compact support to allow for a completely smooth transition between enriched and non enriched subdomains. Here we propose a simple and efficient solution to this basic issue in the XFEM. Note that the present derivation framework allows the piecewise linear field in cut elements (including their nodal values) to be reproduced.

It is then natural to enlarge the application of the enrichment to all system nodes. For this purpose, we only need to use a discrete characteristic function κ(x i ) defined on all mesh nodes such that κ(x i ) = 1 when x i ∈ I * , and κ(x i ) = 0, otherwise. Thus, the displacement field ( 2) is updated as

u(x) = ne i=1 N i (x)a i (x) + ne i=1 N i (x)ξ(x, x i )b i , (65)
where ξ(x, x i ) takes κ(x i )ψ(x, x i ) for strong discontinuities and κ(x i )ϕ(x, x i ) for weak discontinuities. The above treatment shows significant advantages: i) for strong discontinuous problems, the discontinuous interpolation remains linear; ii) in numerical implementations, additional unknowns are needed only for enriched nodes where κ(x i ) = 1; and iii) it is of unified nature, that is to say, all the obtained enrichment functions work well. Let us consider a 3D elasticity problem. Numerically, the above treatment is realized as follows:

Denote n s as the number of system degrees of freedom and N k as the number of degrees of freedom cumulated from the 1 st node to the k th node. For an enriched node indexed in the system as I, the set of its nodal unknowns takes the form [a I1 , a I2 , a I3 , b I1 , b I2 , b I3 ] T . In the column matrix of system unknowns, these components occupy the positions N I-1 + 1, N I-1 + 2, N I-1 + 3, N I-1 + 4, N I-1 + 5 and N I-1 + 6, respectively. When node I is non enriched, i.e. κ(x I ) = 0, the last three components (related to nodal variables b) in the above list will be replaced by n s + 1, such that N I-1 + 1, N I-1 + 2, N I-1 + 3, n s + 1, n s + 1, and n s + 1. Accordingly, we should enlarge the dimension of the column matrix of system unknowns u from n s to n s + 1 and set in prior u(n s + 1) = 0. That is the small cost caused by the proposed treatment.

NUMERICAL EXAMPLES

In this section, we will validate the derived enrichment functions by applying them to heterogeneous problems with imperfect interfaces. For simplicity, numerical calculations are limited to the case of one curved interface. Global convergence rates are assessed via the benchmark of a spherical inclusion coated by coherent imperfect interface and embedded in an infinite matrix. We also present displacement field and stress distribution for a 2D discontinuous problem with one embedded spring-layer interface.

A spherical inclusion with coherent interface

In the recent work of Zhu and his co-workers [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF] concerning 3D numerical modelling via XFEM of spring-layer imperfect interfaces effects, quite satisfactory convergence rates have been reported for various combinations of material constants, in which use has been made of the enrichment function ψ = 1 2 sign( φx )(1 -sign( φx )sign(φ r )) for describing strong discontinuities and of Table I for imposing external boundary conditions. In this part, we will assess computational convergences for a weakly discontinuous benchmark with coherent imperfect interfaces.

Ð (2) Ð (2) Ð (1) Ð (1) S I S I R0 R0 O O Figure 18
. Illustration of a spherical inclusion Ω (2) coated by imperfect interface and embedded in an infinite matrix Ω (1) .

As depicted in Figure 18, the spherical inclusion Ω (2) enrobed by an imperfect coherent interface is embedded in an infinite medium Ω (1) and subjected to a remote uniform hydrostatic strain ε ∞ , such that ε ∞ 11 = ε ∞ 22 = ε ∞ 33 = ε 0 and null for other components. For this axisymmetric problem, the exact displacement field was given by [START_REF] Sharma | Effect of surfaces on the size-dependent elastic state of nano-inhomogeneities[END_REF]

u 0 (r) = Aε 0 r, 0 ≤ r ≤ R 0 ε 0 r + R 3 0 r 3 (A -1)ε 0 r, r ≥ R 0 (66)
with the coefficient A dependent on material constants and the radius of the inclusion

A = 3κ 1 + 4µ 1 3κ 2 + 4µ 1 + 4(λ s + µ s )/R 0 (67)
Above, κ 1 and µ 1 are the bulk and shear moduli of the matrix, respectively, and κ 2 is the bulk modulus for the inclusion. In isotropic case, λ s and µ s are Lamé constants characterizing the imperfect interface. Note that with the help of the Hadamard relation, the Taylor expansion as well as the approach of asymptotic analysis, the coherent interface model can be derived rigorously through an equivalent replacement of a rigid and thin interphase having the elasticity tensor C 0 by a mechanically imperfect interface [START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modelling of thin interphases as imperfect interfaces[END_REF]. In that context, when the interphase is uniform in thickness, denoted by h, and isotropic with two Lamé's constants λ 0 and µ 0 , the interface parameters λ s and µ s are determined by

λ s = 2µ 0 λ 0 + 2µ 0 hλ 0 , µ s = hµ 0 . ( 68
)
It follows that

λ s µ s = 2ν 0 1 -ν 0 ; λ s + µ s = E 0 h 2(1 -ν 0 ) (69)
Imposing the usual conditions that -1 ≤ ν 0 ≤ 1/2 and E 0 > 0, we have

-1 ≤ λ s /µ s ≤ 2, λ s + µ s > 0, and µ s > 0. ( 70
)
Then, in the absence of body forces in bulk phases and eigenstrains in the inclusion, the weak formulation of the elasticity problem with weak discontinuities characterized by coherent imperfect interfaces is given by [START_REF] Yvonnet | An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites[END_REF] Ω

δε ij C ijkl ε kl dΩ + SI δε ij C s ijkl ε kl dS = ∂Ωt ti δu i dS ( 71
)
ZHU QI-ZHI in which the interface stiffness tensor is expressed in terms of the moduli λ s and µ s C s ijkl = λ s P ij P kl + µ s (P ik P jl + P il P jk ) (72)

where P ij denotes the interface tangential operator P ij = δ ij -n i n j with n i being the unit norm vector defined at a point on S I . δ ij is the second order identity tensor.

Here we present briefly the discretization of the coherent imperfect interface problem under investigation. Applying the standard finite element procedure to the weak formulation (71), we obtain

(K + K s )u = ∂Ω N T tdS, ( 73
)
where u is the column of all system unknowns. The system stiffness matrix contains two parts.

The first part due to the bulk phases is expressed as

K IJ = Ω B T I C i B J dΩ, i = 1, 2. ( 74
)
where B I is a matrix of shape function derivatives related to the XFEM discretization scheme.

The second part K s of the system stiffness matrix due to the interface is provided by

K s IJ = SI B T I C s B J dS, (75)
The computational global convergence is assessed via the L 2 -norm indicator:

ℓ 2 = u -u 0 L 2 (Ω) u 0 L 2 (Ω) = 1 |Ω| Ω u -u 0 2 dΩ 1/2 1 |Ω| Ω u 0 2 dΩ 1/2 , (76)
where |Ω| denotes the volume of Ω and u 0 is calculated by the analytical solution (66). Before presenting numerical results, remark that the displacement field given in Equation (66) is linear in the inclusion but nonlinear in the matrix. Note that because of its feature of outer enrichment, the function ϕ 4 = 1 2 tanh(ρ φx ) sign( φx ) + 1 should be more suitable for describing the current discontinuous problem and therefore will be applied to numerical computations. For completeness, comparisons via global convergences predicted by different enrichment functions will be presented at last stage.

6.1.1. Comparative study on the coefficient ρ We first assess the influence of the coefficient ρ on global convergence. The unit cube under consideration is at nanolength scale, that is, the side length is set as 1nm. The radius of the inclusion is taken to be R 0 = 0.37nm. For simplicity, we use the Young's modulus E = 10MPa and the Poisson's ratio ν = 0.3 for the bulk phases (matrix and inclusion), and λ s = µ s = 5N/m for the interface. The remote hydrostatic strain is prescribed as ε = [1 1 1 0 0 0]. For comparisons, five constant values of ρ, i.e. ρ = 1, 10, 50, 100 and 200 are taken into account. Additionally, since regular meshes are adopted [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF], we take φ c = 1/N c as the characteristic length, N c being the number of layers of subcubes along one dimension. Global convergence rates are reported and compared in Figure 19 for a series of mesh densities with N c ranging from 10 to 50. It is seen from this nanoscale modelling that when ρ is constant and increases, the convergence rate becomes bigger while numerical accuracy gets worse for relatively lower-density meshes. For the cases that ρ = 1 and 10, the convergences decrease progressively towards unacceptable levels. However, the choice ρ = 1/φ c = N c whose value evolves with the refinement of the mesh provides a quite stable convergence.

The two choices ρ = 1 and ρ = N c are then evaluated under various combinations of material constants. Precisely, for the bulk phases, we account for the cases of a rigid inclusion (µ 1 = 10MPa and µ 2 = 200MPa) and of a soft inclusion (µ 1 = 200MPa and µ 2 = 10MPa). For the interface moduli (λ s , µ s ), we consider three combinations (-5, 5), (0, 5) and [START_REF] Belytschko | A review of extended/generalized finite element methods for material modelling[END_REF][START_REF] Belytschko | A review of extended/generalized finite element methods for material modelling[END_REF], among which the first set (λ s = -5 and µ s = 5) gives the relation λ s + µ s = 0, which according to (67), renders the displacement field size-independent. The numerical results shown in Figure 20 confirm that ρ = N c can provide quite satisfactory and stable convergences and also a good compromise between numerical accuracy and computation convergence. Once again, the case of ρ = 1 leads to progressive decrease in convergence, as shown in Figure 21. In conclusion, with an appropriate value of ρ, the smoothing term tanh(ρ φx ) in place of φx can provide significant improvements.

-

Comparisons between enrichment functions

In the previous tests, only the outer enrichment function ϕ 4 is used. We now compare the enrichment functions we derived with those proposed by Moës et al. [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] and Farsad et al. [START_REF] Farsad | An extended finite element/level set method to study surface effects on the mechanical behavior and properties of nanomaterials[END_REF]. Again, we take E = 10MPa and ν = 0.3 for the bulk phases (matrix and inclusion), and λ s = µ s = 5N/m for the interface. The remote hydrostatic strain is set as ε = [1 1 1 0 0 0]. Moreover, ρ = N c is used. The convergence rates are plotted and compared in Figure 22. It is shown that global convergences strongly depend on the choice of enrichment function. More precisely, the functions ϕ 2 = 1 2 tanh(ρ φx ) sign( φx ) -sign(φ r ) , ϕ 3 = 1 2 tanh(ρ φx ) sign( φx ) -1 and ϕ(x, x r ) = φxφr predict, at least on nanoscale, quite low convergences (about R=1) and relatively poor accuracies. However, both functions ϕ 1 = 1 2 tanh(ρ φx ) sign( φx ) + sign(φ r ) and ϕ 4 = 1 2 tanh(ρ φx ) sign( φx ) + 1 are more efficient. Although function sign( φx ) -sign( φr ), the shifted form of sign( φx ), works quite well for strong discontinuities, the shifted form of function φx shows poor capability of achieving satisfactory convergence and accuracy in the weakly discontinuous problem under investigation.

A 2D heterogeneous problem with spring-layer interface

The second numerical application aims at assessing predictive capacity of the discontinuous enrichment functions for modelling heterogeneous problems with strong discontinuities. As shown in Figure 23, the materials occupying the unit square domain consist of a matrix with the elasticity tensor C 2 and a circular inclusion with the elasticity tensor C 1 . The circular interface S I between the bulk phases is assumed to be imperfect and described by the springlayer imperfect interface model. The spring-layer model states that the displacement vector u suffers a jump across S I whereas the traction vector t is continuous across the same interface and proportional to the displacement jump: t = K. u where K denotes the second-order interfacial stiffness tensor. The weak formulation corresponding to the above elastic spring-layer interface problem is given by

Ω δε ij C ijkl ε kl dΩ + SI δu i K s ij u j dS = ∂Ωt ti δu i dS, (77)
where t denotes stress traction on boundary Ω t . For more details on mathematical derivation of the model, formulation of the above weak form as well as its XFEM-framed numerical implementation, one can refer to [START_REF] Zhu | Three-dimensional numerical modelling by XFEM of springlayer imperfect curved interfaces with application to linearly elastic composite materials[END_REF] . Numerically, the geometric domain is first discretized by a regular mesh of triangular linear elements. The mesh density is controlled by the number of segments n along one side. To achieve a good compromise between numerical accuracy and display convenience, we adopt here n = 60. The interface is introduced into the regular mesh through the level set function

f (x, y) = x 2 + y 2 -r = 0, (78)
where r represents the radius of the circular inclusion centered at (0, 0) and takes here r = 0.3. Accordingly, there exist four nodes (0, -0.3), (0, 0.3), (-0.3, 0) and (0.3, 0) locating on the interface, which generate some blending elements of the first type. Therefore, this example allows us to test all the four types of elements (standard elements, fully cut elements and two types of blending elements). The structure is subjected to a uniform displacement u 0 = 0.001 on its side x = 0.5. Both the cases of a rigid inclusion (µ 1 /µ 2 = 15) and a soft inclusion (µ 1 /µ 2 = 1/15) are taken into account. The displacement fields over the deformed configurations are reported in Figure 24. As expected, we observe in the loading direction a more pronounced displacement jump in the case of a rigid inclusion. The smoothness of the crescent-shaped zone of displacement discontinuity implies that the enrichment function and the relevant numerical treatment allow ZHU QI-ZHI to describe correctly the discontinuous displacement field of the spring-layer interface problem in question (otherwise, the field of displacement jump will be saw-toothed). In Figure 25, the distributions of the normal stress are depicted along some chosen lines which pass the center and are inclined at different angles θ with respect to the x-axis. In all the tested cases, the requirement on the continuity of stress traction across the interface is completely satisfied. It is also pointed out that there is no visible difference between numerical results predicted by the enrichment functions ψ i , i = 1, 2, 3, 4, 5. Such observations are not surprising because all these functions allow reproducing the piecewise linear field in question.

SUMMARIES AND CONCLUDING REMARKS

In this paper, we elaborated a consistent two-step procedure for mathematical derivation of enrichment functions in the extended finite element method, in which use has been made of the level set method, the domain decomposition and reproduction technique as well as characteristic functions. The main results and original contributions achieved in this work are summarized as follows:

• We derived mathematically the enrichment functions

ψ i = 1 2 sign( φx )P i ( φx , φ r ), i = 1, 2, 3, 4, 5 (79)
for describing strong discontinuities, and proposed the functions

ϕ i = 1 2 tanh ρ φx P i ( φx , φ r ), i = 1, 2, 3, 4, 5 (80)
for numerical modelling of weak discontinuities, in which the terms P i ( φx , φ r ) dedicated to characterizing relative position between integration point x and node x r are given by

P 1 = 1 + sign( φx )sign(φ r ), (81a)
P 2 = 1 -sign( φx )sign(φ r ), (81b)
P 3 = 1 -sign( φx ), (81c)
P 4 = 1 + sign( φx ), (81d)
P 5 = 1. (81e)
Most of these functions have never been reported in the literature. The basic enrichment functions often used in the XFEM have been found, corrected or improved in this work. Moreover, with the help of the property tanh ρ φx = tanh ρ φx sign( φx ), the following relationship between the two types of functions is established

ϕ i = tanh ρ φx ψ i , i = 1, 2, 3, 4, 5 (82)
The functions ψ i and ϕ i are illustrated and compared in Figure 26 via a 1D element IJ. • For strong discontinuities, two general Sign-type functions and four Heaviside-type functions have also been derived in the same framework.

ZHU QI-ZHI

• The physical meanings of both the basic and complementary nodal variables are now clarified, which can help to impose Dirichlet-type external boundary conditions. • A simple, efficient and unified numerical treatment for the numerical issue of blending elements was proposed, which shows significant advantages in view of the requirement of the partition of unity. Its numerical implementation has also been provided.

• Numerical examples concerning classical discontinuous problems with spring-layer or coherent imperfect interfaces have been presented. Through global convergences and field distributions, we validated and compared the predictive ability of the derived enrichment functions in describing and modelling heterogeneous problems with discontinuities.

From the above-mentioned aspects, we can address the following comments on XFEM enrichment functions:

• The two classes of enrichment functions, which share the terms P i , can be constructed in the same way. Precisely, functions ψ i are obtained by multiplying P i with the characteristic function sign( φx ), while ϕ i are further formulated by multiplying ψ i with a smoothing factor, such as, φx or tanh( φx ). • This study confirms that the Sign function and the Heaviside function themselves as well as their respective shifted forms are all suitable for describing strong discontinuities.

There is in fact no essential difference between these two types of enrichment functions. • For the heterogeneous problem with strong discontinuities, when global convergence and field distribution are concerned, we found no visible differences between the derived enrichment functions ψ i . Nevertheless, the inner enrichment function ψ 3 ( φx ) shows particular interest for dealing with multiple cuttings and/or complex geometries. • For heterogeneous problems with weak discontinuities, at least in the tested benchmark with coherent interfaces at nanolength scale, computational convergences and accuracies strongly depend on the choice of enrichment function. The proposed improvement seems to be necessary to achieve satisfactory convergence and accuracy. • From numerical point of view, the two Heaviside-type step functions used in the Hansbo and Hansbo's work [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] for constructing two one-sided independent interpolations can be merged into one single function ψ 1 (x, x r ).

Although the present work mainly focused on mathematical derivation of basic enrichment functions in the XFEM, some useful extensions and direct applications are realized. Based on these results, further investigations along this line may include the following aspects:

• For multiple cuttings, the formulation (44) will be implemented in the XFEM context.

Validations will be performed via benchmarks and applications to composite materials with a large number of inclusions will be carried out. • Triple (or even more) junctions may take place when material interfaces meet at one point, like that often encountered in polycrystals ( [START_REF] Simone | A generalized finite element method for polycrystals with discontinuous grain boundaries[END_REF] [9] [10] [START_REF] Mousavi | Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons[END_REF]). When limited to the case where no intersections between interfaces occur, a triple junction can be regarded as the asymptotic case of a three-cutting configuration, like that depicted in Figure 9.

Once interfaces are well defined by the level set method, it is then possible to adapt the formulation (44) to model triple junctions. Undoubtedly, local integration for an element containing a triple junction would be more complicated. On this topic, one can refer to the work by Natarajan et al. [START_REF] Natarajan | Numerical integration over arbitrary polygonal domains based on Schwarz-Christoffel conformal mapping[END_REF] and that by Mousavi and Sukumar [START_REF] Mousavi | Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons[END_REF]. • In failure processes due to arbitrary cracking in heterogeneous materials, intersections between static and evolving interfaces often take place. In future work, it will be our particular interest to model competing process between progressive debondings at imperfect interfaces and propagation of distributed cracks in the matrix. Locally anisotropic propagation of microcracks in the matrix can be modelled by using homogenization-based multi-scale methods (e.g. [START_REF] Zhu | Homogenization-based analysis of anisotropic damage in brittle materials with unilateral effect and interactions between microcracks[END_REF]) whereas numerical issue of macrocracks formed by coalescence of microcracks will be addressed within a combined XFEM/LSM context. For the latter, one can benefit from the cutting algorithm used in [START_REF] Richardson | An XFEM method for modeling geometrically elaborate crack propagation in brittle materials[END_REF] for describing multiple interacting and branched cracks.

ACKNOWLEDGEMENT

The author would like to thank ANR (the French National Agency for Research) for the partial financial support received through research project MICROFISS (ANR-07-BLAN-0300). Instructive suggestions from professor Qi-Chang He at Université Paris-Est and helpful discussions with Dr. Shui-Tao Gu during his postdoctoral stay in France (currently at Southwest Jiaotong University of China) are gratefully acknowledged.

APPENDIX

By taking α = 1 in the characteristic functions ( 12) and ( 13) and following the derivation procedure elaborated in Section 3.4, we can derive the following Heaviside-type enrichment functions for describing strong discontinuities ψ1 (x, x r ) = H 1 ( φx) + H 1 (φr) -1, (83) ψ2 (x, x r ) = H 1 ( φx) -H 1 (φr), (84)

ψ3 (x) = H 1 ( φx) -1, (85) ψ4
(x) = H 1 ( φx), (86)
which correspond to domain decompositions (D-I), (D-II), (D-III) and (D-IV), respectively. Obviously, both the Heaviside function itself ψ4 (x, x r ) and its shifted form ψ2 (x, x r ) are derived in the present framework. It is also proved that the above functions ψi lead to the same displacement jump vector (47) To facilitate the imposition of Dirichlet boundary conditions, we summarize in Table II the physical meanings of both the basic and complementary nodal unknowns.

Table II. Imposing Dirichlet boundary conditions for an enriched node x i

Case ψ 1 (x, x r ) ψ 2 (x, x r ) ψ 3 (x) ψ 4 (x) all 4 eqns.

x i ∈ I * φ i > 0 a i = u - i a i = u + i a i = u + i a i = u - i b i = u + i -u - i but x i / ∈ S I φ i < 0 a i = u + i a i = u - i a i = u + i a i = u - i b i = u + i -u - i x i ∈ S I φ i = 0 a i = u - i a i = u + i a i = u + i a i = u - i b i = u + i -u - i

Figure 2 .

 2 Figure 2. (a) A triangular linear element Ωe which is geometrically divided by the interface segment S αβ into two parts, (b) domain decomposition of the linear but discontinuous displacement field over Ωe with natural nodal values u i , u j and u k , (c) basic and complementary nodal variables generated by extending linearly the subfields u + (x) and u -(x) towards the whole triangle.

Figure 4 .

 4 Figure 4. Three-part domain decomposition of the displacement field with a = u - i , u + j , u - k

Figure 5 .

 5 Figure 5. Three-part domain decomposition of the displacement field with a = u + i , u - j , u + k

Figure 6 .

 6 Figure 6. Two-part domain decomposition of the displacement field with a = u + i , u + j , u + k

Figure 8 .

 8 Figure 8. (D-I ) domain decompositions of the displacement field in the case that node x i locates on the interface

Figure 9 .

 9 Figure 9. A triangular linear element cut by three interfaces across which discontinuities occur: a) geometric representation, b) domain decomposition

Figure 10 .

 10 Figure 10. Three dimensional cutting configurations for a linear tetrahedral element

Figure 11 .

 11 Figure 11. Illustration of the enrichment function ψ 1 (x, x r ) on the cut element Ω ijk

Figure 12 .

 12 Figure 12. Illustration of the enrichment function ψ 2 (x, x r ) on the cut element Ω ijk

Figure 13 .

 13 Figure 13. Illustration of the node-independent enrichment function ψ 3 (x), ψ 4 (x) and ψ 5 (x) on the element Ω ijk cut by the interface segment Γ αβ

Figure 17 .

 17 Figure 17. Illustration of a cut triangular element Ω ijk used to analyze the enrichment function proposed in [13].

Figure 19 .

 19 Figure 19. Global convergence rates R for various choices of ρ: calculations performed on the nanoscale with the enrichment function ϕ 4 .

)

  Cube with the side length of 1nm rigid, λ s = -5N/m, µ s = 5N/m (R = 1.59) rigid, λ s = 0N/m, µ s = 5N/m (R = 1.60) rigid, λ s = 5N/m, µ s = 5N/m (R = 1.60) soft, λ s = -5N/m, µ s = 5N/m (R = 1.36) soft, λ s = 0N/m, µ s = 5N/m (R = 1.45) soft, λ s = 5N/m, µ s = 5N/m (R = 1.49)

Figure 20 .

 20 Figure 20. Convergence rates R for ρ = Nc and various combinations of material constants: µ 1 = 10MPa and µ 2 = 200MPa for a rigid inclusion while µ 1 = 200MPa and µ 2 = 10MPa for a soft inclusion.

Figure 21 .

 21 Figure 21. Convergence rates R for ρ = 1 and various combinations of material constants: µ 1 = 10MPa and µ 2 = 200MPa for a rigid inclusion while µ 1 = 200MPa and µ 2 = 10MPa for a soft inclusion.

6 log 10 ( 1 /Figure 22 .

 610122 Figure 22. Comparisons between the enrichment functions through global convergence rate.

Figure 23 .

 23 Figure 23. A two-dimensional spring-layer imperfect interface problem: geometric descriptions, boundary conditions and loading

Figure 24 .Figure 25 .

 2425 Figure 24. Displacement distributions on deformed structures in which elements in red denote the cut or blending elements and elements in blue are the standard elements

ÁFigure 26 .

 26 Figure 26. Illustrations of the derived enrichment functions, which are performed on the 1D element IJ. The functions ψ 1 , ψ 2 , ϕ 1 and ϕ 2 are node-dependent while the others are node-independent.

Copyright c 2011 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2011) Prepared using nmeauth.cls DOI: 10.1002/nme

ZHU QI-ZHI or more explicitly ϕ 1 ( φx , φ r ) = 1 2 φx sign( φx ) + sign(φ r ) , (58a) ϕ 2 ( φx , φ r ) = 1 2 φx sign( φx ) -sign(φ r ) , (58b)

Remarkably, enrichment functions for weak discontinuities are linked to those for strong discontinuities via Equations (57). Above, function ϕ 5 is further written in the form

which can also be obtained from linear combination of ϕ 1 and ϕ 2 or that of ϕ 3 and ϕ 4 .

Remarks

• Referring to the particular case depicted in Figure 16, the edge S 12 locating on the interface is shared by the elements Ω 124 and Ω 128 . It is shown that the enrichment function proposed by Moës et al. [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] ϕ

takes zero value for integration points inside Ω 124 and Ω 128 . Consequently, the enrichment term will be removed for these two elements and function (60) may fail in modelling weak discontinuity across the interface segment S 12 . Remark that the enrichment functions ϕ i , i = 1, ..., 5 we elaborated allow remedying this shortcoming. • Recently, in order to model surface effects on the mechanical behaviors and properties of nanomaterials, Farsad and his co-workers [START_REF] Farsad | An extended finite element/level set method to study surface effects on the mechanical behavior and properties of nanomaterials[END_REF] proposed another enrichment function, which reads as the shifted form of the characteristic function φx