Convergence behaviour of the enriched scaled boundary finite element method

In this work, a very efficient numerical solution of three‐dimensional boundary value problems of linear elasticity including stress singularities is discussed, focussing on its convergence behaviour. For the employed scaled boundary finite element method, a discretization is only needed at the boundary, while the solution is considered analytically in a scaling coordinate. This presents a major advantage for two‐dimensional problems, when the scaling center is placed at a stress singularity. Unfortunately, three‐dimensional problems usually do not only include point singularities but also line singularities, which results in singular gradients in the boundary coordinates and thereby diminishes the method's original advantages. To alleviate this drawback, this work discusses an enrichment of the separation of variables representation with analytical asymptotic near fields of the line singularities. In contrast to previous works, besides the near‐field functions with λ=0.5, also those with λ=1.5 were determined and used for enrichment. This leads to a high accuracy and it is shown that this approach is required to recover the convergence properties of smooth boundary value problems without singularities when using quadratic Lagrange shape functions. In order to recover the convergence rates for higher order shape functions, near‐field functions with higher singularity exponent have to be included for enrichment.


Summary
In this work, a very efficient numerical solution of three-dimensional boundary value problems of linear elasticity including stress singularities is discussed, focussing on its convergence behaviour. For the employed scaled boundary finite element method, a discretization is only needed at the boundary, while the solution is considered analytically in a scaling coordinate. This presents a major advantage for two-dimensional problems, when the scaling center is placed at a stress singularity. Unfortunately, three-dimensional problems usually do not only include point singularities but also line singularities, which results in singular gradients in the boundary coordinates and thereby diminishes the method's original advantages. To alleviate this drawback, this work discusses an enrichment of the separation of variables representation with analytical asymptotic near fields of the line singularities. In contrast to previous works, besides the near-field functions with = 0.5, also those with = 1.5 were determined and used for enrichment. This leads to a high accuracy and it is shown that this approach is required to recover the convergence properties of smooth boundary value problems without singularities when using quadratic Lagrange shape functions. In order to recover the convergence rates for higher order shape functions, near-field functions with higher singularity exponent have to be included for enrichment.

KEYWORDS
convergence, enriched base functions, enriched scaled boundary finite element method (enrSBFEM), scaled boundary finite element method (SBFEM), three-dimensional elasticity are independent of the polynomial degree of the employed shape functions. 3,[7][8][9] In order to achieve better accuracy and rates of convergence, an appropriate adaptation or a skillful choice of the discretization and element formulation within the FEM is necessary. Different strategies have been developed, so far. In other works, 4,10-17 the employment of a domain rescaling in the vicinity of the singularity, such that the singular near field is transformed into a smooth one, is proposed. In order to achieve a better representation of the solution, eg, in other works, 8,[18][19][20][21][22][23][24][25][26][27][28][29] the finite element space is chosen in such a way that not only the commonly used polynomial functions are in this space but also additional functions which should represent the displacement field in the vicinity of the singularity (enriched FEM, XFEM). Another option is to perform the discretization only in the circumferential direction, while the solution of the BVP is considered analytically in a radial direction toward the singularity (FEM eigenanalysis techniques [30][31][32][33][34][35][36][37][38][39][40][41] ; scaled boundary FEM (SBFEM) [42][43][44][45][46][47][48][49][50][51][52][53] ). The SBFEM was successfully applied to fracture mechanics problems, like the analysis of asymptotic near fields, the evaluation of generalised stress intensity factors [54][55][56] and automatic crack growth algorithms. [57][58][59][60] The application of the method to other physical problems is, eg, presented in other works. [61][62][63][64][65][66][67][68] The SBFEM is an appropriate method for the treatment of two-dimensional (2D) crack problems, when the singularity is located entirely within the considered domain. However, in three-dimensional (3D) problems in fracture mechanics, stress singularities can still occur at the point, where the crack front intersects the discretized boundary. Because the solution at this point is determined numerically, similar problems as in the standard FEM arise. In order to solve this problem, other authors 69-77 employed a mesh refinement toward the singular point in the circumferential coordinates. Apel et al 78 employed the graded mesh approach such that the convergence orders of a smooth problem were obtained again for displacement exponents as well as deformation modes. However, while such graded meshes yield the desired convergence rates, they generally are not very efficient regarding the necessary degrees of freedom. A more efficient approach, concerning the degrees of freedom, has been proposed in the works of Hell and Becker 79,80 where stress singularity exponents were determined using an enriched formulation of the SBFEM (enrSBFEM) and then compared, within the frame of a convergence study, to the results obtained from a standard formulation.
This work ties in with these forgoing ones and the previous convergence studies are extended to include the 3D deformation modes and the solution of full BVPs. This formulation is described in detail in Section 2.3 after prerequisites of the standard FEM using Lagrange shape functions in all spatial dimensions are briefly revisited in Section 2.1 (including convergence issues in the presence of stress singularities). The basics of the standard 3D SBFEM are revisited in Section 2.2. Numerical results concerning the convergence properties of the enriched SBFEM formulation are thoroughly studied and compared to the standard SBFEM formulation in Section 3. In Section 3.1, the implementation and numerical models (including geometry, enriched subdomains, etc.) are described in detail. The results in Section 3.2 refer to an example of a single crack in a homogeneous isotropic continuum. As there exists an analytical solution for this case, also the error in the associated deformation modes is studied. Finally, some conclusions on the enrSBFEM and its performance are drawn in Section 4.

