Crack analysis using numerical manifold method with strain smoothing technique and corrected approximation for blending elements

Engineering Analysis with Boundary Elements 113 (2020) 402–415

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Crack analysis using numerical manifold method with strain smoothing technique and corrected approximation for blending elements Feng Liu a,b, Kaiyu Zhang a, Dongdong Xu c,∗ a

State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China Department of Civil Engineering, Monash University, Clayton, VIC 3800, Australia c Key Laboratory of Geotechnical Mechanics and Engineering, Ministry of Water Resources, Yangtze River Scientific Research Institute, Wuhan 430010, China b

a r t i c l e

i n f o

Keywords: Numerical manifold method Edge-based strain smoothing Blending elements Stress intensity factor Crack propagation

a b s t r a c t Although numerical manifold element (NMM) is particularly suitable for crack analysis, the accuracy of standard NMM on traditional triangular mesh is not good enough. On the one hand, the stiffness of NMM is overly stiff, leading to inferior results, similar to the linear triangular FEM models. On the other hand, blending elements induced issues have seldom been taken into account in the framework of NMM. In this paper, these two issues are separately addressed. For the first issue, the edge-based strain smoothing technique (originated from the ES-FEM, Liu et al., 2009) is tailored for crack analysis and incorporated to improve the accuracy of NMM (abbreviated as ES-NMM). In this way, the stiffness of NMM can be much softer and much more accurate results can be obtained without increasing the DOFs. To study the second issue, different kinds of enrichment schemes are considered, such as one-layer enrichment, two-layer enrichment and corrected enrichment (enlightened by corrected XFEM, Fries, 2008). A number of numerical examples are carried out to validate the effectiveness and accuracy of the present ES-NMM. It is found that, by taking advantages of ES-FEM, corrected XFEM and NMM, the corrected ES-NMM yields the most accurate results and thus is suggested.

1. Introduction It is of practical significance to study the crack initiation, propagation and coalescence of solid structure in order to predict its life span. Due to the limitation of analytical and experimental methods for crack analysis, great efforts have been devoted to establish different numerical methods, including the finite element method (FEM) [3–5], the boundary element method (BEM) [6,7], the meshfree method [8,9], the extended finite element method (XFEM) [10,11], the numerical manifold method (NMM) [12,13], just to name a few. FEM has been widely used to solve crack problems. In FEM, the finite element mesh needs to coincide with the crack surfaces. To simulate crack propagation, a local remeshing [14] may be required, which makes the simulation complicated and time-consuming. Based on FEM, different methods have been proposed and used for crack problems, including polygonal FEM [15], virtual element method [16], generalized FEM [17], and so on.

∗

BEM has also been developed to solve crack propagation problems. A big advantage is that it could reduce the dimensions of the problem and simplify the complexity. However, it is difficult for BEM to solve nonlinear problems. By combining BEM with FEM, the scaled boundary finite element method [18,19] proposed by Song and Wolf seems to be suitable for crack problems. The meshfree method is an excellent alternative, which only utilize nodes to discretize the problem domain and thus is free from mesh related problems. This makes meshfree method more suitable for crack analysis. Lots of important work has been done to apply meshfree methods to crack simulations [8,9,20–25]. However, most meshfree methods suffer from difficulties in accurate numerical integration and high computational cost. In order to eliminate the remeshing procedure of FEM, the XFEM was proposed [10,11], in which the Heaviside functions and the asymptotic crack-tip functions are incorporated into FEM to account for the cracks. The main feature of XFEM is that the mesh does not need to conform

Corresponding author. E-mail address: [email protected] (D. Xu).

https://doi.org/10.1016/j.enganabound.2020.01.015 Received 17 October 2019; Received in revised form 1 January 2020; Accepted 29 January 2020 Available online 14 February 2020 0955-7997/© 2020 Elsevier Ltd. All rights reserved.

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

with the crack surfaces, and the same mesh can be used during crack propagation. Several improvements have been made to enhance the accuracy and convergence rate of XFEM [2,26–29]. Up to now, XFEM becomes a powerful tool and has been applied to many complex crack problems [30–34]. Other similar methods for crack problems include Phantom-node method [35], extended isogeometric analysis [36], and so on. The NMM, which can be viewed as a combination of the FEM and the discontinuous deformation analysis (DDA) [37], was first proposed by Shi [38] to solve both continuum and discontinuum problems in a unified way. By adopting dual cover systems, namely, the mathematical cover and physical cover, NMM can handle complex cracks in a nature way. NMM has been successfully applied to different fields, e.g. crack problems [39–47], seepage flow [48,49], plate problems [50], slope stability analysis [51,52], fluid-structure interaction [53], to name a few. In most studies, NMM constructs its cover system based on regular triangular mesh due to its simplicity. However, such a practice leads to poor accuracy since the strain is constant in each manifold element (for definition of manifold element, see Section 2.1). To this end, different methods have been developed to enhance the accuracy of NMM. By using high-order cover functions, NMM can be easily extended to high order [54] or Hermitian from [55] for better accuracy. However, this will significantly increase the number of degrees of freedom (DOFs). Another feasible way is to incorporate the “FE-Meshfree” technique into NMM to develop continuous nodal stress NMM [56]. Instead of finite element mesh, the influential domains of moving least squares (MLS) nodes can also be adopted to form the cover system, leading to MLS-based NMM [12,48,57]. While the MLS-based NMM inherits the advantage of NMM and meshfree method, it also has the drawbacks of meshfree method. Another aspect to be mentioned is the mesh refinement with NMM, which is quite useful when steep gradients or strong singularities involved. Several works have been done, such as multilayer covers method [58], domain decomposition method [59], local refinement [60,61] and h-adaptive analysis [62]. Recently, an edge-based smoothed NMM (ES-NMM) [63] has been proposed and applied to static, free and forced vibration analyses. It is found that, based on triangular mesh, the ES-NMM has superconvergence and ultra-accuracy without increasing DOFs. The core of ES-NMM is the strain smoothing technique, which is first proposed by Chen et al. [64] in meshfree method. Based on Chen et al.’s work [64], Liu et al. further developed the strain smoothing a lot in both FEM scheme [65–69] and meshfree scheme [70,71]. More importantly, they developed the fundamental theory, i.e. the G space theory and weakened weak formulation [70,72,73]. More details on strain smoothing technique and its applications in FEM scheme can be found in Liu and Nguyen’s book [74] and reviews [75,76]. In this study, to make full use of the high accuracy of ES-NMM and the advantage of NMM in dealing with discontinuity, the ES-NMM is further applied to linear fracture problems. Problems arising from blending element, seldom touched in NMM, are also discussed and examined. The proposed corrected ESNMM, which is free from problems in blending elements, can greatly enhance the accuracy of NMM for crack problems. Compared with the extended finite element method with edge-based strain smoothing [27], the smoothing domains in corrected ES-NMM do not need to be split into triangular sub-domains before applying the strain smoothing technique. Furthermore, inheriting the unique cover systems from NMM, corrected ES-NMM is more suitable for complex cracks. Compared with phantomnode method with edge-based strain smoothing [77], in which the crack tips must fall on the edges, the crack tips can locate everywhere in the proposed method. The rest of the paper is organized as follows. The fundamental of NMM with edge-based strain smoothing technique for crack analysis is introduced in Section 2. To study the blending elements related problems, different enrichment schemes are present in Section 3. Other numerical implementations are presented in Section 4. Six numerical ex-

amples will be carried out in Section 5 to validate the effectiveness and robustness of the proposed method. Some conclusions and discussions are made in Section 6.

