An enriched Scaled Boundary Finite Element Method for 3D cracks

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

An enriched Scaled Boundary Finite Element Method for 3D cracks Sascha Hell , Wilfried Becker ⁎

Technische Universitaet Darmstadt, Franziska-Braun-Strasse 7, 64287 Darmstadt, Germany

ARTICLE INFO

ABSTRACT

2018 MSC: 00-01 99-00

In this work, a direct enrichment approach for the Scaled Boundary Finite Element Method (SBFEM) is proposed and implemented. While the SBFEM has proven its effectiveness in the treatment of 2D boundary value problems with stress singularities, it suffers from reduced convergence rates and accuracy in 3D cases with line singularities. This is due to the then unpreventable occurrence of singularities in the discretized boundary coordinates of the SBFEM. It is shown that enriching the characteristic separation of variables approach of the method directly with the 2D asymptotic near-fields of the line singularities recovers the optimal convergence rates in a uniquely DOF-efficient way. It will also be seen that, at the same time, the presented enrichment approach for the SBFEM can yield a very high accuracy – even for rather coarse discretizations.

Keywords: SBFEM enrSBFEM Three-dimensional stress singularity Enrichment 3D crack problem Linear elasticity

1. Introduction The efficient and accurate solution of boundary value problems (BVPs) which contain singularities is especially demanding [1–4]. Singularities occur at points or along lines which typically represent some kind of discontinuity, like discontinuities of geometry (e.g. cracks and notches [5–7]), material (e.g. multi-material junctions [8–12]) or loading. These singularities include theoretically infinite values in the derivatives of the basic problem variables. In this work, we focus on problems of linear elasticity where the basic physical quantities are the displacements. The singularities in the first derivative accordingly occur for mechanical strains and stresses. But other physical relations as e.g. heat conduction or even coupled physics like piezoelectric problems could equally be addressed [4,13–15]. 1.1. Prerequisites We distinguish between line singularities (2D stress singularities) and point singularities (3D stress singularities) [16–20]. Line singularities describe singular stresses along a line (e.g. notch or crack front, bimaterial edge, etc.) and are characterized by a twodimensional decay behaviour perpendicular to this line of singular stresses. On the other hand, point singularities e.g. occur at discontinuities of line singularities, like at kinks of cracks [21–24] or cracks meeting a domain’s boundary [24–30], and include a particular three-dimensional decay of the stresses (cf. Fig. 1). This decay with increasing distance from the singularity typically is of a 1, so that the resulting asymptotic nearpower-law type and particularly controlled by the so-called stress singularity exponent field can respectively be written as:

⁎

Corresponding author. E-mail address: [email protected] (S. Hell).

https://doi.org/10.1016/j.engfracmech.2019.04.032 Received 17 March 2019; Accepted 26 April 2019 0013-7944/ © 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Sascha Hell and Wilfried Becker, Engineering Fracture Mechanics, https://doi.org/10.1016/j.engfracmech.2019.04.032

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Nomenclature

ne nS N N p p, p , p PsS

Abbreviations COS coordinate system DES differential equation system DOF degrees of freedom (enr)SBFEM (Enriched) Scaled Boundary Finite Element Method (G)SIF (generalised) stress intensity factor LES linear equation system MWF modal weighting factors SBCOS scaled boundary coordinate system

P r R S t u u u Wa, Wi x , y, z x

Indices and operators operator for virtual quantities Euclidean vector norm value at outer/inner boundary quantity associated to positive/negative part of the spectrum of eigenvalues (·)2D , (·)3D quantity associated to 2D/3D situation (·) ext extrapolated value from Richardson extrapolation (·)st , (·) enr associated to standard formulation / enrichment (·) derivative in scaling coordinate

· (·)a , (·) i (·) p , (·) n

Greek symbols

u,

Latin symbols

E Ei F G h I J Jr KI, K L

m

t

a,

i

s,

ns

R

derivative from matrix of shape functions vector of free constants in SBFEM solution, MWFs, GSIFs stiffness matrix relative error in the approximation of the displacement exponents / SBFEM eigenvalues Young’s modulus system matrices in SBFEM equations (i = 0, 1, 2 ) matrix of enrichment functions shear modulus characteristic element length identity matrix Jacobian of SBCOS at = 1 Jacobian of conical COS at = 1 KII, KIII stress intensity factors of classical crack modes matrix of coefficients of SBFEM-DES differential operator edge length of considered cubical domain convergence rates of displacement exponents

B c, c u , c C e

number of elements number of scaling rays with singular stress matrix of shape functions set of scaling rays order of Lagrange shape function load vector point, where scaling ray with singularity meets boundary set of scaling rays in enriched domain radial coordinate in polar, conical or spherical COS radial coordinate in cylindrical COS scaling center surface tractions (boundary condition) displacement vector approximate solution for displacements displacement functions in along scaling rays external/internal work Cartesian coordinates position vector

,

1,

uj ,

j

u

2

j

boundary of domain boundary with displacement/traction boundary condition outer/inner boundary discretized boundary (with/without singularity) boundary along scaling rays (e.g. crack faces) strain tensor (in vector notation) approximated strains angular coordinate in spherical COS material parameter displacement exponent / eigenvalue diagonal matrix of displacement exponents / eigenvalues Poisson’s ratio coordinates in SBCOS rescaled radial coordinate stress tensor (in vector notation) approximated stresses matrix of angular functions of deformation modes eigenvector from solution of SBFEM-DES matrix of eigenvectors j displacement part in matrix of eigenvectors angular coordinate in polar, conical or spherical COS

Fig. 1. Exemplary geometrical configurations with stress singularities. 2D: (a) V-notch with notch opening-angle (for = 0 : crack situation); 3D: (b) Fichera problem with three V-notch edges meeting in one point, (c) wedge-shaped crack with a kink in the course of the crack front, (d) surfacebreaking crack with a crack front ending at a traction-free boundary. 2

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

u r,

c u2D j r

=

2D j

2D uj (

),

r,

j =1

u r, ,

c 2D j r

=

2D 1 2D ( j

j

),

(1)

j=1

c u3D j r

=

3D j

3D uj

,

,

j =1

r, ,

c 3D j r

= j=1

3D 1 3D j

j

,

.

(2)

The asymptotic near-field with respect to a point Pasymp (Fig. 1) fulfils the field equations as well as the local boundary conditions in the vicinity of Pasymp . In problems of essentially two-dimensional nature, the asymptotic near-field is given in polar coordinates r, . In inherently three-dimensional problems with a point singularity, spherical coordinates r, , are used. While the stress singularity exponents j 1 describe the decay behavior of the stresses from Pasymp , the exponents j describe the increase of displacements from Pasymp . For a clear distinction between them, we will refer to j as the displacement exponents and j 1 as the stress exponents. The angular tensor functions uj and j are denominated deformation modes and stress modes respectively, as they mainly describe the shape of the displacement and stress field of each term j in the power-law series. Like the displacement exponents, deformation and stress modes only depend on the local boundary conditions. Compatibility with far-field boundary conditions is ensured by the modal weighting factors c uj (MWF) or the generalized stress intensity factors c j (GSIFs). The MWFs and GSIFs are linearly dependent and therefore equivalently describe the intensity of each term j in the asymptotic near-field. To ensure purely real-valued displacements and stresses, all three factors within a term in the power-law series are either real-valued or occur in complex conjugated pairs. Of course, the 2D problems generally result from a simplification of the actual 3D case and are therefore associated with line singularities with a local cylindrical r , , s coordinate system whose polar plane is always perpendicular to the line of discontinuity. Then, the MWFs and GSIFs become a function of s and can vary along the line of discontinuity. 1.2. Methods for the analysis of stress singularities As there is no general analytical procedure to treat 3D BVPs with point and line singularities, numerical approaches like finite differences methods (FDM) and finite element methods (FEM) need to be employed. Unfortunately, numerical methods face special challenges in the presence of singularities (e.g. [1,31–35]). They struggle with representing the very high (up to infinite) gradients in the vicinity of a singularity. This is because they are bound to a finite discretization and can therefore only resolve finite gradients. High differences in the density of the discretization within a considered domain can also lead to numerical difficulties (large condition numbers of the involved matrices). All of this results in a generally reduced accuracy as well as convergence of the BVP solution. Several strategies have been developed to face this challenge. Some are:

• The employment of a domain rescaling in the vicinity of the singularity (graded meshes/discretization [1,36,37]; rescaling within • • •

one row of finite elements [3,38–41]). It involves a rescaling of the vicinity of the singularity such that the discontinuous singular 1) , the associated displacement field behaves like near-field is transformed into a smooth one: if the stress exponent is Re( u r Re( ) , and with a rescaling of the radial coordinate = r Re( ) the singularity disappears in the scaled coordinates: u . The use of displacement representations which are enriched with the asymptotic near-fields of the singularity (enriched FEM [42–49], eXtended FEM a.k.a. XFEM [50–54]). This approach has the advantage that modal weighting factors (or GSIFs) can be directly included as degrees of freedom (DOF) of the actual finite element analysis and is generally very DOF-efficient. A major issue of this approach is the accuracy of the numerical integration within an enriched element as the integrand includes singular terms from the enrichment. The exclusion of the vicinity of the singularity from the discretized domain, using a certain set of analytically determined deformation modes near the singularity and finally, coupling the MWFs to the discretized domain (singular hybrid finite elements [3,55,56]). The analytical consideration of the BVP in a radial direction towards the singularity while a discretization is only performed in the circumferential direction (FEM eigenanalysis techniques [24,57–67]; Scaled Boundary FEM a.k.a. SBFEM [20,68–75]).

Except for the strategy mentioned last, all other methods require the a priori-knowledge of the existing asymptotic near-fields (1) and (2) in some way. Generally, this makes the FEM eigenanalysis techniques, like the SBFEM, indispensable for the accurate consideration of 3D stress singularities. FEM eigenanalysis techniques use a separation of variables approach for the displacements ui in the BVP of linear elasticity (here 3D):

ui = r

i(

, ),

(3)

i = 1, 2, 3.

This ansatz strongly resembles the terms in the power-law series of the asymptotic near-fields which is no coincidence: This method was developed in [24] to be able to determine the asymptotic near-fields at point singularities in 3D elasticity problems, as analytical methods seemed to mostly fail in this case. However, it was not until the (independent) development of the Scaled Boundary Finite Element Method (SBFEM, [68–70,76–78]) that the field of application of these techniques was extended from the determination of asymptotic near-fields to the solution of full BVPs which must also include the fulfilment of far field boundary conditions and thereby the determination of MWFs. In addition to the advantage of being independent of the supply of asymptotic 3

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

near-fields from other sources, the major advantage of the SBFEM is that it avoids a discretization in the radial direction, where the high gradients arise, and confines the discretization to the circumferential coordinates. In doing so, the convergence challenge is simply and completely circumvented – at least in the case of two-dimensional problems where the SBFEM was successfully applied to fracture mechanics problems ranging from the analysis of asymptotic near-fields to the evaluation of GSIFs (e.g. [79–81]) and even automatic crack growth algorithms (e.g. [82–85]). A late comprehensive overview of fracture mechanics applications is given in [75]. The application of the method to other physics is e.g. presented in [13,14,86–91]. 1.3. The SBFEM for the analysis of 3D stress singularities Despite the mentioned advantages of the SBFEM and its successes, it is to be pointed out that in three-dimensional problems where point and line singularities arise, singularities inevitably also exist in the discretized circumferential coordinates. The consequences are again reduced convergence properties and accuracy. The proposal of a solution to this specific problem is going to be the focus of this work. The first successful attempt to address this problem was delivered in [26] using a domain scaling approach – although not in conjunction with the application of a finite element discretization but a finite difference scheme. Bažant and Estenssoro [24] and Somaratna and Ting [92] tried to beneficially adapt the number of Gauss points while using a spherical coordinate system with its pole at the line singularity. But this approach was not generally expedient and also quite limited with respect to the location of the line singularities. The works [93–101] employed a simple mesh refinement towards the line singularities in the circumferential coordinates using biquadratic shape functions and an isoparametric element formulation. In Ref. [102] an adaptive meshing algorithm based on an error indicator that works with globally smoothed stresses is employed. Apel et al. [19] were the first to apply the domain rescaling approach in the FEM eigenanalysis technique. In comparison to [93,94], they employed a mesh refinement following the domain rescaling scheme such that the convergence orders of a smooth problem were obtained again for displacement exponents as well as deformation modes. In particular, they also gave detailed mathematical error estimates to their graded mesh approach. However, while such graded meshes yield the desired convergence rates, they generally are not very efficient regarding the necessary degrees of freedom. The reason why a DOF-efficient approach is important, is that the solution process of the SBFEM and FEM eigenanalysis techniques includes an eigenvalue problem. The computational effort for the solution of an eigenvalue problem grows cubically with the number of degrees of freedom. Consequently, while the considered 3D domain only needs to be discretized in the circumferential coordinates reducing the problem dimension by one, this advantage can be consumed by the quickly growing effort for solving the eigenvalue problem if the DOF are not confined to a minimum. This is especially true when the complete set of eigenvalues and eigenvectors, representing the displacement exponents and deformation modes, need to be determined which becomes necessary when full BVPs are to be solved. Then, more efficient partial eigenvalue solvers like iterative Arnoldi methods [19,103] cannot be used anymore. In reference [104], a quarter-point element formulation (according to [39]) was applied in a FEM eigenanalysis to solve the problem of a surface-breaking crack in an isotropic body. This approach involved the placement of the midside-nodes of the underlying biquadratic Lagrange shape functions (in the boundary coordinates 1, 2 ) to the 1/4-position. A similar rescaling of the element domain has also been employed by [105] but in a case where the displacement exponent of the line singularity was = 0.708 0.5. Consequently, they placed the midside-nodes at the 3/8-position instead of the 1/4-position leading to an element domain rescaling exponent ln(1/2)/ln(3/8) 0.7067 ( 0.708) . The authors claim having obtained good results although not giving results of a thorough convergence analysis. While the fact that the displacement exponent does not exactly meet the rescaling exponent seems to be rather negligible, their approach also includes abandoning the completeness of the employed shape functions. This is because a constant stress distribution can no longer be exactly represented within such rescaled elements. Additionally, the domain on which the line singularity is represented appropriately by the shape functions is confined to the single row of elements adjacent to the singularity. Consequently, the size of the rescaled domain varies with element size, thereby somewhat limiting the convergence order. Of these mentioned strategies, so far, only general mesh refinement schemes as employed by Gharemani [93,94] have been used to solve full BVPs, as e.g. by [20]. Other works, as [73], do not especially consider the impact of line singularities on accuracy and convergence properties. But particularly in these cases, more DOF-efficient approaches are in demand. The adaptation of the enrichment approach1 (enrichment of the displacement representation with the appropriate 2D near-fields of present line singularities) seems to be promising due to the following points: 1. Properly formulated, only very few additional DOF need to be introduced. 2. The additional computational cost due to the numerical integration of singular functions will quickly be outweighed by the savings in DOF.

1 Note that such a procedure is to be explicitly distinguished from the ones described by Natarajan et al. [106,107] to which they refer as eXtended Scaled Boundary Finite Element Method (xSBFEM). Instead of including an enrichment of the actual displacement representation in the SBFEM, their xSBFEM is based on a coupling of subdomains where a 2D XFEM is employed and subdomains where a 2D SBFEM is employed. The role of the SBFEM in their procedure is to cover the singular behaviour at the crack tip.

4

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

3. As the enrichment includes the exact 2D near-fields, results should be very accurate already for a rather coarse discretization in the circumferential coordinates. First efforts in this direction have been made in [108,109] 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. Using a decomposition of the appropriate analytic asymptotic near-field and a simple higher-order Gauss integration scheme within the enriched elements, the authors were able to confirm a significantly improved performance of the enriched formulation. The present work ties in with these forgoing ones. The method is detailed further, its efficiency improved and extended to easily include arbitrary 2D line singularities different from the classical crack singularity: This different, more DOF-efficient, direct formulation that directly uses the 2D asymptotic near-fields of the involved line singularities, is used in conjunction with an improved numerical integration scheme originally proposed in [110] (for an XFEM framework). The direct formulation of the enrSBFEM is described in detail in Section 2 after the basics of the standard 3D SBFEM have been briefly revisited. The convergence properties of the enriched formulation are studied and compared to the standard formulation in Section 3. Two major examples are considered here: an essentially two-dimensional case of a single crack in a homogeneous isotropic continuum, for which an analytical solution exists, and an inherently three-dimensional case of two meeting cracks in a homogeneous isotropic continuum. Finally, some conclusions on the enrSBFEM and its performance are drawn in Section 4. 2. Theory 2.1. Scaled Boundary Finite Element Method (SBFEM) The formulation of the Scaled Boundary Finite Element Method in the framework of linear elasticity theory is briefly revisited. We consider a domain with two similar, scalable boundaries: an outer boundary a and an inner boundary i (cf. Fig. 2). Similar and scalable means that the inner boundary i conforms to a scaled version of the outer boundary a . This scaling happens with respect to a scaling center S by a factor of i/ a < 1. The scaling center S must be located such that any straight line (scaling ray) drawn from S would cross each of the two boundaries i and a once. This is the scalability requirement. This definition explicitly includes that scaling rays themselves may present a part R of the overall boundary , so that

=

=

a

(4)

R.

i

For convenience and according to the scalability requirement, a scaled boundary coordinate system with coordinates , in 2D or , 1, 2 in 3D is introduced (Fig. 2). The scaling coordinate is running from the scaling center to the boundaries (along the scaling rays with a value of i at the inner boundary i and a value of a at the outer boundary a ). The boundary coordinates 1, 2 are running along the boundary (for 2D domains, just one boundary coordinate is sufficient). In this work, we concentrate on 3D problems of linear elasticity and the following derivations are given correspondingly. The coordinate transformation from a Cartesian coordinate system x , y , z into the scaled boundary coordinate system is given by

x = x 0 + x ( 1,

(5)

2)

with x 0 being the Cartesian coordinates of the Scaling Center. The associated Jacobian J ( 1, x

y

z

x

y

z

x 1

=

1

x 2

1

y 2

1

y

2

z

2)

is

x

=J

1,

2

y

z 2

(6)

z

Fig. 2. Scaled boundary coordinate system in (a) 2D and (b) 3D domain. 5

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

and the volume differential can be rewritten as

dV = d

= det(J )

2d

(7)

d 1 d 2.

In linear elasticity with surface tractions t on a boundary sector t (volume forces neglected) and displacement boundary conditions u on a boundary sector u , the principle of virtual work holds (e.g. [34,69,70,111]). Its displacement-based formulation makes use of the differential operator x

LT

0

0

=

0

y

0

0

0

0

z

z

y

z

0 x

y x

0

(8)

and the elasticity matrix C and can be written as

uTLTC L u d

Wi =

=

t

uTt d = Wa.

(9)

Herein, u is the displacement vector and u is the vector of virtual displacements. An approximating separation of variables approach is employed for the displacements u as well as the virtual displacements u

u

u( ,

1,

2)

= N ( 1,

2) u (

),

u

u( ,

1,

2)

= N ( 1,

2)

(10)

u ( ).

The dependence on the scaling coordinate is confined to the vector functions u( ) and u( ) . The dependence on the boundary coordinates 1, 2 is confined to the approximating shape functions N( 1, 2 ) making u( ) and u( ) shape function coefficients which can vary along the scaling rays. Consequently, the discretization of the domain and the approximation of the field functions is limited to the boundary coordinates 1, 2 (cf. Figs. 3 and 4). As approximating shape functions, Lagrange polynomials are usually employed (e.g. [20,112]) but also other basis functions like Gauss-Lobatto polynomials and Fourier series could be used (e.g. [64,71,113]). Eq. (10) are employed in the virtual work Eq. (9). Expansion of the products, partial integration in to eliminate the derivative of the virtual displacements and integration over the boundary coordinates 1, 2 leads to the following representation a i

[ u ( )]T (E0 2 u ( ) + [2E 0 + E1T

E2 ] u ( ))d + [[ u ( )]T ( 2E 0 u ( ) + E1T u ( ))] a

E1] u ( ) + [E1T

i

[ u ( a )]T pa (11)

[ u ( i )]T pi = 0.

For brevity, boundary sectors R are considered to be unloaded and unconstrained and (·) denotes the derivative in the scaling coordinate . The system matrices E 0 , E1, E 2 and the boundary load vectors pa and pi at the outer and inner boundary respectively are defined by

E0 = E1 = E2 =

pa =

2 a

1

2

1

2

1

2

1

BT CB det(J ) d 1 d 2 ,

(12)

BT CB det(J ) d 1 d 2 ,

(13)

BT CB det(J ) d 1 d 2 ,

2

NT ta

x

× 1

x 2

(14)

d 1d 2,

Fig. 3. Examplary SBFEM discretization of a 2D domain scaling ray boundary R .

(15)

, which fulfills the scalability requirement, with inner boundary i , outer boundary

6

a

and

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Fig. 4. Examplary SBFEM discretization of a 3D rectangular solid with admissible choice of the scaling center S and only an outer boundary

pi =

2 i

1

x

NT ti

2

x

× 1

d 1 d 2.

2

a.

(16)

Herein, t a, t i denote the surface tractions at the outer and inner boundary respectively. The matrices B and B result from the kinematic equations for the strains = L u :

,

1,

=LN

2

1,

u( )

+ L

2

N ( 1, 1

2)

+L

1

N ( 1, 2

2)

2

u( )

= B

1,

u( ) 2

+ B

1,

u( ) 2

.

(17)

Requiring that this equation holds for arbitrary virtual displacements u( ) , the following two SBFEM equations can be derived:

E 0 2 u ( ) + [2E0 + E1T 2 a E0 u ( a ) + 2 i E0 u ( i )

E1] u ( ) + [E1T

T a E1 u ( a ) T i E1 u ( i )

=

(18)

E 2 ] u ( ) = 0,

pa pi .

(19)

Eq. (18) is a homogeneous Cauchy-Euler differential equation system (DES) of second order. It corresponds to the requirement of fulfilling the field equations of linear elasticity (static equilibrium, kinematics, material law) and the boundary conditions at the boundary sectors R (Fig. 2). The linear equation system (19) (LES) enforces the boundary conditions at the outer and inner boundary a and i respectively. The DES (18) can be transformed into a system of ordinary differential equations of first order with constant coefficients. This is done by a variable substitution t = ln( ) and the introduction of auxiliary functions v ( ) = u ( ) . The ordinary DES can be converted into an eigenvalue problem

e j t (K

j j)

j

= 0.

(20)

with

E 0 1 [E 0

K=

E1 + E1T] I

E0 1 [E1T 0

E2 ]

(21)

and leads to the DES solution u

u( ) =

(22)

c.

Here, is a diagonal matrix consisting of the eigenvalues j with j = 1, …, n and is a diagonal matrix consisting of the terms j . The matrix u consists of the displacement part of the eigenvectors j . Vector c consists of constants which are to be determined from LES (19). For 3D problems, the eigenvalue spectrum is symmetric to = 1/2 (for 2D problems, symmetric to = 0 ) and can be 1, k = 1, …, n /2 . Accordingly, the whole DES divided into two parts: a positive part with pk 0 and a negative part with nk solution can also be separated into two parts: u p

u( ) =

pc

p

+

u n

nc

(23)

n.

This representation is very similar to the asymptotic near fields in the vicinity of stress singularities (cf. Eqs. (1) and (2)). Accordingly, the partial eigenvectors in the matrices up , un can similarly be interpreted as deformation modes from a mechanical point of view. The eigenvalues j are decay exponents which govern the decay of deformation modes from the boundaries ( up decay from the outer boundary, un decay from the inner boundary) and we denote them displacement exponents (in contrast to the stress exponents 1). Finally, the constants in the vectors c p, c n are modal weighting factors (MWF) of the positive and the negative deformation modes. As already mentioned, they are determined from the LES (1) and (2))