Prerequisites
This section includes a short introduction into the established standard FE-approach which makes use of a discretization in all spatial dimensions and Lagrange shape functions. The corresponding estimates of discretization errors are also shortly revisited.

Problem.
Many physical processes can be modelled using partial differential equations (PDEs). In linear elasticity theory for static problems, there also exists a system of elliptic PDEs. Including boundary conditions, it can be written as with Ω = Γ N ∪ Γ D , Γ D ≠ ∅, where Ω ⊂ R 3 is the region covered by the solid, f are the body force intensities, t the surface tractions, n is the outward surface normal, and u ∶ Ω → R 3 the displacements. Finally, ∇u ∶ Ω → R 3×3 is the gradient of the displacements and (u) = 1 2 (∇u + ∇ T u) =∶ ∇ s u the linearised symmetric strain tensor, which is valid for small deformations. The stress tensor can be expressed by the constitutive law (u) = C ∶ (u), where C ∶ R 3×3 s → R 3×3 s represents the elasticity tensor. ‡ Finite element methods provide a numerical approach for the solution of system (1)- (3). The established standard FE approach, with a spatial discretization in all spatial dimensions and the use of Lagrange shape functions, as well as its variations are based on a weak formulation of system (1)-(3). In the following, H m (Ω) denotes the Sobolev-space where is a multiindex = ( 1 , … , n ), | | = ∑ n i=1 i is the length of the multiindex, and D v is the 'th weak derivative of v. For further details, please refer to the literature. 81,82 Suppose V ∶= V(Ω) ∶= {v ∈ H 1 (Ω), v| Γ D = 0} and V ∶= V 3 , then the weak formulation of (1)-(3) reads:

Standard Finite Element Discretization.
A common finite element space V h ⊂ V as an N-dimensional subspace of V, N < ∞, is used for the discretization in each spatial dimension, so that the Cartesian product Then, the Galerkin approximation reads: Find Let With the notations a i ∶= a( i , ), A ∶= (a i ) ∈ R 3N×3N and b ∶= l( ), b = (b ) ∈ R 3N , we get the linear equation system Often a piecewise polynomial ansatz is employed as follows: where  p (T i ) is the space of polynomials of degree p in the subdomain T i . Most commonly, Lagrange polynomials are used as basis functions, ie, a basis function i has to be in V (p) h and fulfil i (x j ) = ij with supports at x j , j = 1, … , N, where ij denotes the Kronecker delta.