2. Fundamental of ES-NMM for crack problems The fundamental of edge-based smoothed NMM (ES-NMM) will be discussed here. Note that, since cracks are involved, the governing equations and related implementation details in this paper are quite different from [63]. First, singular physical patches are introduced here to describe the singularity of crack tips. Second, the strain smoothing technique is only applied on some of the manifold elements. To this end, ES-NMM for crack analysis is briefly introduced in this section. 2.1. Cover systems of NMM NMM introduces dual cover systems, called the mathematical cover (MC) and the physical cover (PC), respectively. The dual cover systems are the unique feature of NMM that distinguishes it from FEM. Mathematical cover consists of a number of simply-connected domains, each of which is called a mathematical patch (MP). Note that MPs are independent of the geometrical details of the problem domain Ω, and can overlap each other. The only requirement is that all of them together must cover the entire domain Ω. The physical patches (PPs) are generated by the intersection of each MP with the components of the problem domain, including the boundaries, holes, material interfaces, cracks, and so on. In this process, one MP can be divided into more than one PPs. However, if a MP locates totally outside Ω, no PP will be generated. All the PPs constitute the PC, which will exactly cover the whole problem domain Ω. A manifold element (ME) can be formed by the intersection of several nearby PPs. The MEs are used for integrating the weak form of the problem in traditional NMM. This, however, is not the case when the edge-based strain smoothing technique is incorporated. Details will be introduced in the following section. To illustrate the above concepts more clearly, we take a domain with a crack shown in Fig. 1 as an example. The problem domain Ω consists of the physical boundaries (denoted by blue thick lines) and a crack (denoted by red thick lines). The regular triangular mesh, which is called the mathematical mesh here, is used to generate the mathematical cover. Each node in the mathematical mesh is called a star. A MP consists of all elements linked at a given star in the mathematical mesh. In most cases, a MP is a regular hexagon containing 6 triangles. Several typical stars (marked as small solid circles) along with their MPs (marked as colored regular hexagons) are also plotted in Fig. 1. Since MP1 , MP2 and MP3 locate in Ω, only one PP can be obtained from each of them. These three PPs are denoted by PP1-1 , PP2-1 and PP3-1 , which are exactly the same as MP1 , MP2 and MP3 , respectively. The overlap of PP1-1 , PP2-1 and PP3-1 forms the manifold element E1 . MP4 is cut by the crack completely, thus two PPs, PP4-1 and PP4-2 are obtained. Similarly, PP6-1 and PP6-2 are obtained from MP6 . MP5 is also cut by the crack, however, it hasn’t been divided into two pieces, thus only MP5-1 is obtained. The overlap of PP4-1 , PP5-1 and PP6-1 forms E2 , while the overlap of PP4-2 , PP5-1 and PP6-2 forms E3 . By the intersection of MP7 with Ω, only one physical patch PP7-1 inside Ω is generated. Similarly, PP8-1 and PP9-1 are generated from MP8 and MP9 , respectively. The overlap of PP7-1 , PP8-1 and PP9-1 forms manifold element E4 . There are two types of PPs, the nonsingular patches and the singular patches (also termed as enriched patches). Most PPs far from crack tips are nonsingular patches, say PP1-1 and PP 8 –1 . PPs near the crack tips are singular patches. This is different from the practice in NMM literature [40,42,44], where only the patches containing one or more crack tips are defined as singular patches. In this paper, the singular patches depend on the enrichment schemes introduced in Section 3. For example, PP5-1 is always a singular patch in all enrichment schemes. However, for PP6-1 and PP 6 –2 , sometimes they are singular patches, sometimes not. 403

F. Liu, K. Zhang and D. Xu

MP1(PP1-1)

Engineering Analysis with Boundary Elements 113 (2020) 402–415

MP2(PP2-1)

DOFs. The accuracy of present method will be improved by the edgebased strain smoothing technique discussed later. For a singular patch PPm , the local approximation function is enriched with the asymptotic crack-tip functions, with the following form

MP3(PP3-1)

Physical boundary MP4

E1

PP4-1

𝒖𝑚 (𝐱) = 𝑻 (𝐱)𝒅 𝑚 +𝑭 e (𝐱)𝒈𝑚 ,𝐱 ∈ 𝑃 𝑃m

where gm is the enriched DOFs with eight components, and Fe (x) is the matrix of singular bases, [ ] Ψ1 0 Ψ2 0 Ψ3 0 Ψ4 0 𝑭 e (𝐱 ) = (7) 0 Ψ1 0 Ψ2 0 Ψ3 0 Ψ4

PP4-2 MP5

PP5-1

[ Ψ1 E2 E3

MP6

Crack

MP8

PP8-1

MP7

] [√ 𝜃 Ψ4 = 𝑟 cos 2

Ψ3

√

𝑟 sin

𝜃 2

√

𝑟 cos

3𝜃 2

√

𝑟 sin

3𝜃 2

] (8)

𝑛

PP6-2

PP9-1

Ψ2

in which (r, 𝜃) are the polar coordinates in the local coordinates system with the origin at the crack tip. The final global approximation over a given point x can be expressed as the product of local approximation function and the partition of unity,

PP6-1

E4

MP9

(6)

𝒖ℎ (𝐱) =

𝑝 ∑

𝑚

PP7-1

𝝓𝑚 (𝐱)𝒖𝑚 (𝐱), for 𝐱 ∈ Ω

(9)

here np is the total number of physical patches that cover x. Considering that triangular mesh is used and suppose that three physical patches associated with x are PPi , PPj and PPk , the global approximation can also be rewritten in matrix form E1 =PP1-1 ∩ PP2-1 ∩ PP3-1

E2 =PP4-1 ∩ PP5-1 ∩ PP6-1

𝒖ℎ (𝐱)=𝑵 𝒑

(10)

E4 =PP7-1 ∩ PP8-1 ∩ PP9-1

E3 =PP4-2 ∩ PP5-1 ∩ PP6-2

where p is the vector of generalized DOFs [ ] 𝒑𝑇 = 𝒑𝑇𝑖 𝒑𝑇𝑗 𝒑𝑇𝑘

(11)

and N is the matrix of shape function [ ] 𝑵 = 𝑵𝑖 𝑵𝑗 𝑵𝑘

(12)

Fig. 1. Illustration of the problem domain and the cover systems. (Modified from [63]). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.2. The partition of unity and approximation

For nonsingular patch, 𝒑𝑚 = 𝒅 𝑚 ,

Over each mathematical patch MPi , a partition of unity (also called weight function or shape function) 𝜙i (x) is defined 𝜙𝑖 (𝐱) = 0, 𝑖𝑓 𝐱 ∉ 𝑀 𝑃𝑖 ;

(1)

0 ≤ 𝜙𝑖 (𝐱) ≤ 1, 𝑖𝑓 𝐱 ∈ 𝑀 𝑃𝑖 ;

(2)

𝑛 ∑ 𝑖=1

𝜙𝑖 (𝐱) = 1, 𝑖𝑓 𝐱 ∈ Ω

𝜺 = 𝐿𝑑 𝒖

1 0

0 1

𝑥 0 0 𝑥

𝑦 0 0 𝑦

[ 𝐿𝑑 ≡

𝜕 𝜕𝑥

0

0 𝜕 𝜕𝑦

𝜕 ]𝑇 𝜕𝑦 𝜕 𝜕𝑥

(16)

Substituting Eq. (10) into Eq. (15) leads to 𝜺 = 𝑩𝒑

(4)

] ... ...

(15)

with

(17)

with compatible B matrix