(E 0 (E 0

u p

p u p

+ E1T p

+

u p+ I p) a +I E1T pu ) i p

(E 0 (E 0

u n

n u n

+ E1T n

+

E1T

u n+I n) a u n+ I n) i

p cp = pa cn i 7

(24)

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

and related to the boundary displacements by

u ( a) ua = = ui u ( i)

u p p a

u n n a

u p i

u n n i

p

cp . cn

(25)

2.2. Direct Enrichment of the SBFEM (enrSBFEM) The SBFEM formulation is extended to incorporate an enrichment (enriched SBFEM; abbreviation: enrSBFEM). Different from the description in [108] a direct enrichment is presented, in which the separation of variables ansatz (cf. Eq. (10)) is directly enriched with the asymptotic near-field solutions of the line singularities. In the following, the differences, which arise from the enrichment of the separation of variables ansatz (10), to the standard SBFEM formulation are highlighted. 2.2.1. Scaled Boundary Coordinate Systems (SBCOS) In addition to the scaled boundary coordinate systems , 1e , 2e (SBCOS, Eq. (26)) of the standard formulation, one for every element e, further SBCOS with boundary coordinates 1s , 2s are to be defined, one for each enrichment s = 1, …, nS , i.e. for each line singularity along a scaling ray (Fig. 5). As the asymptotic near-fields of line singularities are typically given in 2D polar coordinates, it is useful to also introduce an auxiliary polar COS with radial coordinate rs and angular coordinate s in the 1s , 2s -planes. Together with the scaling coordinate , the radial and angular coordinates rs and s form a conical COS. Accordingly, the coordinate transformations of the enrSBFEM to the Cartesian x , y , z -COS is:

x = x 0 + x (e) ( 1e ,

x = x 0 + x (rs) (rs,

e 2)

s)

(26)

SBCOS,

(27)

conicalCOS.

The Jacobian Jr of the conical COS is given by x

y

z

x r

y r

z r

x

y

z

x r

=

r

r

r

y

x

= Jr

y

z

r

. (28)

z

The associated differential operator L becomes

L = Lx

x

+ Ly

y

+ Lz

z

=L

+ Lr

1 r

1 r

+L

(29)

and the volumetric and surface differential can be written as:

dV =

Jr r ,

2r

d dr d ,

dS =

xr r

×

xr r

2r

dr d .

(30)

A simple example is given in Appendix A, for clarity.

Fig. 5. (a) A plane, isoparametric scaled boundary finite element represents a pyramidal subdomain. The scaling center S is at = 0 (here located at the origin of the Cartesian x , y, z -COS). The elemental SBCOS , 1e , 2e and a conical COS , rs, s are given with rs = 0 along the line singularity s (here: crack front). Scaling rays meet the boundary at the (four) nodes of the discretization of the boundary coordinates. (b) Boundary mesh of a cubical domain with crack. The scaling center S is located at the center of the cubical domain (on the crack front). The crack (crack faces lightly shaded, implemented by double nodes) meets the discretized boundary at the points PsS (s = 1, …, nS with nS = 2 ). The base side of elements containing enriched nodes (highlighted by red circles) form the partial boundary s (remaining boundary denominated ns ). The element of (a) can be found in (b), darkly shaded and rotated by 90° around vertical axis. 8

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

2.2.2. Enrichment of the Separation of Variables Ansatz The fundamental precondition for the applicability of the SBFEM is that the considered BVP fulfills the scalability requirement (Section 2.1), so that a separation of variables approach can be employed. This separation of variables approach is the basis for the semi-analytical character of the method. Accordingly, an enrichment of the displacement field also needs to comply with this approach. For the displacements in an enriched scaled boundary finite element e, we consequently write:

u (e )

,

e 1,

e 2

= N ( 1e ,

e 2 ) ue (

) + F (s ) (

e 1,

e 2D(s ) ( 2)c

). (31)

e e uenr s ( , 1, 2)

e e ust e ( , 1, 2)

, The first summand represents the standard SBFEM separation of variables representation with the displacements along the scaling rays in the vector function u e ( ), and the matrix of shape functions N( 1e , 2e ) to interpolate between them. The u est (

e 1,

e 2)

S second summand u enr s ( , 1, 2 ) is the enrichment associated to the line singularity along the scaling ray SPs . This enrichment consists e e ( s ) of the matrix of enrichment functions F ( 1 , 2 ) in the discretized boundary coordinates and the associated coefficient functions in the vector c 2D(s) ( ) . In this work, only asymptotic near-fields of line singularities are considered as enrichment functions. While in [108,109,114] decompositions of such asymptotic near-fields were employed, here, the asymptotic near-fields are directly used (without being preliminarily decomposed), which we denote direct enrichment method, or in short direct method (in contrast to the formerly presented decomposition method). Consequently, in the direct method, the vector c 2D(s) ( ) of coefficient functions, which represent the additional 2D DOFs due to the enrichment, can directly represent the MWFs c u2D j or GSIFs c j of the line singularities (Eq. (1)) or at least a value directly proportional to them. Compared to the decomposition method, the direct method is somewhat simpler in its implementation. But its major advantages are that

(1) it is more DOF-efficient than using the decomposed near-fields because only one additional DOF per enrichment and term in the asymptotic near-field is needed, and (2) it facilitates the possibility of considering more than only one singularity exponent for the enrichment, which becomes numerically unstable in the decomposition method. The enrichment u enr (Eq. (31)) of the SBFEM must fully comply with the requirements of any Finite Element Method to its s displacement ansatz functions (cf. e.g. [42,115] for the enrichment of a standard FEM). For brevity, the superscripts will be neglected in the following illustrations. According to the interpolation condition, the functions u( ) should keep on representing the displacements along the scaling rays (at the nodes in the discretized boundary coordinates). For the enriched separation of variables approach, this is only possible if the enrichment vanishes at the scaling rays (at the nodes in the discretized boundary coordinates). Accordingly, the values of the enrichment functions along the scaling rays are multiplied with the shape functions N( 1, 2 ) and subtracted from the enrichment functions. Futhermore, typically not the whole domain is to be enriched with the same enrichment functions (see e.g. the cubical cracked domain in Fig. 5b), but is confined to the close vicinity of the line singularity. To ensure the compatibility of enriched with nonenriched parts of the domain, the enrichment needs to be reduced to zero. In the XFEM, this is called blending and can simply be done by multiplying the enrichment with shape functions N( 1, 2 ), then denominated blending functions, which take the value of 1 at enriched nodes and zero at non-enriched nodes. This procedure can be transferred to the enrSBFEM, where the blending functions need to take the value of 1 at the enriched scaling rays and 0 at non-enriched scaling rays. With these two modifications, the “Partition of Unity” condition is also fulfilled and there is no problem with the completeness of the shape functions (as e.g. occurs for quadratic elements with midside nodes located at a position different from the 1/4 or 1/2 position, e.g. [105]). The complete enriched separation of variables representation for the displacements ui (x), i = x , y, z , can then be written as nS

ui (x) =

Nk (x ) uik ( ) k N

+

nF

Nl (x ) l P

2D Nk (x ) Fism (xr k) csm ( ) .

Fism (xr ) s=1 m =1

k N

(32)

The first summand of Eq. (32) still conforms to the standard SBFEM formulation and the second summand is the enrichment. Set N is the set of all scaling rays (or all nodes of the discretization in the boundary coordinates 1, 2 ). Set P is the set of enriched scaling rays. N ( 1, 2) are the approximating shape functions and Fism (r , ) are the enrichment functions containing the considered 2D n terms of the asymptotic near-fields r j u2D j ( ) of the line singularities (Eq. (1)). The number of enrichment functions is given by F and nS is the number of scaling rays with line singularities. The enrichment can only be non-zero between the scaling rays because the interpolation Nk (x ) Fism (xr k) of the enrichment functions is subtracted from the enrichment functions Fism (r , ) themselves. Ac2D ( ) are the 2D MWFs along cordingly, the functions uik ( ) indeed represent the displacements along the scaling rays. The functions csm the scaling rays. The enrichment is reduced to zero from the enriched scaling rays to the non-enriched ones using the blending functions Nl . 2.3. Further implementation: numerical integration, assembly, solution Just as in the standard formulation of the SBFEM, the separation of variables representation for the displacements is to be 9

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Fig. 6. Quasi polar integration: decomposition of the quadrilateral domain into two triangular domains; then, mapping of a quadrilateral domain and its Gauss points to each of the two triangular domains (by collapsing one side of the quadrilateral), so that the Gauss points near r = 0 (red point) move closer together. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