Error estimates.
The following norms and seminorm are defined (cf, eg 82,83 ): In this section, we assume Ω ⊂ R 2 and use L 2 as a short notation of L 2 (Ω). Furthermore, a simplified case of the weak form (5) is considered using only one-dimensional (1D) displacements u ∶ Ω → R, vanishing surface tractions t ≡ 0 and the identity tensor as an elasticity tensor C = I. Accordingly, the weak formulation here reads For the error of a FEM with finite-dimensional trial spaces V h ⊂ V, the following estimates hold: where p is the polynomial degree of the basis functions of the trial space (see, eg 82 ). Accordingly, at a uniform refinement of the discretization, the displacement error in L 2 -norm converges to zero with the order p + 1. This holds true as long as solutions u ∈ H 1 (Ω) are in H p+1 (Ω) and therefore ||∇ p+1 u|| L 2 (Ω) is bounded. This is not generally the case, eg, when Ω exhibits notches or cracks and the displacement solutions contain singularities. Then, restriction (12) is weak as the right-hand side is unbounded and thus does not yield an actual upper bound for the error anymore. With Equation (12), we have error estimates for the FEM with piecewise polynomial shape functions at hand. In order to achieve a reasonable estimate for linear shape functions, u has to be at least in H 2 (Ω), because the right-hand side otherwise becomes unbounded. The leading term of the asymptotic near-field solution, which is not a rigid body motion, contains a factor √ r . For Ω ⊂ R 2 , √ r certainly is in H 1 (Ω), but not in H 2 (Ω). Consequently, the standard error estimate is useless. In this case, a different estimate showing that the rate of convergence is bounded to a measure related to the strength of the singularity has to be considered (cf other works [7][8][9] ). There also is an estimate for an enriched FEM which makes use of a displacement representation including piecewise polynomial shape functions extended by enrichment functions. It states that, for an appropriate choice of enrichment functions, the optimal rates of convergence, being those of a problem without singularities, can be recovered. In the following, these error estimates are briefly revisited.
Each component of the displacement solution in the vicinity of a crack tip can be represented by an appropriate asymptotic expansion (see, eg, other works [84][85][86] ) for i = x, y, z, whereû i0 is a rigid body motion. For increasing n, the summands in the solution representation are getting smoother. ∶= 1 = 0.5 is designated as the leading singularity. Equation (13) can also be expressed as a finite sum and a smooth remainder U 8 The differentiability conditions of U depend on the choice of S (the larger S, the smoother U Without enrichment. For a FEM with piecewise polynomial trial spaces V h (see Equation (8)), in the case of occurring singularities with leading singularity exponent 0 < 1 = < 1, the following error estimates, as proven in the work of Blum, 8 hold: It is worth noting that the rates of convergence are independent of the polynomial degree of the employed shape functions.
With enrichment. For a FEM with piecewise polynomial trial spaces with polynomial degree p, which is extended by enrichment functions up to a singularity exponent S , the following error estimates hold 8 : The estimates (17) result from a generalisation of the estimates (16). According to them, singular functions up to an exponent of S+1 > p have to be included into the trial space to retrieve optimal orders of convergence.

Scaled boundary finite element method
For the SBFEM a star-convex domain is presupposed. In the star center the so-called scaling center O is placed and, there, a new coordinate system is introduced: the scaled boundary coordinate system. It consists of a radial coordinate running from the scaling center toward the domain's boundary, and (in the 3D case) two surface coordinates 1 , 2 . In contrast to the standard FE approach, there is no need to discretize the total domain Ω. For the SBFEM, we only have to introduce a discretization in the surface (or boundary) coordinates of the considered domain. Consequently, the problem dimension is reduced by one and the number of nodes is thereby reduced drastically compared to a standard FEM. Application of the FEM in 1 -, 2 -direction leads to an ordinary PDE in variable , which can be solved analytically.
We define the following differential operator: Using Voigt notation, here denoted as V , V , C V , the representation V (u) = u can be obtained for the strains and In the following, the Voigt notation will be employed and the index · V omitted for brevity.
In order to apply the SBFEM, a transformation T from the global x-, y-, z-coordinate system to the -, 1 -, 2 -coordinate system is required. The scaling variable is zero at the scaling center O = (x 0 , y 0 , z 0 ) and one at the boundary. The boundary coordinates 1 and 2 are tangential to the surface of domain Ω ( Figure 1). The corresponding transformation reads where [x , y , z ] ∶ R 2 → R 3 is a parametric surface representation of the domain boundary. The transformation T is characterized by Therefore, the differential operator of Equation (18) can be decomposed in the following manner: Here, the entries of the matrices L x , L y , and L z consist of zeros and ones and b , b 1 , and b 2 are obtained from Equation (20) (ie, from the inverse of J). The definitions

FIGURE 1
Transformation from cartesian to scaled boundary coordinate system, with scaling variable and boundary variables 1 , 2 and transformation from the boundary variables to the reference element then lead to the following representation: Only the boundary of the region Ω is discretized so that the shape functions i only depend on the two boundary coordinates 1 , 2 (3D case; in the 2D case, there naturally is only one boundary coordinate ). The resulting equations will be solved analytically with respect to the scaling variable . To achieve this, a separation of variables representation is employed. Together with the basis functions with coefficient functionsū i ∶ R → R. Using an analogous approach for virtual displacements δu h ∈ V h , which serve as test functions, leads to and accordingly ∫ The separation of variables representations of u h and δu h as well as Equation (23) are applied in Equation (26). After that, under transformation (19), by use of the matrices where J = J( 1 , 2 ) and after partial integration, we obtain The assumption of vanishing body forces f = 0 and again the application of transformation (19) for the linear form yields Using Equations (28) and (29), the following relations can be derived from the weak formulation (26): (30) represents a Cauchy-Euler differential equation system of second order. Its solution yields a number of deformation modes fulfiling the field equations as well as boundary conditions along the boundary faces which originate at the scaling center. The linear equation system (31) serves to additionally enforce the boundary conditions on the discretized boundary.
Substituting = e t and introducing the functionsv(t) ∶=ū ,t (t) in the differential equation system (30), we get a system of ordinary differential equations (ODEs) of first order with constant coefficients: To solve this, the exponential where the equivalence holds as e t ≠ 0 for all t ∈ R. In conclusion, we have to solve an eigenvalue problem to determine the solutionw(t). Suppose that n is the number of eigenvalues of . Furthermore, it is assumed that the geometric multiplicity is equal to the algebraic one for each eigenvalue and that the eigenvectors are linearly independent. This is a valid assumption in 3D problems. Then, the solution of the ODE-system (32) is the linear combination of the fundamental solutions. After a back-substitution t = ln , the solutionū( ) of Equation (30) is with coefficients c i ∈ C to be determined from Equation (31). The eigenvectors ( u ) i describe the nodal displacements at the boundary = 1. We also denote them deformation modes as they contain the essential deformation information.
The respective eigenvalues i can be interpreted as their associated decay (Re( i ) < 0) or growth rates (Re( i ) > 0) in the scaling direction . 49 This makes representation (34) similar to asymptotic expansions usually employed in the analysis of stress singularities. The spectrum of eigenvalues is symmetric, in the 3D case with respect to −0.5. There are as many eigenvalues in every half of the spectrum as there are DOF due to discretization. However, as long as only an outer boundary is considered, ie, the scaling center is located within or on the considered body, the coefficients c i corresponding to ≤ 0 need to be zero to keep the displacement bounded at the scaling center (regularity condition). According to Equation (23), a derivative of the displacements with respect to , respectively their division by is needed to calculate the stresses. Consequently, in the representation of the stresses, the exponent of is reduced by one, ie, stress singularities occur for exponents Re( i − 1) < 0. In the context of the current work, we will refer to i as the displacement exponents and i − 1 simply as the stress exponents. § In summary, solution (34) is obtained by application of the SBFEM, in which the eigenvalues i and the eigenvectors ( u ) i follow from the differential equation system (30), whereas the coefficients c i follow from the linear equation system (31) for enforcing the boundary conditions.

Enrichment of the SBFEM
Especially in the 2D case, the SBFEM has proven to be a powerful method for the treatment of cracks and stress singularities (see, eg, the work of Song et al 48 ). For 3D problems, the method is still advantageous because of the accompanying § Re(| i − 1|) are commonly called stress singularity exponents. Note that the spectrum of lambda is symmetric to −0.5. Deformation modes with lambda smaller than −0.5 only become active, when an inner scaled boundary exists (eg, at = 0.1). This case is not considered within the scope of this work. information about the eigenvectors and eigenmodes of the BVP and its efficiency caused by the reduced number of degrees of freedom since not the whole domain but only the boundary is discretized. However, in 3D problems of fracture mechanics, stress singularities do not only occur at points but also along lines (eg, along crack fronts and notch fronts). Points of stress singularities with a true 3D decay behaviour then usually occur at discontinuities of such singularity lines. As the SBFEM approach can only account for point singularities (by placing the scaling center there), the singularity lines remain to be represented by the shape functions in the boundary coordinates. Consequently, in such cases, the 3D SBFEM must suffer from similar drawbacks as the standard FEM approach.
The objective of the current work is to retrieve the improved accuracy and rate of convergence with as few additional degrees of freedom as possible. The idea of enrichment is to choose the trial space in such a way that not only the commonly used polynomial shape functions are in this space but also additional trial functions which should represent the displacement in the vicinity of singularity lines. In order to get appropriate additional trial functions for cases involving cracks in homogeneous isotropic materials, we use the classical crack tip near field with displacement exponent = 0.5. With respect to the crack tip coordinates r and from Figure 2, the classical crack tip near field reads 86 in the x ′ -, y ′ -, z ′ -coordinate system ( Figure 2B), wherein denotes the shear modulus. Each of its components can be represented by a linear combination of the basis functions 26 The functions f 1 (r, ), … , f 4 (r, ) have to be represented in terms of 1 and 2 . Accordingly, an additional transformation from r-, -to 1 -, 2 -coordinates is needed. In addition to Equation (35), which represents the deformation modes with respect to the eigenvalue = 0.5, further deformation modes can be determined by means of the method of complex potentials. 86 The deformation modes with respect to the eigenvalue = 1.5 arê K IV cos 2 ( + 2 − (2 + 1) cos ) + K V sin 2 ( + 4 + (2 + 5) cos ) K IV sin 2 (− + 2 + (−2 + 1) cos ) + K V cos 2 ( − 4 + (−2 + 5) cos ) K VI sin 2 (2 cos + 1) For every point R, where the crack front meets the discretized boundary, nodes in its neighbourhood are chosen to be enriched. Let I = {1, … , N} be the index set of the nodes of the finite element discretization. Then, I R enr ⊂ I denotes the index set of the nodes, which should be enriched, close to the point R. We assume that I s enr ∩ I t enr = ∅ for s ≠ t, so that a node cannot belong to more than one enriched domain at the same time. Suppose { 1 , … , N } is the Lagrange basis of a piecewise polynomial function space V (p) h ( Ω) (8) defined in the boundary coordinates 1 , 2 . A very general approach for the enrichment in the standard FEM (cf, eg, the work of Mohammadi 29 ) is employing a set of enrichment coefficients for every enriched node. The transfer of this approach to the SBFEM leads to the following enriched separation of variables representation for a component u h of u h : Here, n R is the number of singular points to be enriched, is the scaling coordinate, and = ( 1 , 2 ) are the boundary coordinates.
The enrichment approach (38) introduces a large number of additional free functions in . This number can be reduced by a factor of n R ∕ ∑ n R R=1 |I R enr | by supposing that the enrichment coefficient functionsā R i1 ( ), … ,ā R i4 ( ) are equal for every enriched node i ∈ I R enr : Let (1, 1 ), … , (1, N ) be the nodal coordinates of the discretized boundary. In order to keep the interpolation property This approach was developed and first employed 79,80 and later also applied 87,88 The effect of the term ∑ i∈I R enr i is called blending (see, eg, the work of Mohammadi 29 ). Because of the choice of i as Lagrange basis, ∑ i∈I R enr i ( ) = 1 will hold if (1, ) is located within an element which has all of its nodes enriched, and ∑ i∈I R enr i ( ) = 0 will hold if (1, ) is located within an element that has none of its nodes enriched. In elements containing enriched nodes as well as not enriched nodes, the function decreases from one (enriched node) to zero (not enriched node) with the smoothness of the chosen polynomial shape functions of degree p (Figure 3).
A set of basis functions of the enriched function space defined on the discretized boundary (for all spatial directions) V enr,

FIGURE 3 Exemplary illustrations for
blending of an enrichment with dimension 3N + 12n R . Therefore, the number of DOF in the eigenvalue problem is increased by four per spatial dimension and crack, which results in a total increase of only 12n R . Thereby, an ansatz for the complete discrete solution u h reads The enrichment with a decomposition of the crack tip field was first introduced by Fleming et al 26 for an element-free Galerkin formulation. In this work, we will refer to this enrichment method as decomposition method. Alternatively, the enrichment can be conducted by directly using the complete crack tip near field (35) as originally proposed by Benzley 18 and Gifford and Hilton 19 for a 2D FEM. We use the definitions and again the transformation from r-, -to 1 -, 2 -coordinates. Together with we get the following basis: V enr, To this enrichment method, we will refer as the direct method in the following. Here, the dimension is 3N + 3n R , which means that, in comparison with the decomposition method, even less additional DOF are obtained in the eigenvalue problem (only 3n R instead of 12n R ).
Additionally, this reduction makes it feasible to also include further crack modes in the enrichment. In the case of an enrichment up to = 1.5, the enrichment functions are additionally appended to Equation (44) (see Equation (37)), so that the number of DOF in the eigenvalue problem then amounts to 3N + 6n R . Let generally n enr denote the number of enrichment functions. ¶ Then, the enriched separation of variables representation of the direct method reads Due to the direct incorporation of the crack modes as enrichment functions, the weighting functionsā ( ) can be interpreted as the corresponding weights of the crack modes and, in case of j < 1, even as generalised stress intensity factors.
we again obtain the common form of the SBFEM separation of variables approach (cf Equation (24)), which can be employed in accordance with the SBFEM procedure, described in Section 2.2.

Implementation and numerical modelling
In the following, the SBFEM with and without enrichment will be compared. A homogeneous isotropic 3D continuum # is considered and a cubic subdomain Ω with corners {[x, y, z]|x, y, z ∈ {−1, 1}} is modelled using the SBFEM. The cube exhibits one crack, whose crack front meets the discretized boundary at x 1 = [1, 0, 0] and x 2 = [−1, 0, 0] ( Figure 4). The crack is modelled by doubling all nodes at the crack faces and assigning them to a crack side (cf Figure 5). The scaling center is placed at x 0 = [0, 0, 0]. For every point x R , R = 1, 2, where the crack front meets the discretized boundary, all nodes were enriched (ie, all nodes on a corresponding side of the cube except for the ones on the cube edges). In order to determine the rate of convergence, a uniform discretization (equally sized quadrilateral elements) is employed in each case. Discretizations n e = 2, 4, … , 16, 20, … , 32 with n e being the number of elements per edge of the considered cubic domain were applied. With an edge length of 2 mm, this corresponds to mesh sizes h = 2 mm n e . ¶ Enrichment up to = 0.5: n enr = 12n R using the decomposition method and n enr = 3n R using the direct method. Enrichment up to = 1.5: n enr = 24n R using the decomposition method and n enr = 6n R using the direct method. # The material data of steel E = 210 GPa, = 0.3 are employed, but due to the assumption of homogeneous material properties, the results are independent of Young's modulus E and in fact only moderately dependent on the chosen value of Poisson's ratio . We considered various modes and respective displacements at the boundary. The convergence of the following norms will be investigated: to which we will simply refer as the L 2 -norm and where || · || is the Euclidean vector norm, and which we will simply denote the energy norm. As only the boundary is discretized, the norms are also determined by only integrating over the boundary coordinates 1 , 2 .
Error estimates for the standard FEM on subsets of R 2 (and the case of 1D displacements) were revisited in Section 2.1. The resulting expected rates of convergence are shown in Table 1 for linear as well as quadratic Lagrange shape functions without and with enrichment, ie, with extension of the trial space using deformation modes up to a displacement exponent of S = 0.5, 1.5. According to these error estimates, the order of convergence without enrichment is independent of the polynomial degree of the employed shape functions limited to (h) in the L 2 -norm and (h 1∕2 ) in the energy norm. When using linear shape functions with enrichment, the optimal order of convergence can already be recovered by an enrichment with deformation modes up to the strongest singularity with displacement exponent S = 0.5. When using quadratic shape functions, an enrichment with singular functions up to a displacement exponent of S = 1.5 is necessary.

No enr.
With enr. up to S = 0.5 With enr. up to S = 1.5    Table 1 for comparison with the rates of convergence that we obtained in our convergence studies, measured in norms (52) and (53). As will be seen from our observations, this assumption seems to be reasonable.

Numerical Integration.
To determine the matrices (27) on the one hand and the norms in postprocessing on the other, a Gaussian quadrature rule has been employed for numerical integration. While this is the standard procedure for nonenriched displacement representations (working very well in this case), a problem arises in the case of using an enrichment.
The matrix of shape functions reads The integrands of the matrices E 1 and E 2 in Equation (27)  . Consequently, near the crack front (r → 0), a singularity occurs so that Gaussian quadrature is a rather suboptimal integration procedure. The determination of the energy norm (53) in postprocessing involves a very similar problem. An integration scheme, which is well suited in this case and also easy to implement, is a quasi polar integration method as proposed by Laborde et al. 89 This is because the integration in polar coordinates includes an area differential rdrd , whose factor r cancels out the factor √ r −1 √ r −1 from the matrices B in E 2 . Quasi polar integration is applied in the elements adjacent to the crack front. That means, if an element node lies on the crack front, the integrals on the quadrilateral element will be transformed to two integral domains of triangular shape, so that the tip of the triangles coincides with the crack front (cf Figure 6). Then, in each of the triangular domains, the same number of integration points (10×10) as in a quadrilateral domain is considered (thus, the domain of an isoparametric quadrilateral is mapped to a triangular domain). In all the other elements which contain enriched nodes (fully enriched as well as blending elements), the occurring gradients are not as severely high as in the very vicinity of the crack front, so that a simple 10×10 Gauss integration seemed sufficient. Finally, for all the other elements, only a 3×3 Gauss integration was employed.

Convergence of the deformation modes in the BVP
In contrast to the displacement exponents , the convergence behaviour of the deformation modes ( u ) i cannot be studied directly. As eigenvectors, they generally are freely scalable and, consequently, need to be appropriately normalised before further investigation. This, eg, can be done by solving a full BVP, ie, also requesting the fulfilment of the boundary conditions at the discretized boundary (Equation (31)), thereby determining the weighting factors c i (in (34)), and finally evaluating the product c i ( u ) i instead of ( u ) i alone. Of course, when only one deformation mode shall be studied, the boundary conditions at the discretized boundary must be defined appropriately and this can easily be implemented, when the exact analytical representation of the deformation mode is known. Then, the boundary conditions are chosen such that the resulting displacements represent the analytical deformation mode of interest, which in this case is one of the crack modes.

TABLE 2
Obtained rates of convergence of the error ||u − u h || measured in || · || L 2 ( Ω) -and || · || a( Ω) -norm for C4-elements without (SBFEM) and with (enrSBFEM) enrichment (with 2D = 0.5) Using the analytical solution of the crack modes I-VI as a reference, in the following, the obtained orders of convergence of the error norms (52) and (53) (L 2 -norm and energy norm) are studied. The results with and without enrichment are compared for the three crack opening modes associated to the eigenvalue I-III = 0.5 (see Figure 4, first row), as well as to IV-VI = 1.5 (see Figure 4, second row). The respective analytical solutions are given in Equations (35) and (37). In order to produce a plane strain condition, for modes I, II, IV, and V, homogeneous Dirichlet boundary conditions are specified at x = −1 and x = 1 in x-direction (u x = 0). For modes III and VI, which only include a displacement in x-direction, homogeneous Dirichlet boundary conditions are specified in y-and z-directions (u y = u z = 0). Additionally, to completely inhibit any displacement in x-direction for modes I, II, IV, and V and a displacement in y-and z-directions for modes III and VI, the following boundary conditions for the enrichment coefficients a j (see Equation (48)) have to be specified a = 0 , = 3, 6, … , n enr for modes I, II, IV and V, as well as a = 0 , = 1, 2, 4, 5, … , n enr − 2, n enr − 1 for modes III and VI.
At the nodes of the remaining four cube faces, the displacement of the analytical solution is prescribed. Thus, in this area, the error will constitute the classical interpolation error. Consequently, optimal convergence properties will prevail there automatically, independent of an enrichment. This means that only the two-dimensional displacements at the cube faces where a crack front meets the discretized boundary need to be studied. For symmetry reasons, we only determined the error norms on one of these two cube faces (at x = −1). Detailed numerical experiments indicated that the results obtained using the decomposition enrichment method hardly differ from those obtained using the direct enrichment method. This is why we choose the more efficient direct method for further considerations.
Bilinear shape functions on discretized boundary (C4). Figure 7 shows the convergence of the displacement solutions with increasing mesh density in a double-logarithmic plot. Particularly, the results for the six considered crack modes I-VI in the L 2 -norm and energy norm are depicted. The corresponding numerical calculations were performed using bilinear shape functions (C4-elements, p = 1). Black data points represent results from computations without enrichment and red data points those with an enrichment using the singular crack modes with = 0.5. The straight lines again originate from regression calculation with exponential approximation functions. For this regression, the two finest discretizations were employed. In both norms, a significant improvement of both the accuracy and rate of convergence can be observed due to the enrichment (see also Table 2) for the modes I-III (solid lines). Without enrichment, the rate of convergence for all three crack modes is m st L 2 ,I-III ≈ 1 in the L 2 -norm and m st a,I-III ≈ 0.5 in the energy norm. With enrichment, the rates of convergence improve to m enr L 2 ,I-III ≈ 2 in the L 2 -norm and to m enr a,I-III ≈ 1.5 in the energy norm. These values correspond to the optimal values which would also be expected from the error estimates of a standard FEM with and without enrichment (cf Table 1). However, the rate of convergence of the relative error in the energy norm even exceeds the expected value of 1. For the crack modes IV-VI, which contain the weaker singular factor r 3/2 (dashed lines in Figure 7 Table 1.

Biquadratic shape functions on discretized boundary (C8)
As we have seen, the rates of convergence expected for the classical FEM (Table 1) coincide with those obtained from numerical calculations with the SBFEM using bilinear shape functions, at least for the canonical example described above. Now, the question arises whether this can also be observed when using higher order shape functions. To answer this question, biquadratic shape functions at the discretized boundary are considered. The respective SBFEs contain eight nodes ‖ at the boundary and we denote them C8-elements in the following. Analogously to Figure 7 for C4-elements, Figure 8 shows the results of the convergence study for C8-elements. Figures 8A and 8C show the convergence of the relative error in L 2 -norm respectively energy norm with the same enrichment which was used for the C4-elements (crack modes with = 0.5). For crack modes I-III (solid lines), the accuracy significantly improves in the L 2 -norm as well as in the energy norm due to the enrichment. Since the crack modes I-III do not contain any further singularity than the one associated to = 0.5, the expected orders  ||u − u h ||, measured in || · || L 2 ( Ω) -and || · || a( Ω) -norm for C8-elements (biquadratic Lagrange shape functions in 1 , 2 ) without enrichment (SBFEM) and with enrichment (enrSBFEM, enriched with two-dimensional crack modes 2D = 0. 5, 1.5) of convergence with enrichment are 3 for the L 2 -norm and 2 for the energy norm, already when using an enrichment with deformation modes up to = 0.5. Indeed, when using this enrichment, the rates of convergence can be significantly improved from m st a,I-III ≈ 2.8 in the energy norm. These convergence orders are even higher than the expected order of 3, respectively, 2 and remain unchanged for a further enrichment with crack modes up to = 1.5 (see also Figure 8B and 8D). In contrast to that, the deformation modes IV-VI (dashed lines) do contain a singularity associated to = 1.5, but none associated to = 0.5. Consequently, the expected rates of convergence without enrichment are 2.5 in the L 2 -norm and 1.5 in the energy norm. Indeed, with m st L 2 ,IV-VI ≈ 2.4 and m st a,IV-VI ≈ 1.5 (see Table 3), higher rates of convergence than for deformation modes I-III without enrichment are achieved. An enrichment with crack modes with displacement exponent = 0.5, here, does not lead to the desired improvement. On the other hand, an enrichment up to a displacement exponent of = 1.5 yields the rates of convergence of m enr 1.5 L 2 ,IV-VI ≈ 3.6, respectively m enr 1.5 a,IV-VI ≈ 2.6 (cf Table 3 and Figure 8B and 8D). Thus, with an appropriate enrichment, the same orders of convergence as they are from the known error estimates of the standard FEM are obtained.

Distribution of the relative error
In order to gain a deeper insight into the impact of enrichment on the displacement solution, we consider the error distribution as well as the local convergence in certain patches at the discretized boundary. In the following, we denote the cube faces containing points where the crack front meets the discretized boundary (here at x = ±1) Γ s and the remaining cube faces Γ ns , so that Ω = Γ s ∪ Γ ns . Figure 9 shows the distribution of the relative error measured in the L 2 -norm (A,B) and energy norm (C,D), without (A,C) and with (B,D) enrichment, when considering the example of the BVP for crack mode I. This example uses a discretization with 14 C4-elements per cube edge (ie, bilinear shape functions) and, in this case sufficient, the enrichment with 2D crack modes up to = 0.5. Even without enrichment, the L 2 -error at the boundary Γ ns already is very small ( Figure 9A). This is since the analytical solution is directly prescribed at the boundary nodes here making the error simply represent the interpolation error. Because  of the fact that at Γ ns only the smooth angular function of the crack mode has to be reproduced with bilinear shape functions, this error is small in comparison. In contrast, at the boundary Γ s , the error is expectedly large when no enrichment is considered since there the √ r -function with an infinite gradient at the crack front has to be reproduced. With enrichment, the error at Γ s decreases to approximately the level of the interpolation error at Γ ns , which means a drastic reduction. When considering the distribution of the relative error measured in the energy norm ( Figure 9C,D), it becomes apparent yet from the different scale that the error in the energy norm exceeds the error in the L 2 -norm by more than an order of magnitude. Without enrichment, the energy error mainly concentrates in the vicinity of the crack faces, ie, also at Γ ns where only the interpolation error exists ( Figure 9A). This results from the fact that the stress vector at the crack faces reduces to exactly zero in the analytical solution, which is difficult to reproduce for the approximate solution. The resulting error in the stresses enters the integrand of the energy norm quadratically. If there is an enrichment with the exact 2D solution ( Figure 9D), obviously, the stress boundary conditions are fulfiled exactly, so that the relative error, particularly in fully enriched elements, is drastically reduced. In blending elements, the error is reduced so that it approximately cor-

FIGURE 10
Patches on the discretized boundary of the considered subdomain Ω: numbered for Table 4     solving full boundary value problem (C4-elements) to represent crack mode I (cf Figure 10) responds to the interpolation error. The energy error close to the crack faces at Γ ns (interpolation error) naturally remains unchanged as there is no enrichment in effect.

Local convergence.
For the analysis of local convergence, patches E1, … ,E10 were defined ( Figure 10) and the discretizations with K = 4, 8, 12, 16, 20 elements per cube edge were considered. The discretizations with C4-elements are chosen such that the nodes of the coarsest mesh also are nodes of the finer meshes. Figures 11A and 11B show the convergence of the relative error of the displacements in the L 2 -norm for these patches at the example of crack mode I. Figure 11A contains the results for the standard formulation of the SBFEM. As expected, the error is particularly high in the patches E1 and E2, which are located directly adjacent to the crack front, whereas the error in patches E9 and E10 describing the interpolation error is particularly small. The corresponding rates of convergence are m st L2 ≈ 1, respectively, m st L2 ≈ 2 ( Table 4). The rate of convergence as well as the accuracy for patches E3-E8 moves between these two extremes but also shows that the negative influence of the singularity is not limited to its immediate vicinity. This indicates that the use of simple quarter point elements might not be sufficient here to recover the optimal rate of convergence and that an enrichment of a constant sized domain is required. Figure 11B shows the respective results of the enrSBFEM. Here, the relative error of all considered patches converges with the same optimal order. Figures 11C and 11D contain the results for the energy norm. Without enrichment, the relative error is particularly high in patches E1 and E2 close to the singularity. In the further patches E3 and E10, which are located in the vicinity of the crack faces, the error also is comparably high. Nevertheless, the error only converges in patches E1 and E2 with the reduced order (h 0.5 ) and otherwise with (h). With enrichment, the accuracy again is significantly improved, especially in patches E1 and E2, where the error now even converges quadratically and is more than two orders of magnitude smaller for all discretizations. However, the rate of convergence is also consistently improved on patches E3-E8 to a value of about m enr a ≈ 1.5. Only on patches E9 and E10 the error still converges with the expected rate of the interpolation error and, consequently, remains unchanged.

CONCLUSIONS
The SBFEM is already established as an appropriate method for 2D crack problems when the crack tip is entirely located within the considered domain. However, in 3D, the crack front describes a line so that there still exists a singularity at the point where the crack front intersects the discretized boundary. The impact of an additional use of enrichment (enrSBFEM, newly introduced in the work of Hell and Becker 79 ) on the convergence properties of deformation modes has been studied and discussed for the example of a single straight crack in a homogeneous isotropic continuum. The results have been compared to the known improvements on convergence rates achieved for an enrichment of a standard FEM for a problem including a singularity. According to numerical investigations, the enrichment of the SBFEM with singular asymptotic near fields seems to have a similar effect on the convergence properties as a corresponding enrichment in a standard FEM, such that the accuracy and rates of convergence could be retrieved. In contrast to previous works, besides the asymptotic near-field functions with = 0.5 also those with = 1.5, 2.5, … were used for the enrichment, which is required to recover the optimal rates of convergence when Lagrange shape functions of higher order than one are used.