here dm is the DOFs of physical patch PPm , and T(x) is the matrix of basis functions. If polynomial basis is used, T(x) has the following form [

(13)

here pm is a vector of 10 × 1, and Nm is a matrix of 2 × 10. The compatible strain 𝜺 can then be computed from the displacement u,

Each physical patch PPi-j inherits the partition of unity function from the corresponding mathematical patch MPi , denoted as 𝜙i − j (x). For the sake of convenience, all the physical patch PPi-j and corresponding partition of unity function are renumbered with single indices and denoted as PPm and 𝜙m (x), respectively. Over each nonsingular patch PPm , a local approximation function (also called cover function), denoted by um (x), can be expressed as

𝑻 (𝐱) =

𝑚 = 𝑖, 𝑗, 𝑘

where dm is the displacement DOFs, and I2 is the 2 × 2 identity matrix. For singular patch, [ ] [ ] 𝑇 𝒑𝑚 = 𝒅 𝑇𝑚 𝒈𝑇𝑚 , 𝑵 𝑚 = 𝜙𝑚 𝑰 2 𝜙𝑚 𝑭 e (𝐱) , 𝑚 = 𝑖, 𝑗, 𝑘 (14)

(3)

𝒖𝑚 (𝐱) = 𝑻 (𝐱)𝒅 𝑚 ,𝐱 ∈ 𝑃 𝑃m

𝑵 𝑚 = 𝜙𝑚 𝑰 2 ,

𝑩 =𝐿 𝑑 𝑵

(18)

2.3. The edge-based smoothing domains and strain smoothing technique (5) Among the smoothed finite element method, ES-FEM is the most outstanding one, which is much more accurate than the linear FEM using the same triangular mesh and even more accurate than the FEM using quadrilateral elements. In order to simulate crack problems, singular ES-FEM [78] and extended ES-FEM [27] have been developed.

It is convenient to construct higher order local approximation function in NMM just by increasing the order of T(x). However, higher order T(x) leads to much more DOFs. In this paper, we take T(x) as the second order identity matrix, which means each nonsingular patch has only two 404

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

as follows,

In this section, the strain smoothing technique will be introduced in the framework of NMM. The construction of edge-based smoothing domains in NMM is quite similar to that of ES-FEM. The only difference is that the manifold elements instead of FEM elements are regarded as the basic units. Two steps are involved. First, we divide each manifold element into a couple of sub-triangles by linking its nodes to the centroid of the element. The centroid can be computed as the average coordinates of the manifold element nodes. Then, the smoothing domain corresponding to a certain edge is formed by the sub-triangles sharing the same edge. In this process, the manifold element can be in any shape with arbitrary numbers of edge. Take a square domain with a crack as an example, the typical edgebased smoothing domains (denoted as red dash area) for both interior edge (say edge l2 ) and boundary edge (say edge l1 ) are illustrated in Fig. 2. The empty blue circles represent the centroids of the manifold elements. It is obvious that for boundary edges, the smoothing domain consists of one sub-triangle, while for interior edge, the smoothing domain consists of two sub-triangles. The smoothing domain of edge l3 will be discussed in Section 2.4. After constructing the smoothing domains, the smoothed strain field now can be evaluated using the displacements. First a smoothing operator is performed to the gradient of displacement on the smoothing domain Ω(k) ∫Ω(𝑘)

∇𝑢(𝐱)Φ(𝑘) (𝐱)𝑑 Ω

𝑢(𝐱)𝒏(𝐱)Φ(𝐱)𝑑 Γ

(23)

⎡𝑏̄ 𝑖𝑥 0 ⎤ ̄ 𝑖 = ⎢0 𝐁 𝑏̄ 𝑖𝑦 ⎥ ⎢̄ ⎥ ⎣𝑏𝑖𝑦 𝑏̄ 𝑖𝑥 ⎦

(24)

1 𝑏̄ 𝑖ℎ = 𝜙𝑖 𝑛ℎ 𝑑 Γ(𝑘) , ℎ = 𝑥, 𝑦 𝐴(𝑘) ∫Γ(𝑘)

(25)

In NMM, the MC does not need to conform to the physical details. On the one hand, this feature allows the mesh independent of boundaries, cracks, etc., simplifying the mesh generation to a great extent. On the other hand, it gives rise to the difficulty of imposing displacement boundary conditions. In this study, by using the penalty function to fulfill the displacement boundary conditions, the potential energy can be written as 𝜋(𝒖) =

1 𝑇 𝜺 𝐃𝛆 𝑑Ω − 𝒖𝑇 𝐛 𝑑Ω − 𝒖𝑇 𝒕̄𝑑𝑆 ∫Ω 2 ∫Ω ∫Γ𝑡 +

1 𝑘 (𝒖 − 𝒖̄ )𝑇 (𝒖 − 𝑢̄ )𝑑𝑆 ∫Γ𝑢 2 𝑝

(26)

where 𝜺 is the strain vector, D is the elasticity matrix, b is the body force, kp is the penalty factor, Γt represents the stress boundary, on which the traction vector 𝒕̄ is loaded, and Γu is the displacement boundary on which the displacement is given by 𝒖̄ . Before deriving the final discretization, attentions should be paid to the effectiveness of the strain smoothing technique. As pointed out by Bordas et al. [79] that, strain smoothing in enriched approximation is only beneficial when the enrichment functions are polynomial. For nonpolynomial enrichment of near-tip functions, strain smoothing leads to inferior results compared to the standard enriched method. This can be understood that the strain (and also the stress) field near the crack tips changes rapidly, which does not conform to the strain smoothing technique, where the strains are assumed to be constant in each smoothing domain. A solution for this issue was proposed in [27], where the smoothing domains are further split into triangular sub-domains before applying the strain smoothing technique, and more integration points are used to perform the boundary integration. This, however, make the implement quite complicated and tedious. In this paper, the strain smoothing technique is only applied on manifold elements which are not influenced by the crack tip enrichment (nonenriched elements). For enriched elements, e.g. the gray areas as shown in Fig. 2, the compatible strain in Eq. (15) is used. Therefore, the smoothing domain of interior edge on the outer boundary of enriched elements consists of only one sub-triangle, similar to the boundary edge. This is illustrated by the smoothing domain of edge l3 in Fig. 2. This treatment does not cause problems, since the displacement is compatible on all the elements. The number of enriched elements depends on the enrichment scheme discussed later. The union of all the enriched elements is

(19)

here A(k) is the area of Ω(k) . Then, integration by parts of Eq. (19) with the smoothing function defined in Eq. (20) for the right-hand side leads to ∫Γ(𝑘)

̄ is the smoothed strain matrix, where 𝐁 [ ] ̄= 𝐁 ̄ 𝑖, ⋯ 𝐁 ̄𝑚 ̄ 1, ⋯ 𝐁 𝐁

2.4. Variational formulation and its discretization

where Φ(k) (x) is the smoothing function and can be simply defined as [1] { 1∕𝐴(𝑘) , 𝐱 ∈ Ω(𝑘) Φ(𝑘) (𝐱)= (20) 0, 𝐱 ∉ Ω( 𝑘 )

∇𝑢ℎ (𝐱) =

(22)

here m is the number of physical patches associated with the smoothing domain Ω(k) , and 𝜙i is the shape function value of physical patch PPi . Normally, m = 4 for interior edge and m = 3 for boundary edge. This is still illustrated in Fig. 2. Five physical patches can be obtained from the six stars numbered from 1 to 5. It is obvious that the smoothing domain of boundary edge l1 is covered by physical patch 2, 4 and 5, while the smoothing domain of interior edge l2 is covered by physical patch 1, 2, 5 and 3. Instead of using Eq. (22), the smoothed strain can also be obtained by the area weighted averaging of compatible strain of the two elements related to the edge.

Fig. 2. Illustration of edge-based smoothing domains in NMM. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

∇𝑢ℎ (𝐱) =

̄𝒑 𝜺̄ (𝑘) = 𝐁

(21)

where Γ(k) is the boundary of Ω(k) with outward unit normal n(x). This step changes the area integration in Eq. (19) into line integration in Eq. (21) and eliminates the calculation of the gradient of displacement. Without considering the singular patches (as will be seen later, the strain smoothing technique is not applied to singular patches), from Eq. (21), the smoothed strain on smoothing domain Ω(k) can be obtained 405

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

blending elements

blending elements

reproducing elements

enriched PPs

reproducing elements

enriched PPs represent that two PPs are obtained from this star

Fig. 3. The first enrichment scheme: one-layer enrichment.

Fig. 4. The second enrichment scheme: two-layer enrichment.

denoted as Ωen , while the non-enriched elements constitute Ωnon . In this way, by substituting Eqs. (10), (17) and (22) into Eq. (26) and putting 𝛿𝜋(u)= 0, the ES-NMM system of equations can be obtained, 𝑲𝒉 = 𝒒

(27)

where h is the global vector of DOFs, including displacement terms and enriched terms, K is the global stiffness matrix, 𝑲=

∫Ω𝑒𝑛

𝑩 𝑇 𝑫 𝐵 𝑑Ω +

∫Ω𝑛𝑜𝑛

𝑇 𝑩̄ 𝑫 𝑩̄ 𝑑Ω + 𝑘𝑝

∫Γ𝑢

𝑵 𝑇 𝑵 𝑑Ω

(28)

and q is the global generalized force vector, 𝒒=

∫Ω

𝑵 𝑇 𝒃𝑑Ω +

∫Γ𝑡

𝑵 𝑇 𝒕̄𝑑𝑆 + 𝑘𝑝

∫Γ𝑢

𝑵 𝑇 𝒖̄ 𝑑𝑆

(29)

Note that, in above equations, the strain smoothing technique is only used for computing the stiffness matrix of Ωnon . 3. Different enrichment schemes Blending elements are the elements that blend the enriched areas with the rest of the domain [2]. In the blending elements, the partition of unity is not satisfied, and thus unwanted terms are introduced in the approximation, degrading the accuracy of XFEM significantly. The problem arising from blending elements has been widely discussed in XFEM and several modifications have been made to improve the performance [2,28,29]. In NMM, blending elements related problems for material discontinuities has been discussed by An et al. [80]. However, to the authors’ knowledge, blending elements involving crack-tip enrichment have not been studied in the framework of NMM. In this paper, three different enrichment schemes are considered to verify the blending element problems. The first enrichment scheme is called one-layer enrichment, which is commonly used in NMM [40,44]. In this scheme, the PPs that cover the crack tips are enriched, as shown in Fig. 3. The black solid circles represent the enriched PPs. Normally, three PPs are enriched when triangular mesh is adopted. Following the description in [2], the reproducing elements are those elements on which the enrichment functions in Eq. (8) can be reproduced exactly. In blending elements, the situation seems the same as XFEM. It will be verified through numerical examples.

blending elements

reproducing elements

type 1 PPs (the enriched PPs) type 2 PPs Fig. 5. The third enrichment scheme: corrected enrichment.

The second enrichment scheme is called two-layer enrichment, in which two-layer PPs are enriched, as shown in Fig. 4. It is expected that the second enrichment scheme surpasses the first one since more enriched DOFs are involved. It should be noted that there are 14 rather than 12 PPs are enriched. This is because two PPs are obtained from each of the star denoted by blue circle. As previously mentioned, some PPs that do not cover the crack tips are also enriched. The third enrichment scheme is called the corrected enrichment, which is first proposed by Fries [2] to deal with the blending elements. In this method, two special types of PPs are defined (see Fig. 5). 406

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

The type 1 PPs are the enriched PPs similar to those of the second enrichment scheme. A ramp function [2] is defined based on the type 2 PPs (denoted as I□ ) 𝑅 (𝐱 ) =

∑ 𝑖∈𝐼 □

𝜙□ 𝑖 (𝐱 )

(30)

where 𝜙□ 𝑖 (𝐱) is the shape function value of type 2 PPs. It can be seen that, R(x) equals to 1 in the reproducing elements, varies from 1 to 0 in the blending elements, and equals to 0 otherwise. The matrix of singular bases in Eq. (6) is then modified as e 𝑭̃ (𝐱) = 𝑅(𝐱)𝑭 e (𝐱)

A

A3 A2 A1

(31)

As a result, in the blending elements, the modified enrichment function varies continuously between standard enrichment functions to 0. The number of enriched DOFs of this scheme only depends on type 1 PPs, and equals to that of the two-layer enrichment, thus, the second scheme can be used to examine the effectiveness of the corrected enrichment. According to previous definitions, Ωen is made up of the blending elements and the reproducing elements. Comparing Fig. 3 to Fig. 5, we can conclude that the second scheme has the largest Ωen , and that the third enrichment scheme increases the DOFs without increasing the area of Ωen . It should be noted that a correct enrichment requires the crack nearby to be straight. Hence, a small Ωen will be better for kinked cracks or crack propagation problems. In this aspect, the third scheme overcomes the second one. By the way, the problem of non-straight cracks can be solved by special treatments, such as the so-called crack mapping method [10] or the new sign convention method [46]. However, these special treatments are not employed in this paper.

Fig. 6. Illustration of the domain form of interaction integrals.

with radius R, as shown in Fig. 6. The radius R is taken as 𝑅 = 𝑅𝑑 ℎ

(32)

where h is the edge length of the mathematical mesh, Rd is a dimensionless factor which determines the size of the circle. The interaction integral should be path independent in theory. However, it is path dependent due to the numerical error. So the selection of Rd will be studied. The second factor is the computation method of actual stresses. The stresses can be computed with or without strain smoothing technique. If the strain smoothing technique is applied, each element contains several sub-parts of the smoothing domains, e.g., for element A, three subparts A1 , A2 and A3 are included, also shown in Fig. 6. Since the strain smoothing technique is not employed for enriched elements, the integral domain must not contain any enriched elements. This means a relatively larger Rd should be chosen. Conversely, if the strain smoothing technique is not used, there is no limitation of Rd .

4. Numerical implementation 4.1. Numerical integration There are two kinds of numerical integrations involved, the line integration defined in Eq. (25) and the conventional domain integration. For the line integration, since the strain smoothing technique only applies to non-enriched elements, one Gaussian point (the midpoint) for each boundary segment is sufficient to integrate Eq. (25) [63]. For domain integration, where the crack tip enrichment functions are involved, high order integration is required. First, consider the case when the elements don’t contain any crack tips. If an element is not a triangle, it is first partitioned into sub-triangles. Then 13 Gauss points are used on each sub-triangle. For tip elements, the integration scheme proposed by Zheng and Xu [46] is adopted, which can treat the 1/r singularity of integrand. As a contrast, in ESm-XFEM [27], the smoothing domains in the vicinity of crack tips are first spilt into small smoothing cells, and then high order Gauss point is used on each boundary segment of the smoothing cells. This makes the implement unnatural and complicated.

4.3. Crack propagation criterion To simulate the process of crack propagation, a crack propagation criterion needs to be specified to determine the crack growth direction. The maximum hoop-stress criterion [81] is adopted in this paper. This criterion advocates that the crack will propagate in the direction perpendicular to the maximum circumferential stress. Cracks will propagate with a given length. 5. Numerical examples In this section, several numerical examples are carried out to assess the efficiency and accuracy of the proposed ES-NMMs. Since regular meshes can be used without considering both the problem boundaries and cracks, regular triangular mesh will always be adopted in all the examples. Always using regular mesh is one of the features of NMM. The element size is defined as the side length of the triangular element. The penalty factor has some influence on the results. If the penalty factor is too small, the constraints cannot be imposed well. However, an overlarge penalty factor can lead to ill-conditioned stiffness matrix. In all the examples, the penalty factor is set to 103 times the Young’s modulus E, which is consistent with MLS-based NMM [82]. To evaluate the numerical results, the following four error indicators are used, which are the relative error in displacement norm ed and in energy norm ee , normalized SIFs MI and MII , respectively: √ √ √ ∫ (𝒖𝑟𝑒𝑓 − 𝒖𝑛𝑢𝑚 )2 𝑑Ω 𝑒𝑑 = √ Ω (33) ∫Ω (𝒖𝑟𝑒𝑓 )2 𝑑Ω

4.2. SIFs evaluation Stress intensity factor (SIF) plays a big role in crack analysis. In this paper, the SIFs are calculated by the domain form of interaction integrals [11]. In this method, two states of a cracked body are considered. One is the actual state based on the numerical results, while the other one represents an auxiliary state from the asymptotic fields for modes I and II. The value of SIFs can be obtained based on the interaction integral of these two states. For details see [11]. Herein, the influences of two factors on the SIFs are considered. The first one is the size of integral domain, which is determined by a circle 407

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

σ

Crack a

L

a

L

σ

(a)

(b)

Fig. 7. (a) Infinite plate under uniform tension and (b) typical mesh with element size h = 0.8 m.

√ √1 √ ∫ (𝜺𝑟𝑒𝑓 − 𝜺𝑛𝑢𝑚 )𝑇 𝑫 (𝜺𝑟𝑒𝑓 − 𝜺𝑛𝑢𝑚 )𝑑Ω √ Ω 𝑒𝑒 = √ 2 1 ∫ (𝜺𝑟𝑒𝑓 )𝑇 𝑫 (𝜺𝑛𝑢𝑚 )𝑑Ω 2 Ω 𝑀I =

𝑀II =

𝐾I𝑛𝑢𝑚 𝐾I𝑟𝑒𝑓 𝐾II𝑛𝑢𝑚 𝐾II𝑟𝑒𝑓

𝐾 𝜃 𝜃 3𝜃 𝜎𝑥𝑦 (𝑟, 𝜃) = √𝐼 sin cos cos . 2 2 2 𝑟

(34)

The analytical near-tip displacement fields are ( ) 2(1 + 𝑣) 𝐾𝐼 √ 𝜃 𝜃 𝑟 cos 2 − 2𝑣 − cos2 , √ 2 2 2𝜋 𝐸 ( ) √ 2(1 + 𝑣) 𝐾𝐼 𝜃 𝜃 𝑢𝑦 (𝑟, 𝜃) = √ 𝑟 sin 2 − 2𝑣 − cos2 . 𝐸 2 2 2𝜋

(35)

𝑢𝑥 (𝑟, 𝜃) =

(36)

(38a) (38b)

√ In Eqs. (37) and (38), the stress intensity factor 𝐾𝐼 = 𝜋𝑎, E and v are the Young’s modulus and the Poisson ratio, respectively. A square with edge length L containing a half crack of length a is taken as the computation zone. The analytical near-tip displacements in Eq. (38) are prescribed along the boundary of the square. In the simulation, the following parameters are taken, L = 10 m, a = 5 m, E = 107 N/m2 , v = 0.3, 𝜎=104 Pa, and plane strain conditions are assumed. To compute the convergences rate, different element sizes, such as h = 1.2, 0.8, 0.6, 0.4,0.3, 0.2 m are used. A typical mesh with h = 0.8 m is shown in Fig. 7(b). The convergences rates of different methods in displacement norm are shown in Fig. 8. The convergence rate of standard NMM is 0.91, which is worse than all the ES-NMMs. Meanwhile, the error level of standard NMM is greater than that of the ES-NMMs. This shows that the strain smoothing technique can improve the performance of NMM. Due to the increase of enriched DOFs, two-layer ES-NMM obviously outcomes one-layer ES-NMM. It is interesting to compare the results of two-layer ES-NMM with corrected ES-NMM since they have the same number of DOFs. The corrected ES-NMM performs better than two-layer ES-NMM in both terms of convergence rate and error level, although the enriched area of corrected ES-NMM is smaller than that of two-layer ES-NMM. This confirms the effectiveness of the corrected enrichment. Similar conclusions can be drawn from the rate of convergences in energy norm, as shown in Fig. 9. Although the convergence rates of all the methods are about 0.5, which conform to the results of XFEM [2], ES-NMMs produce better accuracy. In order to achieve an optimal convergence rate, the geometrical enrichment [83] may be

here the superscript ‘‘ref’’ denotes the reference solution or analytical solution, and the superscript ‘‘num’’ indicates a numerical solution. For the purpose of comparison, we will employ the following methods or schemes: (1) standard NMM: strain smoothing technique is not employed, with one-layer enrichment (see Fig. 3). This is the one frequently adopted in NMM references (e.g. [40,44]); (2) one-layer ES-NMM: strain smoothing technique is employed, with one-layer enrichment (see Fig. 3); (3) two-layer ES-NMM: strain smoothing technique is employed, with two-layer enrichment (see Fig. 4); (4) corrected ES-NMM: strain smoothing technique is employed, with corrected enrichment (see Fig. 5). 5.1. Infinite plate under uniform tension Consider an infinite plate containing a straight crack of length 2a and subjected to a remote tension stress 𝜎, as shown in Fig. 7(a). The analytical stress fields in terms of polar coordinates (r, 𝜃) centered at the crack tip are as follows, ( ) 𝐾 𝜃 𝜃 3𝜃 𝜎𝑥𝑥 (𝑟, 𝜃) = √𝐼 cos 1 − sin sin , 2 2 2 𝑟 ( ) 𝐾 𝜃 𝜃 3𝜃 𝜎𝑦𝑦 (𝑟, 𝜃) = √𝐼 cos 1 + sin sin , 2 2 2 𝑟

(37c)

(37a) (37b) 408

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

Fig. 8. Convergence rate in displacement norm for infinite plate under uniform tension.

Fig. 10. Comparison of the efficiency in terms of displacement error norm for infinite plate under uniform tension.

Fig. 9. Convergence rate in energy norm for infinite plate under uniform tension.

needed, which however leads to ill-conditioned systems and much more DOFs [29]. We also record the CPU times to compare the computational efficiency of different methods. The comparison of the efficiency in terms of displacement error norm and energy error norm are shown in Figs. 10 and 11, respectively. The efficiency of one-layer ES-NMM is equal or a bit better than standard NMM. This demonstrates that the increase of computational cost induced by edge-based smoothing operations is acceptable. The computational efficiencies of two-layer ES-NMM and corrected ES-NMM are obviously better than standard NMM and one-layer ES-NMM. As expected, corrected ES-NMM is the most efficient one in both displacement and energy error norm, which again prove the effectiveness of the corrected enrichment. Furthermore, the influence factors on SIFs are investigated. To consider the influence of R defined in Eq. (32), three different values of Rd , such as Rd =2, 2.5, and 3 are taken. The results of Normalized KI values using interaction integrals without strain smoothing technique are tabulated in Table 1. The corrected ES-NMM performs best, then two-layer ES-NMM, followed one-layer ES-NMM. Although the standard NMM is

Fig. 11. Comparison of the efficiency in terms of energy error norm for infinite plate under uniform tension.

not sensitive to Rd , it is the least accurate one. For two-layer ES-NMM, the optimal Rd is 3, and Rd =2, 2.5 leads to larger errors. For one-layer ES-NMM and corrected ES-NMM, Rd =2.5 or 3 leads to nearly the same accuracy. When the strain smoothing technique is used in the interaction integrals, the results are tabulated in Table 2. In standard NMM, no smoothing domain exists, thus the SIFs cannot be obtained in this way. This also occurs for two-layer ES-NMM with Rd =2 or 2.5, where some elements involved in the interaction integrals belong to Ωen without smoothing domain. From the valid results in Table 2, it is found that corrected ES-NMM is slightly better than two-layer ES-NMM, and they both better than one-layer ES-NMM. By comparing Table 1 with Table 2, it seems that using or not using strain smoothing technique has little influence on the SIFs. 409

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

Table 1 Normalized KI values (MI ) for infinite plate under uniform tension using interaction integrals without strain smoothing technique. Different methods

Rd

Mesh h = 1.2 (118 PPs)

h = 0.8 (244 PPs)

h = 0.6 (406 PPs)

h = 0.4 (860 PPs)

h = 0.3 (1434 PPs)

h = 0.2 (3169 PPs)

Standard NMM

2.0 2.5 3.0

1.0229 1.0237 1.0237

1.0174 1.0174 1.0170

1.0143 1.0105 1.0120

1.0093 1.0098 1.0094

1.0059 1.0071 1.0064

1.0053 1.0033 1.0040

One-layer ES-NMM

2.0 2.5 3.0

1.0278 1.0137 1.0183

1.0255 1.0092 1.0088

1.0210 1.0081 1.0089

1.0148 1.0053 1.0071

1.0170 1.0012 1.0036

1.0105 0.9976 1.0032

Two-layer ES-NMM

2.0 2.5 3.0

1.0182 1.0187 1.0128

1.0084 1.0181 1.0128

1.0148 1.0175 1.0088

1.0132 1.0158 1.0081

1.0065 1.0143 1.0078

1.0105 1.0135 1.0037

Corrected ES-NMM

2.0 2.5 3.0

1.0173 1.0040 1.0078

1.0178 1.0005 1.0003

1.0163 1.0050 1.0026

1.0107 1.0018 1.0024

1.0148 0.9985 1.0006

1.0097 0.9974 1.0017

Table 2 Normalized KI values (MI ) for infinite plate under uniform tension using interaction integrals with strain smoothing technique. Different methods

Rd

Mesh h = 1.2 m (118 PPs)

h = 0.8 m (244 PPs)

h = 0.6 m (406 PPs)

h = 0.4 m (860 PPs)

h = 0.3 m (1434 PPs)

h = 0.2 m (3169 PPs)

Standard NMM

2.0 2.5 3.0

– – –

– – –

– – –

– – –

– – –

– – –

One-layer ES-NMM

2.0 2.5 3.0

1.0168 1.0166 1.0181

1.0118 1.0127 1.0124

1.0083 1.0079 1.0088

1.0047 1.0067 1.0065

1.0033 1.0048 1.0043

1.0016 1.0019 1.0019

Two-layer ES-NMM

2.0 2.5 3.0

– – 1.0085

– – 1.0042

– – 1.0036

– – 1.0021

– – 1.0014

– – 1.0009

Corected ES-NMM

2.0 2.5 3.0

1.0059 1.0060 1.0074

1.0026 1.0041 1.0037

1.0088 0.9995 1.0022

1.0004 1.0023 1.0021

0.9999 1.0017 1.0013

1.0002 1.0005 1.0004

“-” represents that the smoothing domains do not exist for some or all the elements involved in the interaction integrals, thus the SIFs cannot be computed.

5.2. An edge crack under shear load

τ

In the second example, we consider a rectangular plate with an edge crack subjected to a uniform shear force 𝜏 on the top, as shown in Fig. 12(a). The bottom of the plate is fixed and plane strain condition is assumed. The parameters are set as W = 7 mm, H = 8 mm, crack length a = 3.5 mm, and 𝜏=1 Pa. The reference SIFs for this problem are [27] 𝐾I𝑟𝑒𝑓

√ = 34.0 𝑃 𝑎 𝑚𝑚,

𝐾II𝑟𝑒𝑓

√ = 4.55 𝑃 𝑎 𝑚𝑚

H (39)

a Three regular triangular meshes with element size h = 1, 0.5, 0.25 mm are used. The typical mesh with h = 0.5 mm is plotted in Fig. 12(b). Again, the two influence factors on SIFs are investigated. The results of Normalized SIFs using interaction integrals without strain smoothing technique are listed in Table 3. As h decreases, the Normalized SIFs approach to 1. In terms of accuracy, corrected ES-NMM is the most accurate method, followed by two-layer ES-NMM, then one-layer ES-NMM, at last standard NMM. For corrected ES-NMM, the influence of Rd on SIFs is not obvious in this case. The results using interaction integrals with strain smoothing technique are listed in Table 4. It can be concluded that the accuracies of these two methods are comparable. Since the interaction integral without strain smoothing technique has fewer limitations, it is used to compute SIFs in the following examples. According to the results of first two examples, Rd = 2.5 is used for standard NMM, one-layer ES-NMM and corrected ES-NMM, and Rd = 3 for two-layer ES-NMM in following examples.

H W

(a)

(b)

Fig. 12. (a) Edge crack under shear force and (b) Typical mesh with element size h = 0.5 mm.

410

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

σ Table 3 Normalized SIFs for edge crack under shear force using interaction integrals without strain smoothing technique. Different methods

Normalized SIFs

Rd

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.8548 0.8556 0.8548

0.9506 0.9489 0.9493

0.9766 0.9767 0.9770

One-layer ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.9181 0.9021 0.9009

0.9807 0.9683 0.9739

0.9982 0.9823 0.9888

Two-layer ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.9442 0.9551 0.9500

0.9968 0.9995 0.9897

1.0054 1.0089 0.9978

Corrected ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.9827 0.9666 0.9639

1.0009 0.9889 0.9932

1.0082 0.9920 0.9982

Standard NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

0.9670 0.9574 0.9640

0.9643 0.9641 0.9684

1.0056 1.0108 1.0040

One-layer ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

0.9746 0.9731 0.9859

0.9681 0.9706 0.9923

1.0517 1.0100 1.0034

Two-layer ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

1.0522 0.9820 0.9456

1.0123 0.9733 0.9768

0.9794 1.0393 1.0219

Fig. 13. (a) A plate with two cracks emanating from a hole and (b) Typical mesh with a = 0.7b and 𝜃=45°.

Corrected ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

0.9738 0.9882 1.0077

0.9740 0.9810 1.0025

1.0596 1.0096 1.0028

Table 5 Normalized SIFs at crack tip A for a plate with two cracks emanating from a hole.

Standard NMM

Mesh h = 1 mm h = 0.5 mm h = 0.25 mm (187 PPs) (637 PPs) (2307 PPs)

h A

2a r B b

σ (a)

Table 4 Normalized SIFs for edge crack under shear force using interaction integrals with strain smoothing technique. Normalized SIFs

Rd

Standard NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

– – –

– – –

– – –

One-layer ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.8999 0.9013 0.9027

0.9717 0.9724 0.9721

0.9857 0.9874 0.9870

Two-layer ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

– – 0.9400

– – 0.9858

– – 0.9948

Corrected ES-NMM

MⅠ MⅠ MⅠ

2.0 2.5 3.0

0.9630 0.9644 0.9659

0.9911 0.9919 0.9915

0.9950 0.9968 0.9964

Standard NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

– – –

– – –

– – –

One-layer ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

0.9813 0.9632 0.9390

0.9661 0.9589 0.9695

1.0120 1.0248 1.0069

Two-layer ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

– – 0.9428

– – 0.9797

– – 1.0059

Corrected ES-NMM

MⅡ MⅡ MⅡ

2.0 2.5 3.0

0.9865 0.9785 0.9503

0.9757 0.9684 0.9794

1.0142 1.0251 1.0067

h

b

Different methods

Different methods

θ

Mesh

(b)

𝜃

a/b = 0.6

a/b = 0.7

MⅠ

MⅡ

MⅠ

MⅡ

Standard NMM

0° 15° 30° 45°

0.9841 0.9819 0.9784 0.9854

– 1.0081 0.9812 0.9927

0.9810 0.9860 0.9837 0.9832

– 1.0262 1.0006 0.9964

One-layer ES-NMM

0° 15° 30° 45°

0.9890 0.9855 0.9819 0.9897

– 1.0057 0.9790 0.9941

0.9895 0.9926 0.9921 0.9866

– 1.0404 1.0081 1.0062

Two-layer ES-NMM

0° 15° 30° 45°

1.0009 0.9986 1.0038 1.0032

– 0.9945 1.0152 1.0158

1.0010 1.0004 0.9980 1.0012

– 0.9893 1.0042 1.0043

Corrected ES-NMM

0° 15° 30° 45°

0.9976 0.9952 0.9937 0.9953

– 1.0031 0.9921 1.0169

0.9995 0.9962 0.9994 0.9899

– 1.0016 0.9991 0.9891

h = 1 mm h = 0.5 mm h = 0.25 mm (187 PPs) (637 PPs) (2307 PPs)

5.3. A plate with two cracks emanating from a hole A finite plate with two cracks emanating from a hole subjected to a uniform tensile force 𝜎 =1 Pa is considered, as shown in Fig. 13(a). In the computation, the following dimensions are used, h = 2b = 4 mm, r = 0.25b = 0.5 m. The reference solutions of SIFs depend on a and 𝜃, and can be found in reference [84]. Fig. 13(b) gives the mesh for this problem with a = 0.7b and 𝜃=45°. Influences of a and 𝜃 on SIFs are investigated. As for a, a/b = 0.6 and a/b = 0.7 are considered. As for 𝜃, four values, 0°, 15°, 30° and 45° are used. The same mesh is employed for different a/b and different 𝜃, as shown in Fig. 13(b). Note that different crack dimensions lead to slightly different numbers of PPs and MEs. The normalized SIFs of tip A and B are, respectively, listed in Tables 5 and 6. Similar conclusions, such as corrected ES-NMM performs best, ES-NMM overcome standard NMM, can be drawn.

“-” represents that the smoothing domains do not exist for some or all the elements involved in the interaction integrals, thus the SIFs cannot be computed.

411

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

Table 6 Normalized SIFs at crack tip B for a plate with two cracks emanating from a hole. Different methods

𝜃

a/b = 0.6

a/b = 0.7

MⅠ

MⅡ

MⅠ

MⅡ

0° 15° 30° 45°

0.9841 0.9827 0.9799 0.9867

– 0.9971 0.9872 0.9836

0.9810 0.9796 0.9804 0.9827

– 1.0065 0.9817 0.9997

One-layer ES-NMM

0° 15° 30° 45°

0.9890 0.9899 0.9845 0.9894

– 0.9902 0.9696 0.9760

0.9895 0.9942 0.9860 0.9877

– 1.0045 0.9831 1.0075

Two-layer ES-NMM

0° 15° 30° 45°

1.0009 1.0042 0.9998 0.9988

– 1.0190 0.9942 1.0188

1.0010 0.9988 1.0014 0.9989

– 0.9798 1.0043 1.0039

0° 15° 30° 45°

0.9976 1.0001 0.9928 0.9977

– 0.9929 0.9901 0.9934

0.9995 0.9952 0.9944 0.9924

– 0.9979 0.9915 0.9972

Standard NMM

Corrected ES-NMM

(a) Step 0

(b) Step 1

(c) Step 3

(d) Step 6

0.13 F

F 61 mm

224 mm

(e) Step 9

61 mm

82 mm

458 mm

458 mm

Fig. 15. The crack path for the Arrea–Ingraffea beam by corrected ES-NMM.

(a)

σ

(b)

h1

Fig. 14. (a) The geometry of Arrea–Ingraffea beam and (b) the mesh with 3479 PPs and 6679 MEs.

a 5.4. Arrea–Ingraffea beam In the first three examples, ES-NMMs have shown their ability to accurately calculate the SIFs of different problems. From this section, several crack propagation problems are used to verify the effectiveness and robustness of corrected ES-NMM, which is the most excellent one. The Arrea–Ingraffea beam [85] is considered as the first crack propagation problem. The geometry and the boundary conditions are schematically shown in Fig. 14(a). The mesh with 3479 PPs and 6679 MEs is shown in Fig. 14(b). In the computation, the Young’s modulus E = 28 000 MPa, the Poisson ratio v = 0.18, and plane strain is assumed. The prediction of crack path by the corrected ES-NMM is plotted in Fig. 15. The crack curves nicely, and becomes parallel to the loading direction, which agrees well with the experiment [85].

r

H

b1

b2

h2

5.5. Crack growth in a plate with a circular hole Fig. 16. (a) The geometry of a plate with a circular hole and (b) the mesh with 2322 PPs and 4380 MEs.

In this example, a plate with a circular hole is considered, as shown in Fig. 16(a). The plate is fixed at the bottom and subjected to a uniform 412

F. Liu, K. Zhang and D. Xu

(a) Step 0

Engineering Analysis with Boundary Elements 113 (2020) 402–415

(b) Step 1

(c) Step 3

Fig. 18. Crack growth path in a plate with a circular hole: (a) predicted by the Trig3-CNS (NMM) [56]; (b) predicted by the extended meshfree Galerkin radial point interpolation method [86].

(d) Step 5

(e) Step 7

(f) Step 9

Fig. 17. The crack path for a plate with a circular hole by corrected ES-NMM.

tensile load 𝜎 at the top. The geometric dimensions are taken as H = 3 m, h1 =1.2 m, h2 =1.5 m, b1 =0.7 m, b2 =0.3 m, initial crack length a = 0.1 m, the radius of the hole r = 0.2 m. The material parameters in the computation are as follows, 𝜎 =5000 Pa, the Young’s modulus E = 3.0 × 107 Pa, the Poisson’s ratio v = 0.2, and the plane strain condition is assumed. Fig. 16(b) describes the initial mesh with 2322 PPs and 4380 MEs. The crack growth paths obtained by corrected ES-NMM are shown in Fig. 17. It can be found that as the crack tip approaches the hole, it turns toward the hole and finally collapses with the hole. This is quite consistent with the results from Trig3-CNS (NMM) [56] and the extended meshfree Galerkin radial point interpolation method [86], see Fig. 18. Fig. 19 shows the load-displacement curves for three meshes (with 1532, 2322 and 3981 PPs, respectively). The curves are plotted following the procedure in the reference [12], where the fracture toughness condition (set to 1 kN m−3/2 here) as well as the equilibrium condition is satisfied at every step by scaling the magnitude of the load. In the load-displacement curves, the displacement is calculated by averaging the displacements of the upper side of the plate, while the load is calculated after the scaling. The Load-Displacement curves exhibit a typical snap-back response similar to [12], and indicate mesh-independence of the present method.

Fig. 19. Load–displacement curves of three meshes for a plate with a circular hole.

the initial mesh with 2806 PPs and 5342 MEs. In the computation, the parameters are taken as F = 100 N, the Young’s modulus E = 28 × 106 Pa, the Poisson’s ratio v = 0.18, and the plane strain condition is assumed. The paths of crack growth are plotted in Fig. 21, from which we can find that two cracks symmetrically propagate in the direction of maximal load on the opposite side of the beam. The predicated crack path agrees well with results by Trig3-CNS (NMM) [56] and MLS-based NMM with structured mesh refinement [60], see Fig. 22.

5.6. Crack growth in a beam under four-point shear loading As the last example, a pre-cracked beam subjected to four-point shear loading is considered. This example is used to validate the effectiveness of corrected ES-NMM for multiple crack problem. The geometry and the boundary conditions are shown in Fig. 20(a). Fig. 20(b) describes 413

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

6. Conclusion

F

In this study, the edge-based smoothed NMM is extended to linear elastic fracture problems. Blending elements related problems are also discussed by considering different enrichment schemes. Some important observations from this work are as follows:

40 mm 40 mm 120 mm

0.2F

320 mm

80 mm 80 mm

(1) Compared to the standard NMM, ES-NMMs produce better accuracy and slightly higher convergence rate without increasing DOFs, due to the use of edge-based strain smoothing technique. (2) Similar to the XFEM, blending elements induced problems exist in NMM. (3) Corrected ES-NMM, which is free from blending elements induced problems, performs best among the ES-NMMs. (4) Although corrected ES-NMM can greatly improve the accuracy of NMM, it cannot achieve optimal convergence rates. To achieve optimal convergence rates, the geometrical enrichment may be needed. (5) Whether or not using strain smoothing technique in the interaction integral has limited influence on the SIFs. Due to the limitations when strain smoothing technique is used, interaction integral without strain smoothing technique is suggested. (6) In the crack propagation analysis, the crack growth paths obtained by corrected ES-NMM agree well with experiments or those predicted by other numerical methods. (7) Inheriting all the advantages from NMM, ES-FEM and corrected XFEM, corrected ES-NMM is an accurate method which is suitable for crack analysis.

320 mm

(a)

(b) Fig. 20. (a) The geometry of a beam under four-point shear loading and (b) the mesh with 2806 PPs and 5342 MEs.

Acknowledgments This work was supported by the National Key R&D Program of China (Grant No. 2018YFC0406800), the China Scholarship Council (Grant No. 201806255077), the National Natural Science Foundation of China (Grant No. 11602165, 51704211). References [1] Liu GR, Nguyen-Thoi T, Lam KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. J Sound Vib 2009;320(4–5):1100–30. [2] Fries TP. A corrected XFEM approximation without problems in blending elements. Int J Numer Methods Eng 2008;75(5):503–32. [3] Chan S, Tuba I, Wilson W. On the finite element method in linear fracture mechanics. Eng Fract Mech 1970;2(1):1–17. [4] Aoki S, Kishimoto K, Kondo H, Sakata M. Elastodynamic analysis of crack by finite element method using singular element. Int J Fract 1978;14(1):59–68. [5] Areias P, Reinoso J, Camanho P, de Sá JC, Rabczuk T. Effective 2D and 3D crack propagation with local mesh refinement and the screened Poisson equation. Eng Fract Mech 2018;189:339–60. [6] Blandford GE, Ingraffea AR, Liggett JA. Two-dimensional stress intensity factor computations using the boundary element method. Int J Numer Methods Eng 1981;17(3):387–404. [7] Saleh A, Aliabadi M. Crack growth analysis in concrete using boundary element method. Eng Fract Mech 1995;51(4):533–45. [8] Belytschko T, Lu Y, Gu L. Crack propagation by element-free Galerkin methods. Eng Fract Mech 1995;51(2):295–315. [9] Zhuang X, Augarde CE, Mathisen KM. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. Int J Numer Methods Eng 2012;92(11):969–98. [10] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45(5):601–20. [11] Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng 1999;46(1):131–50. [12] Zheng H, Liu F, Li CG. The MLS-based numerical manifold method with applications to crack analysis. Int J Fract 2014;190(1–2):147–66. [13] Tsay R, Chiou Y, Chuang W. Crack growth prediction by manifold method. J Eng Mech 1999;125(8):884–90. [14] Bouchard PO, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria. Comput Methods Appl Mech Eng 2003;192(35–36):3887–908. [15] Tang X-h, Wu S-c, Zheng C, Zhang J-h. A novel virtual node method for polygonal elements. Appl Math Mech 2009;30(10):1233–46. [16] Hussein A, Aldakheel F, Hudobivnik B, Wriggers P, Guidault P-A, Allix O. A computational framework for brittle crack propagation based on an efficient virtual element method. Finite Elem Anal Des 2019;159(2019):15–32.

Fig. 21. The crack path for a plate with a circular hole by corrected ES-NMM.

(a)

(b) Fig. 22. Crack growth path in a beam under four-point shear loading: (a) predicted by the Trig3-CNS (NMM) [56] and (b) predicted by MLS-based NMM with structured mesh refinement [60].

414

F. Liu, K. Zhang and D. Xu

Engineering Analysis with Boundary Elements 113 (2020) 402–415

[17] Duarte CA, Hamzeh O, Liszka T, Tworzydlo W. engineering. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Comput Methods Appl Mech Eng 2001;190(15–17):2227–62. [18] 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. [19] Song C, Wolf JP. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method. Comput Struct 2002;80(2):183–97. [20] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61(13):2316–43. [21] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. Int J Numer Methods Eng 2005;63(11):1559–82. [22] Zheng C, Wu S, Tang X, Zhang J. A meshfree poly-cell Galerkin (MPG) approach for problems of elasticity and fracture. Comput Model Eng Sci 2008:149–78. [23] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Comput Methods Appl Mech Eng 2010;199(37):2437–55. [24] Amiri F, Anitescu C, Arroyo M, Bordas SPA, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Comput Mech 2014;53(1):45–57. [25] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196(29–30):2777–99. [26] Fries T-P, Byfut A, Alizada A, Cheng KW, Schröder A. Hanging nodes and XFEM. Int J Numer Methods Eng 2011;86(4–5):404–30. [27] Chen L, Rabczuk T, Bordas SPA, Liu GR, Zeng KY, Kerfriden P. Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Comput Methods Appl Mech Eng 2012;209:250–65. [28] Chahine E, Laborde P, Renard Y. A non-conformal eXtended finite element approach: integral matching XFEM. Appl Numer Math 2011;61(3):322–43. [29] Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. Int J Numer Methods Eng 2005;64(3):354–81. [30] Chen J-W, Zhou X-P. The enhanced extended finite element method for the propagation of complex branched cracks. Eng Anal Bound Elem 2019;104:46–62. [31] Mousavi SE, Grinspun E, Sukumar N. Higher-order extended finite elements with harmonic enrichment functions for complex crack problems. Int J Numer Methods Eng 2011;86(4–5):560–74. [32] Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Methods Eng 2000;48:1741–60. [33] Feng SZ, Bordas SPA, Han X, Wang G, Li ZX. A gradient weighted extended finite element method (GW-XFEM) for fracture mechanics. Acta Mech 2019;230(7):2385–98. [34] Xu Q, Chen J-y, Li J, Xu G. Study on the element with the hole and crack. Acta Mech 2013;225(7):1915–30. [35] Chau-Dinh T, Zi G, Lee PS, Rabczuk T, Song JH. Phantom-node method for shell models with arbitrary cracks. Comput Struct 2012;92-93:242–56. [36] Ghorashi SS, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Comput Struct 2015;147:138–46. [37] Shi GH, Goodman RE. Two dimensional discontinuous deformation analysis. Int J Numer Anal Methods Geomech 1985;9(6):541–56. [38] Shi GH. Manifold method of material analysis. In: Proceedings of the 9th army conference on applied mathematics and computing. Minneapolis, Minnesota: Army Research Office; 1991. p. 57–76. [39] Yang YT, Tang XH, Zheng H, Liu QS, He L. Three-dimensional fracture propagation with numerical manifold method. Eng Anal Bound Elem 2016;72(Supplement C):65–77. [40] Yang YT, Tang XH, Zheng H, Liu QS, Liu ZJ. Hydraulic fracturing modeling using the enriched numerical manifold method. Appl Math Model 2018;53:462–86. [41] Zheng H, Yang Y, Shi G. Reformulation of dynamic crack propagation using the numerical manifold method. Eng Anal Bound Elem 2019;105:279–95. [42] Ma GW, An XM, Zhang HH, Li LX. Modeling complex crack problems using the numerical manifold method. Int J Fract 2009;156(1):21–35. [43] Wu Z, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39:38–53. [44] Zhang HH, Li LX, An XM, Ma GW. Numerical analysis of 2-D crack propagation problems using the numerical manifold method. Eng Anal Bound Elem 2010;34(1):41–50. [45] Zhang HH, Liu SM, Han SY, Fan LF. The numerical manifold method for crack modeling of two-dimensional functionally graded materials under thermal shocks. Eng Fract Mech 2019;208:90–106. [46] Zheng H, Xu DD. New strategies for some issues of numerical manifold method in simulation of crack propagation. Int J Numer Methods Eng 2014;97(13):986–1010. [47] Wu ZJ, Sun H, Wong LNY. A cohesive element-based numerical manifold method for hydraulic fracturing modelling with voronoi grains. Rock Mech Rock Eng 2019;52(7):2335–59. [48] Zheng H, Liu F, Li CG. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method. Appl Math Model 2015;39(2):794–808. [49] Jiang QH, Deng SS, Zhou CB, Lu WB. Modeling unconfined seepage flow using three-dimensional numerical manifold method. J Hydrodyn Ser B 2010;22(4):554–61. [50] Zheng H, Liu ZJ, Ge XR. Numerical manifold space of Hermitian form and application to Kirchhoff’s thin plate problems. Int J Numer Methods Eng 2013;95(9):721–39. [51] Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. Int J Rock Mech Min Sci 2011;48(6):964–75. [52] He L, An XM, Ma GW, Zhao ZY. Development of three-dimensional numerical manifold method for jointed rock slope stability analysis. Int J Rock Mech Min Sci 2013;64:22–35.

[53] Xu Y, Yu C, Liu F, Liu Q. A coupled NMM-SPH method for fluid-structure interaction problems. Appl Math Model 2019;76:466–78. [54] Xu DD, Yang YT, Zheng H, Wu AQ. A high order local approximation free from linear dependency with quadrilateral mesh as mathematical cover and applications to linear elastic fractures. Comput Struct 2017;178:1–16. [55] Liu Z, Zhang P, Sun C, Liu F. Two-dimensional Hermitian numerical manifold method. Comput Struct 2020;229. doi:10.1016/j.compstruc.2019.106178. [56] Yang YT, Zheng H. A three-node triangular element fitted to numerical manifold method with continuous nodal stress for crack analysis. Eng Fract Mech 2016;162:51–75. [57] Zheng H, Liu F, Du XL. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Comput Methods Appl Mech Eng 2015;295:150–71. [58] Liu ZJ, Zheng H. Two-dimensional numerical manifold method with multilayer covers. Sci China Technol Sci 2016;59(4):515–30. [59] Liu ZJ, Zheng H, Sun C. A domain decomposition based method for two-dimensional linear elastic fractures. Eng Anal Bound Elem 2016;66:34–48. [60] Liu F, Xia KW. Structured mesh refinement in MLS-based numerical manifold method and its application to crack problems. Eng Anal Bound Elem 2017;84:42–51. [61] Yang SK, Ma GW, Ren XH, Ren F. Cover refinement of numerical manifold method for crack propagation simulation. Eng Anal Bound Elem 2014;43:37–49. [62] Yu CY, Liu F, Xu Y. An h-adaptive numerical manifold method for solid mechanics problems. Sci China Technol Sci 2018;61(6):923–33. [63] Liu F, Yu CY, Yang YT. An edge-based smoothed numerical manifold method and its application to static, free and forced vibration analyses. Eng Anal Bound Elem 2017;86:19–30. [64] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. Int J Numer Methods Eng 2001;50(2):435–66. [65] Nguyen-Thoi T, Vu-Do HC, Rabczuk T, Nguyen-Xuan H. A node-based smoothed finite element method (NS-FEM) for upper bound solution to visco-elastoplastic analyses of solids using triangular and tetrahedral meshes. Comput Methods Appl Mech Eng 2010;199(45–48):3005–27. [66] Nguyen-Thoi T, Liu GR, Lam KY, Zhang GY. A face-based smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements. Int J Numer Methods Eng 2009;78(3):324–53. [67] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Comput Mech 2006;39(6):859–77. [68] Jiang Y, Tay T, Chen L, Zhang Y. Extended finite element method coupled with face-based strain smoothing technique for three‐dimensional fracture problems. Int J Numer Methods Eng 2015;102(13):1894–916. [69] Francis A, Ortiz-Bernardin A, Bordas SPA, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. Int J Numer Methods Eng 2017;109(9):1263–88. [70] Gui-rong L, Gui-yong Z. Smoothed point interpolation methods: g space theory and weakened weak forms. World Scientific; 2013. [71] Zhang G, Wang Y, Jiang Y, Jiang Y, Zong Z. A combination of singular cell-based smoothed radial point interpolation method and fem in solving fracture problem. Int J Comput Methods 2018;15(08):1850079. [72] Liu G. A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods. Int J Comput Methods 2008;5(02):199–236. [73] Liu GR, Chen LEI, Li M. S-Fem for fracture problems, theory, formulation and application. Int J Comput Methods 2014;11(03):1343003. [74] Liu GR, Nguyen TT. Smoothed finite element methods. CRC Press LLC; 2010. [75] Liu G-R. The smoothed finite element method (S-FEM): a framework for the design of numerical models for desired solutions. Front Struct Civil Eng 2019;13(2):456–77. [76] Zeng W, Liu GR. Smoothed finite element methods (S-FEM): an overview and recent developments. Arch Comput Methods Eng 2018;25(2):397–435. [77] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee CK, Zi G, Zhuang X, et al. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. J Appl Math 2013;2013:1–12. [78] Liu P, Bui TQ, Zhang C, Yu TT, Liu GR, Golub MV. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Comput Methods Appl Mech Eng 2012;233-236:68–80. [79] Bordas SPA, Natarajan S, Kerfriden P, Augarde CE, Mahapatra DR, Rabczuk T, et al. On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM). Int J Numer Methods Eng 2011;86(4–5):637–66. [80] An X, Ma G, Cai Y, Zhu H. A new way to treat material discontinuities in the numerical manifold method. Comput Methods Appl Mech Eng 2011;200(47–48):3296–308. [81] Erdogan F, Sih G. On the crack extension in plates under plane loading and transverse shear. J Basic Eng 1963;85(4):519–25. [82] Liu F, Zhang K, Liu Z. Three-dimensional MLS-based numerical manifold method for static and dynamic analysis. Eng Anal Bound Elem 2019;109:43–56. [83] Béchet É, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the XFEM for stress analysis around cracks. Int J Numer Methods Eng 2005;64(8):1033–56. [84] Establishment CA. Handbook of the stress intensity factor. Beijing: Science Press; 1993. [85] Arrea M. Mixed-mode crack propagation in mortar and concrete. Cornell University, Jan.; 1982. [86] Nguyen NT, Bui TQ, Zhang C, Truong TT. Crack growth modeling in elastic solids by the extended meshfree Galerkin radial point interpolation method. Eng Anal Bound Elem 2014;44:87–97.

415