employed in the principle of virtual work (9). For the virtual displacements a similar separation of variables representation is to be used, but with virtual displacement and virtual coefficient functions u( ) and c 2D ( ) along the scaling rays. For brevity of the further notation, the free functions and their virtual couterparts are combined to

u( ) , c 2D ( )

uc ( ) =

uc ( ) =

u( ) . c 2D ( )

(33)

Analogous to the standard formulation, in the equation of strains, those matrix functions in the enrichment which do not contain derivatives in the boundary coordinates are combined to a matrix function Benr ( 1, 2) , and those with derivatives in the boundary coordinates are combined to a matrix function Benr ( 1, 2) (for details, see Appendix B, Eq. (B.2)). Consequently, the integration over boundary coordinates 1, 2 and the scaling coordinate can still be executed separately. st st Integration over the boundary coordinates results in the well-known system matrices Est 0 , E1 , E 2 from the standard part of the enr enr enr separation of variables representation (31), and further system matrices E0 , E1 , E2 from the enrichment and coupling system matrices Ec0 , E1c1, E1c2 , E2c which contain coupling terms between the standard and the enrichment part. Accordingly, compared to the standard formulation of the SBFEM, seven further system matrices are to be determined for the enriched formulation. Moreover, the integrands of these additional system matrices include singular values due to the singularity in the derivatives of the asymptotic nearfields, which constitutes a certain challenge for an accurate numerical integration. The numerical integration over the boundary coordinates is performed in an element-wise manner using a high-order Gauss integration rule. It is sufficient to just use the SBCOS , 1, 2 , even if the asymptotic near-fields for the enrichment are given in polar coordinates, because the coordinates in a polar/conical COS can easily be converted into those of an SBCOS. However, for the integration of the singular derivatives of the enrichment functions in elements directly adjacent to the line singularity, it seems convenient to use quasi-polar coordinates in the corresponding boundary coordinate plane (Fig. 6, cf. [110]). To this end, the typically rectangular integration domain of such an element is divided into two triangular ones and to each of them an isoparametric square integration domain is mapped in such a way, that the Gauss points near the singularity become very close. This approach makes use of the fact that the surface differential r dr d of a polar coordinate system contains the radial coordinate r as a factor, which potentially mitigates a singular integrand.2. Thus, the otherwise required computational effort for an accurate numerical integration can be reduced significantly. Nevertheless, the necessary integration order still depends on the strength of the singularity in the derivatives of the enrichment functions. For all the other elements, a high-order Gauss integration rule usually gives sufficiently accurate results. In the current implementation, an integration order of nG = 10 is employed (exactly integrating polynomials up to order 2nG 1 = 19).3 The system matrices are assembled to

E0 =

E st0

E c0

(E c0) T

Eenr 0

, E1 =

E1st

E1c1

(E1c2) T

E1enr

, E2 =

E st2

E2c

(E c2) T

E2enr

.

(34)

They are to be employed in the 3D SBFEM equations (here: no volume forces, no loads along scaling rays):

u c ( )(E0

[ u c ( )(E0

2u

c(

2u

c(

) + [2E0

E1 + E1T ] u c ( ) + [E1T

) + E1T u c ( ))] = =

a i

E 2 ] u c ( )) = 0,

= u c ( a ) pa + u c ( i ) pi.

(35) (36)

As before, (·) denotes the derivative in the scaling coordinate . Regarding the boundary load vectors pa , pi of the outer and inner boundary respectively, it is to be considered that, for their determination, the respective integrands now also contain the enriched separation of variables representation of the virtual displacements. Accordingly, also here a higher order Gauss integration is recommendable for the terms containing the enrichment functions. Eq. (35) must hold for arbitrary virtual displacement functions u c ( ) . Consequently, the solution of the resulting DES continues to be obtained by solving an eigenvalue problem. It is to be noted, that every scaling ray (node in the discretization of the boundary coordinates) implies 6 DOF in the eigenvalue 2 3

1/ r ) the integrand even becomes simply constant in the radial coordinate r. E.g. in the case of the classical crack singularity ( For comparison, in [48] a Gauss integration order of nG = 48 was employed for an enriched 3D FEM considering an interface crack. 10

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

problem (the number of dimensions times the order of highest derivative in the DES). For the current enriched formulation (direct method), only 2 additional DOF per enrichment function, i.e. per enriching 2D deformation mode along a singularity line, need to be introduced. In the case of the classical crack singularity, 3 deformation modes with singular stresses are to be expected, which accordingly can be covered by only 6 additional DOFs per crack tip in the discretized boundary coordinates, i.e. only 12 additional DOF for the example in Fig. 5b, independent of the discretization in the boundary coordinates 1, 2 . Thus, the enrichment hardly adds any computational effort at the stage of the eigenvalue problem solution, while it can be expected to considerably enhance the accuracy. The major additional computational cost comes from the numerical integration of the enrichment functions and their derivatives. But as already mentioned, a quasi-polar integration scheme reduces the required effort significantly. Additionally, the numerical integration procedure can be massively parallelized reducing the computation time further. In conclusion, when the discretization is refined, the savings in computation time in the eigenvalue problem solution stage quickly outweigh the additional computational cost from the numerical integration. The solution of the eigenvalue problem, with eigenvalues j and the associated eigenvectors T T

T u j

st j

=

enr j

,

u

respectively

= [(

st ) T ,

(

enr ) T]T ,

(37)

also contains parts from the standard formulation, nodal displacement values and parts from the enrichment, enrichment coefficients jenr . With this, the solution of the displacements of the enrSBFEM results in: st j ,

n

ui (x) =

Nk (x )

st ikj

n jc

+

j

j=1 k N

nS

nF

Nl (x ) j=1 l P

Fism (xr ) s=1 m=1

Nk (x ) Fism (xr k) k N

enr smj

jc

j

.

(38)

The weighting factors cj for the corresponding deformation modes follow from Eq. (36) and the resulting4 LES as well as from the boundary conditions at i (inner boundary) and a (outer boundary). 3. Results and discussion In this section, the convergence properties of the direct enrichment method for the SBFEM will be studied. Two examples are considered (Fig. 7): (1) a straight plane crack in a homogeneous isotropic continuum, (2) two perpendicularly meeting straight plane cracks in a homogeneous isotropic continuum. The first example basically is a two-dimensional one. Here, an analytical solution exists for the associated asymptotic near-fields, which can be used for comparison. The second example is inherently three-dimensional. Accordingly, there is no analytical solution. In fact, its associated asymptotic near-fields can only be determined using FEM eigenanalysis methods like the SBFEM. The focus of this study is on the comparison of the enriched and the standard formulation of the SBFEM – especially regarding the loss of optimal convergence of the latter in the presence of line singularities and their recovery using an appropriate enrichment. For the evaluation of the convergence properties, in particular the displacement exponents j , which have been extracted from the SBFEM eigenvalue problem, are considered. If the analytical solution is known, the relative error of the calculated displacement exponents j compared to the analytical ones j is given as:

e

,j

=

j

m

j

C hm = C

j

ne

.

(39)

For a sufficiently fine discretization, a convergence behaviour following a power law with constant factor C and convergence order m can be expected, as for any Finite Element Method; h = / ne constitutes a characteristic element dimension in the discretized 1, 2 plane. However, in the general 3D case, the analytical displacements j are not known, so that the error measure e , j cannot be determined directly. Nevertheless, by means of a Richardson-Extrapolation [116,117] of the type j

= C hm +

ext j ,

(40)

an extrapolated displacement exponent can be calculated from the eigenvalue problem solutions of 3 different discretizations. This extrapolated displacement exponent can be used as an approximate exact solution, which becomes more and more accurate with increasing refinement of the employed 3 discretizations. To the best knowledge of the authors, mathematical convergence estimates of FEM eigenanalyses like the SBFEM, so far, only exist for the standard formulation with uniform and graded meshes [19]. However, the convergence properties found by Apel et al. [19] are very similar to those of a classical FEM (cf. e.g. [33–36]). Accordingly, we expect that also the results for convergence properties ext j

4 Eq. (36) must be fulfilled for arbitrary virtual displacements / enrichment coefficients u c ( a ), u c ( i ) at the outer and inner boundary respectively.

11

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Fig. 7. Cubical subdomain with scaling centre located at the centre of the subdomain: (a,c) with one crack, crack flanks shaded (b,d) with two perpendicularly meeting cracks; (c,d) with discretized boundary and enriched nodes marked in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

of an enriched classical FEM (e.g. [45]) can be transferred to an enriched formulation of the SBFEM and we will assess this hypothesis here. In the classical Finite Element Method (displacement-based, isoparametric formulation, Lagrange shape functions) the expected order of convergence m for a smooth boundary value problem (i.e. without singularities) depends on the polynomial order p of the shape functions; displacements converge with the order m = p + 1, stresses and strains with m = p and strain energy density with m = 2p . For boundary value problems which include singularities, the order of convergence depends on the lowest displacement exponent min of all active deformation modes (except for rigid body motions). Consequently, in the case of the classical crack singularity, displacements only converge with the order of about m = 2 min = 1.0 , stresses and strains with m = min = 0.5 and strain energy density with m = 2 min = 1.0 . However, in the standard FEM, using an appropriate enrichment restores the original convergence rates of a smooth boundary value problem. In the case of the classical crack singularity, the appropriate enrichment functions would be the 2D asymptotic near-fields of the classical crack modes I, II and III. We shall see in the course of this section that this also seems to be true in the case of the SBFEM. 3.1. Numerical model From a three-dimensional homogeneous isotropic continuum, which contains at least one crack, a cubic subdomain with edge length = 2 mm is considered and it is modelled using the SBFEM. Since the resulting stress fields are independent of Young’s modulus E in this case and the associated displacement fields simply scale inversely proportionally with E (cf. Eqs. (42)–(44)), the results shown in this section are practically independent from the material properties. For Poisson’s ratio, a value of = 0.3 was chosen, here, but this also has a minor effect on the results. A crack front runs from one center of a square face of to the center of the opposing face. The crack faces are presumed to be free (no loading, no displacement constraints) and parallel to a square face of . Fig. 7 shows such a subdomain for a continuum with a single straight crack (Fig. 7a) and for a continuum with two perpendicularly meeting cracks (Fig. 7b). This consideration of a cracked subdomain in a continuum ensures that no boundary of the subdomain is a real, for example traction-free boundary of a structural component, as this would usually imply a distinct 3D stress singularity where the crack would meet the boundary of the structure (cf. surface-breaking crack [24–26,118,119]). However, this is not the structural situation of interest here. Moreover and more importantly, we would like the discretized 1, 2 -plane to only contain 2D singularities or no singularities and to be free of 3D singularities. To properly represent 3D stress singularities with the SBFEM, they should be located at 12

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

the scaling center. The faces of subdomain are described by the boundary coordinates 1, 2 and uniformly discretized with an even number of ne quadrilaterals along each edge length of the cubic subdomain. Then, the cracks can be modelled by double nodes in the boundary discretization, yielding

nDOF =

(6ne2 + 2) + (2ne

1) ·3

(41)

crack

cube

DOF if no enrichment is taken into account. The scaling center is located at the center of the cubic subdomain, such that the crack front is represented by two scaling rays. Linear Lagrange shape functions are employed. In the enriched formulation, all the scaling rays ending at the face of a subdomain where a crack meets a discretized boundary are enriched, except for the scaling rays ending at the edges of the cubic subdomain. This ensures that, independent from the discretization, always an equally sized domain ist enriched. According to [110], this was important for fully recovering optimal convergence in the XFEM and we followed this recommendation. As mentioned in Section 2.3, for the numerical integration of the SBFEM system matrices, a Gauss quadrature of order 10 is employed in enriched elements. At the same time, a quasi-polar integration is performed in the elements directly adjacent to the crack front (Fig. 6). To additionally save computation time, the numerical integration is parallelized. The 3 singular terms of the 2D asymptotic near-field at the crack front in a homogeneous isotropic continuum are employed as enrichment, here. They are e.g. given in [120]. The displacements of the crack modes I, II and III (the 3 deformation modes with displacement exponent = 0.5) are:

uI r ,

uII r ,

uIII r ,

=

KI 2G

r 2

KII 2G

r 2

=

=

2KIII G

cos( )

cos( /2) sin( /2) 0

mode I,

sin( /2)[2 + + cos( )] cos( /2)[2 cos( )] 0 r 2

0 0 sin( /2)

(42)

mode II, (43)

mode III (44)

with = 3–4 under the assumptions of plane strain (the assumptions of plane stress do not yield a solution which fulfills all field equations of 3D linear elasticity and therefore should not be used for an enrichment in a truely 3D boundary value problem).

Fig. 8. Deformation modes for a single crack in a homogeneous isotropic 3D-continuum. 1st row: classical crack modes with 3D = 2D = 0.5 (e.g. [120]). 2nd row: higher crack modes with 3D = 2D = 1.5 and equally 2D displacement field. 3rd row: higher crack modes with a displacement field that is variable in the crack front direction, 3D = 1.5 but 2D = 0.5.

13

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

All computations have been performed in

MATLAB®

using the authors’ own 3D SBFEM implementations.

3.2. Straight crack in a homogeneous isotropic continuum The canonical case of a single straight plane crack in a homogneous isotropic continuum is considered. The relative errors between the eigenvalues j , j = I, …, IX determined by solving the eigenvalue problem (20) and the analytical values I,II,III = 0.5 and IV,V,VI = 1.5 are studied. Here, we identify the 2D deformation modes with = 1.5 with IV, V and VI. They are illustrated in Fig. 8, 2nd row. The 3rd row of this figure shows modes which exhibit a displacement field, which varies along the crack front and for which no analytical solution exists. For these crack modes, in contrast to crack modes I-VI, the associated 2D deformation mode (line 3D 2D = 1.5). The interaction of 2D and 3D = 0.5) is not equal to the 3D deformation mode (point singularity, VII,VIII,IX singularity, I,II,III stress singularities in the enrSBFEM is discussed more deeply in Appendix C. Please note that above and in the following, if 2D = 3D , the index ·2D , respectively ·3D is often omitted for brevity. According to Apel et al. [19], in absence of stress singularities in any discretized plane, an optimal rate of convergence of m = 2p can be expected, where p is the polynomial degree of the employed shape functions. But if there is a 2D-deformation mode with 2D within a discretized plane, the rate of convergence is, independent of the polynomial degree p, limited to m = 2 2D . Fig. 9a and b show the convergence of the relative error of the calculated eigenvalues for crack modes I-IX in double logarithmic plots. Bilinear shape functions ( p = 1) at the boundary were employed in the corresponding calculations. Circular marks describe the results obtained using the enriched formulation of the SBFEM (enrSBFEM, direct method), squares those from a standard SBFEM. It is to be mentioned that, although the same discretizations are considered, computations using the enriched formulation include 12 additional DOF compared to computations with the standard SBFEM. However, as the total number of DOF quickly grows with increasing ne (nDOF (n e = 6) = 687), the effect of the additional 12 DOF also quickly fades out. The lines in Fig. 9a and b result from a

1 , where j were determined by the standard respectively enriched SBFEM using Fig. 9. Comparison of the convergence of the relative error j / j bilinear shape functions. a) Modes I (black), II (blue), III (green); b) Modes IV and VII (black), V and VIII (blue), VI and IX (green). Lines are derived from a regression function of the type Cne m using the two finest discretizations to determine C and m . Solid lines: 2D modes; dashed lines: 3D modes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

14

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

regression calculation using a function Cne m where the two free constants are determined by the two finest discretizations. In the double logarithmic plot, the gradient of these regression lines represents the convergence rate m . By means of the regression lines, it can be argued from a phenomenological point of view whether a further refinement of the discretization might lead to a further change of the apparent convergence rate. With and without enrichment, the relative error with respect to the analytical eigenvalues converges to zero, but there are significant differences in the convergence rate as well as the accuracy in general. Furthermore, as illustrated in Fig. 9a, the calculated error values are spread and get closer to the regression line for finer discretizations. We observe this especially for crack mode II (marked blue). This behaviour has also been observed in Hell and Becker [108], where a simple high order Gauss integration was employed to compute the system matrices and the QR-algorithm and the decomposition method instead of the direct enrichment method was used. In that work, only calculations with mesh sizes up to ne = 20 elements per edge were performed, since the numerical effort by using the QR-algorithm, which determines all eigenvalues and all eigenvectors, quickly becomes rather large. However, as only a few eigenvalues are needed in a simple convergence study of the eigenvalues, we employed the MATLAB®-function eigs() here, which is based on iterative Arnoldi methods [103], to approximate several eigenvalues near the value 0 and the respective eigenvectors. Due to the reduced necessary computational effort, we could apply up to the mentioned ne = 32 elements per edge of the cubic subdomain. In Fig. 9a, crack modes I-III are considered. The significant improvement of accuracy due the enrichment is obvious as the relative error is reduced by several orders of magnitude and, in fact, already is very small for coarse discretizations. Without enrichment, the convergence rates are limited to double the value of the associated displacement exponent 2D = 0.5 (m stI III 2 I III = 1). With enrichment, also the convergence rates are significantly improved (cf. Table 1). For crack mode III, even the convergence rate mopt = 2 expected in a smooth BVP is surpassed (m enr 2.4 ). For crack mode I and II “only” convergence rates of m enr 1.9 and III I m enr 1.7 respectively are reached. But especially for crack mode II, the relative position of the associated regression line to the data II points indicates a probably increasing convergence rate at further refined discretizations. In Fig. 9b crack modes IV-IX ( 3D = 1.5) are considered. The regression lines associated to crack modes VII-IX are dashed, the three regression lines associated to crack modes IV-VI are represented by solid lines. For crack modes IV-VI, the results obtained without and with enrichment are almost identical, giving a consistent convergence rate of mopt 2 . This is because the displacement exponent of the 2D crack mode being active on the discretized boundary of the considered subdomain is not smaller than the polynomial degree of the employed shape functions ( 2D = 1.5 > p = 1). In contrast, crack modes VII-IX exhibit a displacement field which is variable in the crack front direction. Different from the so far considered crack modes I-VI, here, the displacement exponent of the 2D crack mode being active on the discretized boundary of the considered subdomain is not equal to the one of the 3D crack mode anymore ( 3D = 1.5), but only constitutes 2D = 0.5 < p = 1. Consequently, without enrichment, only a convergence rate of mst 1 and a reduced accuracy is found, just as previously with crack modes I-III. With enrichment, the relative error of the calculated 3D displacement exponents is again similarly decreased by more than an order of magnitude and convergence rates of menr 2 are achieved (cf. Table 1). 3.3. Two perpendicularly meeting cracks in a homogeneous isotropic continuum The so far considered case of a single straight crack in an infinite homogeneous isotropic continuum was indeed posed in threedimensional space but actually is of two-dimensional character. This made it possible to use analytical solutions for the evaluation of the convergence properties of the standard SBFEM and the proposed enriched formulation. Now, an inherently three-dimensional crack interaction problem, which was considered by Hell and Becker [99] using a standard SBFEM formulation with a graded discretization, shall be revisited: the case of two perpendicularly meeting straight cracks (in a homogeneous isotropic continuum; cf. Fig. 7b). For this case, no analytical solution is known. Accordingly, instead of an analytical solution, an extrapolated value from a Richardson-Extrapolation (40) is employed as a reference. This structural situation of two perpendicularly meeting cracks includes six deformation modes leading to singular stresses at the crack interaction point (where the two cracks meet): the two symmetric crack opening modes co1/co2, the two crack shearing modes cs1/cs2 and the two crack twisting modes ct1/ct2 (cf. Fig. 11). In Fig. 10a and b the convergence of the displacement exponents j is illustrated. Fig. 10a shows the displacement exponents calculated by the SBFEM with (disks) and without (squares) enrichment for discretizations with ne = 2, 4, …, 16, 20, …, 32 elements along each edge of the cubical domain. The extrapolated displacement exponents ext determined from the Richardson-Extrapolation (40) are very similar for the SBFEM formulations with and without enrichment (Table 2 with a maximum relative deviation <0.1%). They are depicted in Fig. 10a as Table 1 Obtained rates of convergence m of the relative error

without (SBFEM) and with (enrSBFEM) enrichment (with Crack mode SBFEM enrSBFEM

j/ j

1 of the eigenvalues/displacement exponents

j

= 0.5).

using bilinear shape functions

I

II

III

IV

V

VI

VII

VIII

IX

1.0 1.9

1.0 1.7

1.0 2.4

2.0 2.0

2.0 2.0

2.0 2.0

1.0 3.9

1.0 2.0

1.0 2.1

15

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Fig. 10. Convergence of displacement exponents/eigenvalues j of 3D deformation modes co1/co2 (black, solid/dashed), cs1/cs2 (blue), ct1/ct2 (green) for homogeneous isotropic continuum. Comparison between standard formulation (SBFEM) and enriched formulation (enrSBFEM) using bilinear shape functions. Lines are (a) extrapolated displacement exponents ext form a Richardson-Extrapolation (40) or (b) regression lines with base functions Cne m which are determined using the two finest discretizations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. 3D deformation modes of asymptotic near-field of two perpendicularly meeting cracks with displacement exponents

16

j

(Table 2).

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Table 2 Extrapolated displacement exponents ext from a Richardson-Extrapolation (40) and convergence rates m of the associated relative error for deformation modes with 0 < j < 1 using bilinear shape functions (p = 1) without (SBFEM) and with enrichment (enrSBFEM, 2D = 0.5). deformation mode SBFEM, ext enrSBFEM, SBFEM, mst

co1

co2

cs1/cs2

ct1

ct2

0.3824

0.6731

0.5077

0.3388

0.7919

1.99

1.93

1.93

2.09

2.20

0.3824 0.99

ext

enrSBFEM, menr

0.6730 1.13

0.5077 1.00

0.3385 1.09

0.7921 0.95

Table 3 Extrapolated displacement exponents ext from a Richardson-Extrapolation and convergence rates m of the associated relative error for deformation modes with j > 1 using bilinear shape functions ( p = 1) without (SBFEM) and with enrichment (enrSBFEM, 2D = 0.5). Deformation mode SBFEM, ext enrSBFEM, SBFEM, mst

ext

enrSBFEM, menr

14

15

16–19

20–21

22

23

24

– 1.260

1.298 1.298

1.412 ± 0.150i 1.411 ± 0.150i

1.499 1.500

1.696 1.696

1.795 1.797

2.000 2.000

1.99

2.12

2.24

1.85

2.00

1.99

1.98

–

1.09

1.06

1.36

1.54

0.79

2.01

horizontal lines. This diagram shows a clearly improved convergence of the displacement exponents when using the enriched formulation of the SBFEM: While the values obtained using the standard formulation only slowly approach the extrapolated values, the deviation of values obtained using the enriched formulation already is negligibly small for very coarse discretizations. The double logarithmic representation in Fig. 10b shows this even clearer. For the enriched formulation, the displacement exponents converge quadratically towards the extrapolated values ext , and for the standard formulation linearly (see also Table 2). At the same time, the accuracy of the enrSBFEM at a discretization with ne = 4 is not obtained by using the standard formulation until a discretization with ne = 32 . Table 3 gives the results obtained for higher deformation modes, i.e. such deformation modes with displacement exponents j of higher values. Deformation modes with j = 0 and j = 1 (rigid body motions and deformations with homogeneous strains) are not considered, as they are already reproduced exactly in an SBFEM with bilinear shape functions (C4 elements), independent of an enrichment. Here, the extrapolated displacement rates of an enriched formulation again match those obtained using the standard formulation quite well as the relative deviation remains smaller than 0.1% . 5, Also, the enriched formulation continues yielding quadratic convergence. Only for the deformation modes j = 20, 21 the convergence rate is somewhat lower (menr 1.85) but still clearly outperforms the convergence rate obtained using the standard SBFEM formulation. As expected, only for the deformation mode with 24 = 2 quadratic convergence is also achieved for the standard formulation. Accordingly, the enriched SBFEM has been shown to retrieve optimal convergence properties also for three-dimensional boundary value problems with inherently three-dimensional stress singularities. 4. Conclusions The Scaled Boundary Finite Element Method using enriched base functions (enrSBFEM), newly introduced in Hell and Becker [108], has been extended and improved to directly include 2D asymptotic near-fields of line singularities (direct enrichment method). This results in

• a simpler implementation, • an increase of efficiency, • a more concise behaviour when using higher order Lagrange shape functions, • and an extended range of problems which can be considered. The accuracy as well as the convergence properties of displacement exponents being associated with 3D singular asymptotic near-fields have been studied using the examples of a single straight crack and two meeting cracks both in a homogeneous isotropic continuum. It has been shown in numerical experiments that the enrichment with singular asymptotic near-fields in the Scaled Boundary Finite Element Method has a comparable effect on the accuracy and convergence properties like a corresponding 5 An exception only occurs for 14 , for which no consistent extrapolated values could be obtained when using the standard formulation. The deviation with respect to the extrapolated value obtained from the enriched formulation was indeed small ( 0.1%), but failed to show monotonic convergence.

17

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

enrichment in a 2D standard FEM. The implementation of the newly proposed direct enrichment method for the enrSBFEM is more straight-forward than the former decomposition enrichment method [108,109,114] because the analytically determined asymptotic near-fields can be employed directly, without further manipulation. We find an increase of efficiency because instead of the formerly necessary 12 additional DOF, now only about 3 DOF are sufficient per line singularity (2D singularity in the discretized boundary coordinates) in the considered subdomain. Furthermore, the consideration of asymptotic near-field functions with several different displacement exponents can be expected to become numerically stable and therefore feasible in the first place. This is the reason why, besides the asymptotic near-field functions with = 0.5 also those with = 1.5, 2.5, … will become available for the enrichment, which is required to recover the optimal rates of convergence when Lagrange shape functions of higher order than linear ones are employed. This will also extend the range of treatable elasticity problems with 3D stress singularities (point singularities) to include arbitrarily anisotropic material behaviour and multimaterial problems. Concluding, the proposed direct enrichment method for the SBFEM results in a very efficient, accurate and versatile tool for the analysis of linear elasticity problems with point singularities and line singularities. Acknowledgements The authors highly appreciate the financial support of the German Research Foundation (DFG), project BE1090/35-2. Appendix A. A conical coordinate system for the enrSBFEM A comprehensible example is given for a conical coordinate system as it is employed for the use in the enrSBFEM. In accordance with Fig. 5, the coordinate transformation from the conical to a general Cartesian coordinate system is given by

d r sin( ) , r cos( )

x0 y0 + z0

x y = z

(A.1)

with d denoting the distance of the boundary inverse J r 1 are

d = 0 0

Jr r ,

r sin( ) r cos( ) sin( ) cos( ) , cos( ) sin( )

s

at r = 0 from the scaling center S. The corresponding Jacobian matrix Jr and its

Jr 1 r,

1/ d 0 0

=

r /d sin( ) cos( )

0 cos( ) . sin( )

(A.2)

The decomposition of the differential operator L finally results in

L = L x /d

+

Lx r /d

L ysin( ) + Lzcos( )

L

1 = Lx d

1

Lr

r

1 r

+ Ly

sin( )

1 r

1 cos( ) r

r

+

L ycos( )

Lzsin( )

L

+ Lz cos( )

1 r

1 r

sin( )

1 r

.

(A.3)

Appendix B. Strains and stresses in the enrSBFEM The approximated strains result from the application of the differential operator L to the separation of variables representation u( , 1 , 2) for the displacements. The approximated stresses are calculated from the multiplication with the elasticity matrix C :

= L u,

=C

= C L u.

The enrichment in the separation of variables representation (32) of the displacements practically consists of two summands and each is a product of three functions. Accordingly, the product rule of differential calculus has to be applied to each of these two summands. Using the decomposition (29) of the differential operator L the mechanical strains in vector notation are given by the following representation:

18

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

= L N ( 1,

^( ) u

2)

Bst ( 1,

+ L

N ( 1, 2 ) 1

1

2)

+ L NB ( 1,

(

2 )(F (r ,

(L

NB ( 1, 2 ) 1

)

Nenr ( 1,

+L

1

NB ( 1,

2

2

(

2)

2)

1

F (r , ) r

Benr r ( 1,

= Bst

1,

+ Benr

2,

u( ) 2

1,

2,

Nenr ( 1, 2 )

Benr 2 ( 1,

+ L r NB ( 1,

2

+ Bst c2D ( )

2

) (F (r,

r ( 1,

)

2 ),

k ))

c2D ( )

2,

r ( 1,

2),

1,

( 1,

Nenr ( 1, ( 1,

+ L 2 NB ( 1 ,

+ L NB ( 1,

r ( 1,

2 ) F (rk ,

2)

NB ( 1, 2 )

Benr 1 ( 1, 1

2

Bst ( 1, 2 )

Benr ( 1 ,

+ L

^( ) u

N ( 1, 2 )

+L

2)

2),

( 1,

2)

F (r , ) r

2 ) F (rk ,

2 )) Nenr ( 1, 2 ) 2

)

)F (r , k

k ))

k)

c 2D ( )

c2D ( )

2 )) c 2D ( )

.

2 ))

(B.1)

u( ) 2

+ Benr

1,

c2D ( ) 2

(B.2)

The matrix functions N ( 1, 2 ), NB ( 1, 2 ) and Nenr ( 1 , 2) contain the same shape functions N ( 1, 2) but generally are of different dimension as they fulfil different tasks. E.g., the functions NB ( 1 , 2) ensure the “blending” between enriched and non-enriched subdomains. Accordingly, they only lead to gradients unequal to zero in the “blending elements” (elements which contain enriched and non-enriched scaling rays / nodes in the discretization of the boundary coordinates). I.e. the matrix Benr 1 can be neglected in all other elements. The matrix Benr r ( 1, 2 ) typically contains singular terms of the type r challenge in an accurate numerical integration (cf. Appendix 2.3).

2D

1

and, consequently, constitutes the major

Appendix C. Interaction of 2D and 3D stress singularities in the enrSBFEM The question arises, how 2D and 3D stress singularities interact in the framework of the enrSBFEM. This issue shall be briefly addressed, here, at a simplified example. From the displacement solution (38), it becomes immediately apparent that in both the standard and the enrichment part, the decay behavior of the deformation modes towards the scaling center ( = 0 ) is controlled by the displacement exponents (or 3D ). This is especially true for the coefficient functions n 2D csm ( )= j=1

enr smj

3D j c

j,

enr where smj determines (as part of the initially freely scalable eigenvector ju in representation (37)) the ratio at which the enrichment functions at = 1 are addressed6. At the same time, however, there is also a decay behavior of the enrichment in the radial direction r of the conical coordinate system, whereas the actual distance R from the scaling ray at r = 0 , which includes the line singularity, depends on both r and (cf. Fig. C12):

(C.1)

R = r.

At first glance, this results in an unclear situation, in which the 2D and 3D decay behaviour become difficult to separate. This is why a well-known and simplified case, in which only a 2D stress singularity exists, shall be considered here, so that this case could also be analysed by a pure 2D modelling. It is the situation shown in Fig. 7a,c of a straight, plane crack in the 3D continuum, where the cube shown represents only a subdomain which is cut out from the infinite continuum. The analytical solution in polar coordinates r, for the case of an infinite continuum and only one addressed deformation mode with decay exponent 2D is of the following type:

u(r , ) = c 2Dr

2D

2D (

(C.2)

).

This solution shall also be employed as enrichment in the separation of variables representation and the corresponding boundary value problem shall be solved. One way to do this is to prescribe the displacements of the analytical solution on the partial boundary ns , where no enrichment

6

Thus,

enr smj

contains, if directly enriched with the 2D deformation modes of the line singularities, the corresponding 2D GSIFs or 2D MWFs. 19

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

Fig. C12. Coordinates in cylindrical COS (R , , z ) and conical COS ( , r , ) at scaling center S with line singularity on

axis (R = r = 0 ).

takes place, and to lock the displacements perpendicular to the plane (partial) boundary s to create plane strain conditions. If the enrSBFEM displacement solution then converges towards the analytical solution, the displacements uik ( ) along the scaling rays in the enriched area also converge towards those of the analytical solution. Since the latter is just identical to the enrichment 2D ( ) along the scaling rays, in the enriched area between the scaling rays of the displacement solution, only the values Fism (rk , k ) csm following fraction remains to be approximately equalized with the given analytical field (C.2):

u ( , r , ) = F (r , ) c 2D ( ) = r

2D

2D (

)

3D

(C.3)

enrc 3D.

To compare the two fields, a transformation into a cylindrical COS (Fig. C12) is adequate. Then the displacement field (C.3) in the enriched area results in

( ) Rd z

u R, , z = = c 3D

enr

() z d

2D

3D

2D

2D (

R

2D

2D (

)

and consequently, for

3D

=

= c 2D (z ) R

2D

)

() z d

2D (

3D

enrc 3D

)

(C.4) (C.5)

2D

and c 3D

enr

= c 2D, correctly represents the given field (C.2) in the whole subdomain, as requested.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

Grisvard P. Elliptic problems in nonsmooth domains. Monographs and studies in mathematics vol. 24. Boston: Pitman Advanced Publishing Program; 1985. Maugin GA. The thermomechanics of plasticity and fracture. Cambridge University Press; 1992. Kuna M. Numerische Beanspruchungsanalyse von Rissen: Finite Elemente in der Bruchmechanik. Wiesbaden: Vieweg+Teubner; 2008. Yosibash Z. Singularities in elliptic boundary value problems and elasticity and their connection with failure initiation. Springer; 2012. Williams ML. Stress singularities resulting from various boundary conditions. J Appl Mech 1952;19(4):526–8. Karal FC, Karp SN. The elastic-field behavior in the neighborhood of a crack of arbitrary angle. Commun Pure Appl Math 1962;15(4):413–21. Knein M. Zur Theorie des Druckversuchs, in: Über die Grundlagen der Balkentheorie/Die Spannungen und Formänderungen von Balken mit rechteckigem Querschnitt/Stegbeanspruchung hoher Biegungsträger/Zur Theorie des Druckversuchs, Abhandlungen aus dem Aerodynamischen Institut an der Technischen Hochschule Aachen, vol. 7. Springer; 1927. p. 43–62. Karal F, Karp S. Stress behavior in the neighborhood of sharp corners. Geophysics 1964;29(3):360–9. Bogy DB. Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J Appl Mech 1968;35(3):460–6. Dundurs J. Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J Appl Mech 1969;36(4):650–2. Bogy DB. Two edge-bonded elastic wedges of different materials and wedge angles under surface tractions. J Appl Mech 1971;38(2):377–86. Bogy DB. On the plane elastostatic problem of a loaded crack terminating at a material interface. J Appl Mech 1971;38(4):911–8. Rubtsov YK. Scaled-boundary method in electroelasticity of thin plates. Int Appl Mech 2006;42(7):832–41. Rubtsov YK, Borisov E. Numerical-analytical solution of plane problems in thermoelasticity. Int Appl Mech 2007;43(12):1406–10. Mayland W. Untersuchungen zu Spannungssingularitätsordnungen in linear-elastischen und piezoelektrischen Multimaterialkonfigurationen mit der RandFinite-Elemente-Methode. In: Forschungsberichte des Instituts für Mechanik der Technischen Universität Darmstadt, Vol. 25, Studienbereich Mechanik, Technische Universität Darmstadt, Darmstadt, Germany; 2012. Sinclair GB. Stress singularities in classical elasticity–I: Removal, interpretation, and analysis. Appl Mech Rev 2004;57(4):251–98. Sinclair GB. Stress singularities in classical elasticity–II: Asymptotic identification. Appl Mech Rev 2004;57(5):385–439. Aksentian OK. Singularities of the stress-strain state of a plate in the neighborhood of an edge. J Appl Math Mech 1967;31(1):193–202. Apel T, Sändig A-M, Solov’ev SI. Computation of 3D vertex singularities for linear elasticity: error estimates for a finite element method on graded meshes, ESAIM. Math Model Numer Anal 2002;36(6):1043–70. Hell S, Becker W. The Scaled Boundary Finite Element Method for the analysis of 3D crack interaction. J Comput Sci 2015;9:76–81. Bažant ZP. Three-dimensional harmonic functions near termination or intersection of gradient singularity lines: a general numerical method. Int J Eng Sci 1974;12(3):221–43. Morrison JA, Lewis JA. Charge singularity at the corner of a flat plate. SIAM J Appl Math 1976;31(2):233–50. Parihar K, Keer L. Stress singularity at the corner of a wedge-shaped crack or inclusion. J Appl Mech 1978;45(4):791–6. Bažant ZP, Estenssoro LF. Surface singularity and crack propagation. Int J Solids Struct 1979;15(5):405–26. Benthem J. State of stress at the vertex of a quarter-infinite crack in a half-space. Int J Solids Struct 1977;13(5):479–92. Benthem J. The quarter-infinite crack in a half space; alternative and additional solutions. Int J Solids Struct 1980;16(2):119–30. Buchholz F-G, Chergui A, Richard H. Fracture analyses and experimental results of crack growth under general mixed mode loading conditions. Eng Fract Mech 2004;71(4–6):455–68. Heyder M, Kolk K, Kuhn G. Numerical and experimental investigations of the influence of corner singularities on 3D fatigue crack propagation. Eng Fract Mech 2005;72(13):2095–105. Heyder M, Kuhn G. 3D fatigue crack propagation: experimental studies. Int J Fatigue 2006;28(5–6):627–34.

20

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker

[30] Schnack E, Kolk A, Dimitrov K. 3D crack growth by considering re-entrant corners. In: ICCES: international conference on computational & experimental engineering and sciences, vol. 8; 2008. p. 107–14. [31] Johnson Jr. MW, McLay RW. Convergence of the finite element method in the theory of elasticity. J Appl Mech 1968;35(2):274–8. [32] Tong P, Pian THH. On the convergence of the finite element method for problems with singularity. Int J Solids Struct 1973;9(3):313–21. [33] Strang G, Fix GJ. An analysis of the finite element method vol. 212. Prentice-Hall Englewood Cliffs; 1973. [34] Bathe K-J. Finite-element-procedures. 2nd ed. Watertown, MA, USA: K.J. Bathe; 2004. [35] Zienkiewicz O, Taylor R, Zhu J. The finite element method: its basis and fundamentals. 7th ed. Butterworth-Heinemann; 2013. [36] Apel T, Sändig A-M, Whiteman JR. Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains. Math Meth Appl Sci 1996;19(1):63–85. [37] Apel T. Anisotropic Finite elements: local estimates and applications, advances in numerical mathematics. Vieweg+Teubner; 1999. [38] Akin JE. The generation of elements with singularities. Int J Numer Meth Eng 1976;10(6):1249–59. [39] Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. Int J Numer Meth Eng 1975;9(3):495–507. [40] Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. Int J Numer Meth Eng 1976;10(1):25–37. [41] Barsoum RS. Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements. Int J Numer Meth Eng 1977;11(1):85–98. [42] Benzley SE. Representation of singularities with isoparametric finite elements. Int J Numer Meth Eng 1974;8(3):537–45. [43] Gifford LNJ, Hilton PD. Stress intensity factors by enriched finite elements. Eng Fract Mech 1978;10(3):485–96. [44] Chen E. Finite element analysis of a bimaterial interface crack. Theor Appl Fract Mech 1985;3(3):257–62. [45] Blum H. Numerical treatment of corner and crack singularities. In: Stein E, Wendland W, editors. Finite element and boundary element techniques from mathematical and engineering point of view. CISM international centre for mechanical sciences, vol. 301. Springer; 1988. p. 171–212. [46] Kaya AC, Nied HF. Interface fracture analysis of bonded ceramic layers using enriched finite elements. In: Kokini K, editor. Proceedings of the symposium on ceramic coatings, ASME MD, vol. 44; 1993. p. 47–71. [47] Pageau SS, Biggers SB. Enrichment of finite elements with numerical solutions for singular stress fields. Int J Numer Meth Eng 1997;40(14):2693–713. [48] Ayhan AO, Kaya AC, Nied HF. Analysis of three-dimensional interface cracks using enriched finite elements. Int J Fract 2007;142(3):255–76. [49] Attaporn W, Koguchi H. Intensity of stress singularity at a vertex and along the free edges of the interface in 3D-dissimilar material joints using 3D-enriched FEM. Comput Model Eng Sci (CMES) 2009;39(3):237–62. [50] Melenk J, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Meth Appl Mech Eng 1996;139(1–4):289–314. [51] Fleming M, Chu YA, Moran B, Belytschko T, Lu YY, Gu L. Enriched element-free Galerkin methods for crack tip fields. Int J Numer Meth Eng 1997;40(8):1483–504. [52] Duarte C, Babuška I, Oden J. Generalized finite element methods for three-dimensional structural mechanics problems. Comput Struct 2000;77(2):215–32. [53] Fries T-P, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. Int J Numer Meth Eng 2010;84(3):253–304. [54] Mohammadi S. XFEM fracture analyses of composites. John Wiley & Sons; 2012. [55] Tong P, Pian THH, Lasry SJ. A hybrid-element approach to crack problems in plane elasticity. Int J Numer Meth Eng 1973;7(3):297–308. [56] Pian THH. State-of-the-art development of hybrid/mixed finite element method. Finite Elem Anal Des 1995;21(1–2):5–20. [57] Yamada Y, Okumura H. Analysis of local stress in composite materials by the 3D finite element. Proceedings of Japan-US conference on composite materials. 1981. p. 55–84. [58] Dasgupta G. A finite element formulation for unbounded homogeneous continua. J Appl Mech 1982;49(1):136–40. [59] Leguillon D, Sanchez-Palencia E. Computation of singular solutions in elliptic problems and elasticity. John Wiley & Sons; 1987. [60] Gu L, Belytschko T. A numerical study of stress singularities in a two-material wedge. Int J Solids Struct 1994;31(6):865–89. [61] Pageau SS, Biggers SB. A finite element approach to three-dimensional singular stress states in anisotropic multi-material wedges and junctions. Int J Solids Struct 1996;33(1):33–47. [62] Leguillon D, Sanchez-Palencia E. On 3D cracks intersecting a free surface in laminated composites. Int J Fract 1999;99(1–2):25–40. [63] Prukvilailert M, Koguchi H. Stress singularity analysis around the singular point on the stress singularity line in three-dimensional joints. Int J Solids Struct 2005;42(11–12):3059–74. [64] Korepanova T, Matveenko V, Sevodina N. Numerical analysis of stress singularity at singular points of three-dimensional elastic bodies. Acta Mech 2013;224(9):2045–63. [65] Khaji N, Khodakarami M. A semi-analytical method with a system of decoupled ordinary differential equations for three-dimensional elastostatic problems. Int J Solids Struct 2012;49(18):2528–46. [66] Nagai M, Ikeda T, Miyazaki N. Stress intensity factor analysis of a three-dimensional interface crack between dissimilar anisotropic materials under thermal stress. Eng Fract Mech 2012;91:14–36. [67] Ikeda T, Koga Y. Scalar parameter analysis of three-dimensional interfacial corner of jointed dissimilar anisotropic materials. In: Papadrakakis M, Papadopoulos V, Stefanou G, Plevris V, editors. Proceedings of the 7th European congress on computational methods in applied sciences and engineering (ECCOMAS 2016); 2016. [68] Song C, Wolf JP. The scaled boundary finite-element method-alias consistent infinitesimal finite-element cell method-for elastodynamics. Comput Meth Appl Mech Eng 1997;147(3–4):329–55. [69] Deeks AJ, Wolf JP. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Comput Mech 2002;28(6):489–504. [70] Wolf JP. The scaled boundary finite element method. John Wiley & Sons; 2003. [71] He Y, Yang H, Deeks AJ. Use of Fourier shape functions in the scaled boundary method. Eng Anal Bound Elem 2014;41:152–9. [72] Behnke R, Mundil M, Birk C, Kaliske M. A physically and geometrically nonlinear scaled-boundary-based finite element formulation for fracture in elastomers. Int J Numer Meth Eng 2014;99(13):966–99. [73] Saputra AA, Birk C, Song C. Computation of three-dimensional fracture parameters at interface cracks and notches by the Scaled Boundary Finite Element Method. Eng Fract Mech 2015;148:213–42. [74] Klinkel S, Chen L, Simeon B. Solids in boundary representation–a NURBS based Galerkin formulation. PAMM 2016;16(1):211–2. [75] Song C, Ooi ET, Natarajan S. A review of the Scaled Boundary Finite Element Method for two-dimensional linear elastic fracture mechanics. Eng Fract Mech 2018;187:45–73. [76] Song C, Wolf J. Consistent infinitesimal finite-element-cell method: out-of-plane motion. J Eng Mech 1995;121(5):613–9. [77] Wolf JP, Song C. Finite-element modelling of unbounded media. Wiley; 1996. [78] Wolf JP, Song C. The scaled boundary finite-element method – a primer: derivations. Comput Struct 2000;78(1–3):191–210. [79] Chidgzey SR, Deeks AJ. Determination of coefficients of crack tip asymptotic fields using the Scaled Boundary Finite Element Method. Eng Fract Mech 2005;72:2019–36. [80] Song C. Analysis of singular stress fields at multi-material corners under thermal loading. Int J Numer Meth Eng 2006;65(5):620–52. [81] Chen S, Li Q, Liu Y, Xue Z. Mode III 2-D fracture analysis by the Scaled Boundary Finite Element Method. Acta Mech Solida Sinica 2013;26(6):619–28. [82] Yang Z. Fully automatic modelling of mixed-mode crack propagation using Scaled Boundary Finite Element Method. Eng Fract Mech 2006;73(12):1711–31. [83] Yang ZJ, Deeks AJ. Fully-automatic modelling of cohesive crack growth using a finite element-scaled boundary finite element coupled method. Eng Fract Mech 2007;74:2547–73. [84] Dai S, Augarde C, Du C, Chen D. A fully automatic polygon Scaled Boundary Finite Element Method for modelling crack propagation. Eng Fract Mech 2015;133:163–78. [85] Ooi ET, Natarajan S, Song C, Tin-Loi F. Crack propagation modelling in functionally graded materials using scaled boundary polygons. Int J Fract 2015;192(1):87–105.

21

Engineering Fracture Mechanics xxx (xxxx) xxx–xxx

S. Hell and W. Becker [86] [87] [88] [89] [90] [91] [92] [93] [94] [95] [96] [97] [98] [99] [100] [101] [102] [103] [104] [105] [106] [107] [108] [109] [110] [111] [112] [113] [114] [115] [116] [117] [118] [119] [120]

Deeks AJ, Cheng L. Potential flow around obstacles using the scaled boundary finite-element method. Int J Numer Meth Fluids 2003;41(7):721–41. Lehmann L, Langer S, Clasen D. Scaled Boundary Finite Element Method for acoustics. J Comput Acoust 2006;14(4):489–506. Artel J, Becker W. Coupled and uncoupled analyses of piezoelectric free-edge effect in laminated plates. Compos Struct 2005;69(3):329–35. Birk C, Song C. A continued-fraction approach for transient diffusion in unbounded medium. Comput Meth Appl Mech Eng 2009;198(33–36):2576–90. Liu J, Lin G, Wang F, Li J. The Scaled Boundary Finite Element Method applied to electromagnetic field problems. IOP Conf Ser: Mater Sci Eng 2010;10(1). 012245. Li C, Ooi ET, Song C, Natarajan S. SBFEM for fracture analysis of piezoelectric composites under thermal load. Int J Solids Struct 2015;52:114–29. Somaratna N, Ting T. Three-dimensional stress singularities in anisotropic materials and composites. Int J Eng Sci 1986;24(7):1115–34. Ghahremani F. A numerical variational method for extracting 3D singularities. Int J Solids Struct 1991;27(11):1371–86. Ghahremani F, Shih C. Corner singularities of three-dimensional planar interface cracks. J Appl Mech 1992;59(1):61–8. Dimitrov A, Andrä H, Schnack E. Efficient computation of order and mode of corner singularities in 3D-elasticity. Int J Numer Meth Eng 2001;52(8):805–27. Dimitrov A, Andrä H, Schnack E. Singularities near three-dimensional corners in composite laminates. Int J Fract 2002;115(4):361–75. Mittelstedt C, Becker W. Semi-analytical computation of 3D stress singularities in linear elasticity. Commun Numer Meth Eng 2005;21(5):247–57. Mittelstedt C, Becker W. Efficient computation of order and mode of three-dimensional stress singularities in linear elasticity by the boundary finite element method. Int J Solids Struct 2006;43(10):2868–903. Hell S, Becker W. Hypersingularities in three-dimensional crack configurations in composite laminates. Proc Appl Math Mech PAMM 2014;14(1):157–8. Hell S, Becker W. Scaled Boundary Finite Element Analysis of three dimensional crack configurations in laminate structures. Proceedings of the international conference on computational methods (ICCM), vol. 2. 2015. p. 837–46. Hell S, Becker W. Analysis of three-dimensional crack configurations in composite laminates using the Scaled Boundary Finite Element Method. In: Proceedings of the 8th international congress of croatian society of mechanics; 2015. Dimitrov A, Schnack E. Asymptotical expansion in non-Lipschitzian domains – a numerical approach. Numer Linear Algebra Appl 2002;9(6–7):467–92. Arnoldi WE. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quart Appl Math 1951;9(1):17–29. Lee Y, Im S. On the computation of the near-tip stress intensities for three-dimensional wedges via two-state M-integral. J Mech Phys Solids 2003;51(5):825–50. Labossiere PEW, Dunn ML. Fracture initiation at three-dimensional bimaterial interface corners. J Mech Phys Solids 2001;49(3):609–34. Natarajan S, Song C. Representation of singular fields without asymptotic enrichment in the extended finite element method. Int J Numer Meth Eng 2013;96(13):813–41. Natarajan S, Song C, Belouettar S. Numerical evaluation of stress intensity factors and T-stress for interfacial cracks and cracks terminating at the interface without asymptotic enrichment. Comput Meth Appl Mech Eng 2014;279:86–112. Hell S, Becker W. The use of enriched base functions in the three-dimensional Scaled Boundary Finite Element Method. In: Papadrakakis M, Papadopoulos V, Stefanou G, Plevris V, editors. Proceedings of the 7th European congress on computational methods in applied sciences and engineering (ECCOMAS, 2016). 2016. Hell S, Becker W. Determination of singularity exponents in 3D elasticity problems using enriched base functions in the Scaled Boundary Finite Element Method. PAMM 2016;16(1):137–8. Laborde P, Pommier J, Renard Y, Salaün M. High order extended finite element method for cracked domains. Int J Numer Meth Eng 2005;64(3):354–81. Zienkiewicz O, Taylor R, Fox D. The Finite Element Method for solid and structural mechanics. 7th ed. Butterworth-Heinemann; 2014. Deeks AJ, Wolf JP. An h-hierarchical adaptive procedure for the Scaled Boundary Finite-Element Method. Int J Numer Meth Eng 2002;54:585–605. Vu TH, Deeks AJ. Use of higher-order shape functions in the Scaled Boundary Finite Element Method. Int J Numer Meth Eng 2006;65:47–70. Hell S, Becker W. Energy release rates at two perpendicularly meeting cracks by use of the Scaled Boundary Finite Element Method. Proc Struct Integr 2016;2:2471–8. Ventura G, Gracie R, Belytschko T. Fast integration and weight function blending in the extended finite element method. Int J Numer Meth Eng 2009;77(1):1–29. Richardson LF. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philos Trans Roy Soc Lond A: Math Phys Eng Sci 1911;210(459–470):307–57. Richardson LF, Gaunt JA. The deferred approach to the limit. Part I. Single lattice. Part II. Interpenetrating lattices. Philos Trans Roy Soc Lond A: Math Phys Eng Sci 1927;226(636–646):299–361. Bažant ZP, Estenssoro LF. Addendum to the paper: surface singularity and crack propagation. Int J Solids Struct 1980;16(5):479–81. Benthem J. On an inversion theorem for conical regions in elasticity theory. J Elast 1979;9(2):159–69. Gross D, Seelig T. Bruchmechanik: Mit einer Einführung in die Mikromechanik. 6th ed. Springer; 2016.

22