Structured mesh refinement in MLS-based numerical manifold method and its application to crack problems

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Contents lists available at ScienceDirect

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

Structured mesh refinement in MLS-based numerical manifold method and its application to crack problems Feng Liu a,b, Kaiwen Xia a,c,∗ a

State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China c Department of Civil Engineering, University of Toronto, Toronto, Ontario M5S 1A4, Canada b

a r t i c l e

i n f o

Keywords: Numerical manifold method Moving least squares Stress intensity factor Crack propagation Eccentric domain of influence Structured mesh refinement

a b s t r a c t The numerical manifold method (NMM) is suitable for both continuity and discontinuity problems due to its two cover systems, namely the mathematical cover and the physical cover. However, the use of these two cover systems makes it difficult to implement local mesh refinement in conventional NMM. In this paper, a structured mesh refinement is proposed for crack problems within the framework of MLS-based NMM. In order to balance both accuracy and efficiency, an eccentric domain of influence is suggested near the transition region, where the domain size is determined by surrounding elements. The visibility criterion used in element-free Galerkin method is introduced to describe the discontinuity near the crack tip. The main advantage of the proposed method is that it can solve complex crack problems straightforwardly without introducing any extra degree of freedom. Several examples are tested to demonstrate the effectiveness of the proposed method for crack problems. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The numerical simulation of crack problem is vital in solid mechanics and one of the most difficult part of crack simulation is the handling of the singularity at the crack tip. Great efforts have been devoted to overcome this problem over years. Local mesh refinement is the common way to resolve this problem. A variety of local remeshing techniques in finite element method (FEM) have been developed, such as interactive approach for local remeshing [1], automatic remeshing technique [2], adaptive refinement procedure [3], advanced remeshing technique [4], and so on. These remeshing techniques along with other mesh refinement methods [5] achieved excellent results in crack propagation simulations. As meshfree methods have been fast developing, refinements or adaptive procedures in meshfree methods have also been proposed. Belytschko et al. [6] suggested a nodal refinement using a star-shaped nodal array at the crack tip in element-free Galerkin method (EFG) [7]; Liu and Tu [8] proposesed an adaptive procedure based on background cells; Rabczuk and Belytschko [9] implemented adaptivity in a structured model; Schweitzer [10] present a refinement scheme in particle–partition of unity method. Another way to handle the singularity at crack tips is to enrich the standard displacement-based approximation with the near tip asymptotic fields, which was first used by the extended finite element method

∗

(XFEM) [11]. In XFEM, the representation of crack is independent of mesh, so that no remeshing is required in crack propagation. Consequently, significant progress has been made for crack simulation with XFEM [12–15]. The numerical manifold method (NMM), originally proposed by Shi [16], is a combined continuum–discontinuum method. On one hand, NMM can solve continuum problems in the form of FEM; on the other hand, it can deal with discontinuum problems just as Discontinuous Deformation Analysis (DDA) [17] does. Hence, NMM is suitable for problems from continuum to discontinuum. In some studies of crack problems with NMM, crack tips were restricted on the edges of manifold elements [18–20]. Two simple refinements near the crack tip were proposed by Tsay et al. [21,22]. Yang et al. [23] discussed a cover refinement method in detail for the simulation of crack propagation in brittle materials. However, these refinements lead to non-regular meshes which are not consistent with original regular meshes. Recently, Liu and Zheng [24] propose a multilayer covers with NMM for complex problems with heterogeneity or nonlinearity. While this method implements local refinement on uniform meshes, it cannot be readily used for crack propagation. Inspired by XFEM, more and more researchers enrich the NMM with asymptotic field to model crack propagation [25–32]. Using this technique, the discontinuities are completely independent of the mesh, and thus remeshing is avoided. However, it suffers from ill-conditioning and low-convergence rate due to the blending elements which is similar

Corresponding author at: State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China. E-mail address: [email protected] (K. Xia).

http://dx.doi.org/10.1016/j.enganabound.2017.08.004 Received 21 April 2017; Received in revised form 9 July 2017; Accepted 7 August 2017 0955-7997/© 2017 Elsevier Ltd. All rights reserved.

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Each PP inherits a partition of unity function 𝜑i (z) from the corresponding MP. By multiplying local approximation function uI (z) with partition of unity 𝜑i (z), the global approximation is obtained,

to XFEM [33–35]. Moreover, the enriched degrees of freedom are lack of physical meaning. In this paper, we focus on structured mesh refinement in NMM. Structured mesh refinement has already been presented in the finite element context [36]. Rabczuk and Belytschko proposed an adaptive scheme for structured meshfree particle methods in 2D and 3D. Following this technique, a structured mesh refinement in NMM is proposed for crack problems. To avoid the incompatibility between the meshes in transition region of refinement, the MLS-based NMM is adopted [28,37]. Dolbow and Belytschko [38] found that when integration cells align with the shape function supports (or to say nodal influence domains) in EFG, integral precision can be improved significantly. Following this principle, the size and shape of nodal influence domains are carefully designed, and an eccentric weight function is proposed to balance both accuracy and efficiency. The proposed method is suitable for complex crack problems and no extra degree of freedom is introduced. The paper is organized as follows. In the next section, we briefly recall the fundamental of NMM. Section 3 introduces the structured mesh refinement strategy and the eccentric weight function. In Section 4, we discuss several aspects for crack modeling by MLS-based NMM, including the representation of discontinuities, the variational formulation, the integration scheme, evaluation of stress intensity factor and fracture criterion and so on. Some typical numerical examples are analyzed in Section 5, followed by some concluding remarks in the last section.

𝒖(𝒛) =

𝑖

To deal with various kinds of problem domains, NMM introduces two covers, called the mathematical cover (MC) and the physical cover (PC). MC consists of a set of simply-connected domains, each of which is called a mathematical patch (MP). Since the MP does not need to coincide with the geometrical details of problem domain Ω, there are various shapes to choose from. The only requirement is that all the MPs together should cover Ω completely. Physical patches (PPs) are obtained by the intersection of each MP and the problem domain one by one. All the PPs constitute the PC. It is obvious that PC relies on the physical features and covers Ω accurately. Then a manifold element (ME) can be formed by the common area of several adjacent PPs. Manifold elements vary in shape and size, and are mainly used for integration. The detailed description of NMM can be found in literatures [27,39]. Let z represents the position vector (x, y). Associated with the i-MP Mi , a partition of unity 𝜙i (z) is defined: (1.1)

0 ≤ 𝜙𝑖 (𝒛) ≤ 1, if 𝒛 ∈ 𝑀𝑖 ;

(1.2)

𝑛 ∑ 𝑖=1

𝜙𝑖 (𝒛) = 1, if 𝒛 ∈ Ω

𝜙𝑖 (𝒛)𝒖𝐼 (𝒛),𝑓 𝑜𝑟𝒛 ∈ Ω

(4)

3. Structured mesh refinement and eccentric domain of influence In this section, the structured mesh refinement strategy is first introduced, then an eccentric domain of influence on the structured mesh refinement and related content are present.

(1.3) 3.1. Structured mesh refinement strategy

here, we do not call 𝜙i (z) ‘weight function’ as usually does in NMM, since the term ‘weight function’ is used for constructing domain of influence of nodes in later section. Over each PP, a local approximation function (also named as cover function), denoted by uI (z), is expressed as 𝒖𝐼 (𝒛) = 𝑻 (𝒛)𝒅 𝐼

𝐼

In theory, the MC can be constructed arbitrarily. However, almost all the applications of NMM have selected the regular finite element mesh as the MC (denoted by FE-based NMM), where a MP is composed of all elements linked at a given node. As the discontinuity updates, the PCs have to be cut frequently to fit the new features of the problem. On one hand, this could be quite tedious during the simulation of crack growth. On the other hand, the precision of NMM with structured triangular mesh is far from satisfactory, because the strain distribution is constant in each manifold element. Since the MLS shape functions, widely used in meshfree method, satisfy the partition of unity, an alternative choice to construct the MC is to use the domains of influence of these MLS-nodes (denoted by MLSbased NMM). The domain of influence of a node is then used to form PPs the same way as a MP does. In order to get a better understanding of the differences between the FE-based NMM and MLS-based NMM, Fig. 1 shows a comparison of the generation of MP and PP. In Fig. 1(a), a rectangular region is covered by regular finite element mesh. The MP of the red node is the red shaded area, which generate a PP with the same region, while the MP of the blue node is a regular hexagon, which generate a PP of the blue shaded pentagon. In Fig. 1(b), the rectangular region is covered by the rectangular domains of influence of several nodes. The MPs and PPs of the yellow and green nodes are formed in a similar manner. In Fig. 1(a), the approximation is implemented based on the finite element mesh. However, in Fig. 1(b), the approximation is accomplished with the scattered nodes using MLS method, where a finite element mesh is not necessary. It should be noted that, in EFG, the nodes are generally arranged inside the domain, but in MLS-based NMM, nodes can be arranged outside the domain only if their domains of influence overlap the domain, as shown in Fig. 1(b). This can be quite useful for problems with sharp areas [37]. Compared to FE-based NMM, the MLS-based NMM is more like a meshfree method in some sense, and can simplify the geometric operations and achieves better accuracy at the cost of some computational efficiency. It can be said that MLS-based NMM inherits the advantages of both NMM and EFG. More details about MLS-based NMM can be found in [28].

2. Brief review on MLS-based NMM

𝜙𝑖 (𝒛) = 0, if 𝒛 ∉ 𝑀𝑖 ;

𝑛 ∑ 𝑚 ∑

Although MLS approximation is known for its ability to deal with a random cloud of particles, structured mesh has its unique advantages. On one hand, the MC is independent of the physical features in MLSbased NMM, indicating that structured mesh is applicable to all kinds of problems. There is no reason to choose scatter nodes when regular nodes are feasible. On the other hand, structured mesh leads to better accuracy in general case [9,24]. However, due to the presence of stress singularity, uniform meshes are not optimal for crack problems, and a mesh refinement strategy may be needed. One of the most popular algorithms to refine the meshes is based on a recursive decomposition using quadtree meshes [40,41]. Hanging nodes are generated if the adjacent elements are not in the same size. Usually hanging nodes exist in the transition region of the quadtree meshes. In mesh-based method, proper handling of hanging nodes is a

(2)

here, dI is the degrees of freedom vector for the PP, and T(z) is the matrix of basis functions, which are able to reflect the local behaviors of the solution. For two-dimensional elasticity problems, the most conventional choice of T(z) is the polynomial basis with the form [ ] 1 0 𝑥 0 𝑦 0 ... 𝑻 (𝒛) = (3) 0 1 0 𝑥 0 𝑦 ... In this paper, we take T(z) as the second order identity matrix, which means the constant basis is used. 43

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Fig. 1. The sketch of MPs and PPs in (a) FE-based NMM, and (b) MLS-based NMM.

Fig. 2. Structured mesh refinement strategy.

core problem using quadtree meshes. Various techniques have been developed to overcome this issue and the reader is referred to [42] for details. Quadtree meshes can also be used in meshfree method to implement structured nodes refinement [43,44]. In meshfree method, the trouble of hanging nodes does not exist thanks to the meshfree approximations. In this paper, MLS-based NMM adopts similar mesh refinement strategy. Due to the meshfree feature of MLS-based NMM, the trouble of hanging nodes does not occur anymore. The process of structured mesh refinement is as follows. First, a uniform mesh covering the whole problem domain is constructed, each small square of which is termed as a pixel. This is similar to the construction of background mesh in EFG for the purpose of integration. To align the domain of influence with the mesh, the rectangular weight function is adopted, and the size is twice of the span of nodes. Then, a pixel is divided into four small pixels every time when refinement is needed. Other than refining the mesh layer-by-layer as in [9], more layers are involved here. The refinement strategy is illustrated in Fig. 2. The hierarchical mesh is required to be changed gradually, which means no pixel is allowed to share more than two edges with a neighbour pixel.

3.2. Eccentric domain of influence When the crack tip locates at a node, a typical four-layer hierarchical mesh is generated as shown in Fig. 3. Suppose that the nodal influence domains are twice of the nodal spacing. Without regard to the mesh refinement, the domain of influence of the yellow node (the yellow dotted line in Fig. 3) covers too many pixels in the refinement region. So the domains of influence of nodes near the transition region have to be adjusted to reduce the inappropriate selection of nodes. To this end, eccentric domain of influence can be designed. Since the domain of influence is entirely dependent on its corresponding weight function, what we have to do is to construct appropriate eccentric weight function. In [44], an eccentric weight function is constructed by introducing a constant parameter 𝛼. In this paper, the eccentric weight function is extended from the rectangular weight function used in EFG without extra parameter. The eccentric rectangular weight function WI (z) of node

Fig. 3. Four-layer hierarchical mesh refinement and the eccentric domains of influence for different nodes. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

I with coordinates (xI , yI ) is taken as 𝑊𝐼 (𝒛) = 𝑊 (𝑟𝑥 )𝑊 (𝑟𝑦 ) where rx , ry are defined as follow: { (𝑥 − 𝑥𝐼 )∕𝑑 𝑥+ if 𝑥 ≥ 𝑥𝐼 𝑟𝑥 = (𝑥𝐼 − 𝑥)∕𝑑 𝑥− if 𝑥 < 𝑥𝐼 44

(5)

(6a)

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

discontinuity properly. Before dealing with such a case here, we will look into the strategy used in meshfree method. To model the strong discontinuities, the visibility method [6,45] is first proposed in EFG. In the visibility method, the crack boundary is assumed to be opaque. The domain of influence is considered as the field of vision at a node, as shown in Fig. 6(a). In MLS-based NMM, a similar technique can be adopted. Discarding the invisible gray region, the influence of domain of node I is the polygonal domain at the top right, as shown in Fig. 6(b). This technique is also termed “visibility method”. As an alternative, the gray region can be taken as a new independent PP. Such a practice is denoted by “partition criterion”. The performances of these two schemes are compared in Section 5.2. Since two PPs can be obtained from node J and the discontinuities can be naturally captured in NMM, no extra treatment is needed.

Fig. 4. The length of influence domain of node I at different directions.

𝑟𝑦 =

{ (𝑦 − 𝑦𝐼 )∕𝑑 𝑦+ (𝑦𝐼 − 𝑦)∕𝑑 𝑦−

if 𝑦 ≥ 𝑦𝐼 if 𝑦 < 𝑦𝐼

4.2. Variational formulation and its discretization

(6b)

By using structured mesh, NMM allows the mesh not to match the boundaries. Meanwhile, the MLS shape function is lack of Kronecker Delta property. As a result, the essential boundary condition cannot be enforced directly as in FEM. We can apply the Lagrange multiplier method or penalty function method to enforce the essential boundary condition. When the penalty function method is used, the penalized functional writes

with dx + , dx − , dy + , dy − the length of domain of influence of node I at different directions (shown in Fig. 4); and the cubic spline form of W(r) is written by ⎧2∕3 − 4𝑟3 + 4𝑟3 , ⎪ 𝑊 (𝑟) = ⎨4∕3 − 4𝑟 + 4𝑟2 − 4∕3𝑟3 , ⎪0 , ⎩

𝑟 ≤ 0.5 0.5 < 𝑟 ≤ 1 𝑟>1

(7)

The dimension of weight function defines the dimension of domain where WI (z) ≠ 0. Based on the neighbour pixels, the dimension of weight function for each node can be determined. First, the pixels associated with the node are found out. Then the pixels are used to determine the spans of nodes in four directions (see Fig. 4). When more than one span is obtained in one direction, the small one is selected. At last, dx + , dx − , dy + , dy − are obtained as twice of the spans, respectively. Three typical cases are plotted in Fig. 3, the red rectangle, the green rectangle and the blue rectangle. The eccentric weight function of the red node and its derivative are plotted in Fig. 5. We can see that the eccentric weight function and its derivative are smooth enough. It should be pointed out that the domain of influence of a node only covers the pixels of its neighboring layer using such a strategy here. As the layer number of refinement increases, the advantage becomes more prominent. The above mesh refinement strategy along with the eccentric domain of influence has the following advantages: (1) each pixel is influenced by nodes of its neighboring layer, avoiding overmany nodes being involved in the construction of MLS shape function; (2) furthermore, it can guarantee the domain of influence align with the integration domains to the greatest extent. These two advantages indicate that the proposed method can obtain high accuracy without increasing too much computation amount.

𝜋(𝒖) =

1 𝐓 𝜺 𝑫 𝜺𝑑Ω − 𝒖𝑇 𝒃 𝑑Ω − 𝒖𝐓 𝒕̄𝑑𝑆 ∫Ω 2 ∫Ω ∫Γ𝐭 +

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

(8)

where 𝜺 is the strain verctor, D the elasticity matrix, b the body force, and Γu represents the displacement boundary on which the displacement is given by 𝒖̄ ; Γt the stress boundary on which the traction vector 𝒕̄ is loaded; kp the penalty factor. By substituting Eq. (4) into Eq. (8) and putting 𝛿𝜋(u) = 0 (here 𝛿 is the variational symbol), we have the MLS-based NMM system of equations as follows: 𝑲𝑼 = 𝒒

(9)

where U is the vector of unknowns, K the global stiffness matrix, 𝑲=

∫Ω

𝑩 𝐓 𝑫 𝑩 𝐝Ω + 𝑘𝐩

∫Γ𝐮

𝑵 𝑇 𝑵 𝑑Ω

(10)

with 𝑩 = 𝑳Td 𝑵 [ 𝑳d ≡

4. Modeling crack problems by MLS-based NMM

𝜕 𝜕𝑥

0

(11) 0 𝜕 𝜕𝑦

𝜕 ] 𝜕𝑦 𝜕 𝜕𝑥

(12)

Several issues on crack modeling by MLS-based NMM are discussed in this section, including the representation of discontinuities in Section 4.1, variational formulation and its discretization in Section 4.2, integration scheme in Section 4.3 and the evaluation of stress intensity factors and fracture criterion in Section 4.4.

and q the generalized force vector,

4.1. Representation of discontinuities

In FE-based NMM, the simplex integration [16] over manifold elements can be used for integration. The simplex integration is analytical and only available when the integrands are polynomials. Since the MLS shape function is rational, numerical integration scheme has to be introduced to calculate Eqs. (10) and (13). The aforementioned pixels are used to execute the numerical integration. When a pixel locates outside the problem domain 𝛀, it should be ignored; when a pixel intersects with the boundary of 𝛀, the internal part is triangulated and 12 points Gauss quadrature is applied on these

𝒒=

∫Ω

𝑵 𝐓 𝒃𝑑Ω +

∫Γ𝐭

𝑵 𝐓 𝒕̄𝑑𝑆 + 𝑘𝑝

∫Γ𝑢

𝑵 𝑇 𝒖̄ 𝐝𝑆

(13)

4.3. Integration scheme

In fracture problems, the discontinuities along the crack should be accurately captured. By dividing the MPs into the PPs and defining the approximation over each PP independently, NMM can represent the discontinuities straightforward. This is a significant advantage comparing to XFEM, where Heaviside functions or discontinuous ‘junction’ functions are needed for representation of different kinds of discontinuities [12]. However, when a PP is cut partially, NMM cannot describe the 45

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Fig. 5. The eccentric weight function of the red node (in Fig. 3) and its derivative, (a) weight function; and (b) its partial derivative with regard to x. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. The visibility criterion (a) in EFG and (b) in present paper.

sub-triangles; when a pixel locates completely inside of 𝛀, 4 × 4 points Gauss quadrature is applied; when a pixel intersects with cracks, first it is cut into several parts by the cracks and each part is triangulated and 12 points Gauss quadrature is applied on these sub-triangles. A typical demonstration is shown in Fig. 7.

Fig. 7. A typical demonstration of the integration scheme in present paper.

4.4. Fracture criterion The direction of crack propagation is determined by the maximum circumferential tensile stress criterion [46], which advocates that the crack will propagate in the direction perpendicular to the maximum circumferential stress. The propagation angle 𝜃 c relative to the current crack segment is calculated as follow [46]: √ ( )2 ⎞ ⎛ 1 − 1 + 8 𝐾II ∕𝐾𝐼 ⎟ ⎜ 𝜃𝑐 = 2tan−1 ⎜ ( ) ⎟ 4 𝐾II ∕𝐾𝐼 ⎜ ⎟ ⎝ ⎠

Fig. 8. The polar coordinates of the crack tip and propagation angle 𝜃 c .

5. Numerical examples In order to verify the proposed mesh refinement algorithms, first we consider the problem of half space under concentrated force. Then several typical crack examples, including single and multiple cracks as well as propagating cracks, are tested to illustrate the effectiveness of the proposed method.

(14)

where KI and KII are stress intensity factors (SIF) of mode-I and mode-II, respectively, the polar coordinates of the crack tip and the propagation angle 𝜃 c are shown in Fig. 8. Several methods have been proposed for the calculation of SIFs such as virtual crack extension method [47], contour integral method [48], interaction integral method [11]. The domain form of interaction integral [11] is most widely used and is adopted in this paper.

5.1. Half space under concentrated vertical force Let us consider a concentrated vertical force P acting at the origin of a horizontal half space, as shown in Fig. 9. The analytical solutions to 46

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Fig. 9. Half space under concentrated vertical force and the computational region. Fig. 11. The relative errors in the displacement norm with different layer number of refinement.

bottom edges are displacement boundaries, which are constrained using the exact displacements above. It is obvious that stress concentration occurs at the origin due to the concentrated force F. So we implement mesh refinement only around the origin. Based on the initial 20 × 10 mesh, two layer refinement meshes are shown in Fig. 10. To study the effectiveness of the proposed method, the relative errors in the displacement norm is defined as √ ∫Ω (𝒖𝑒𝑥 − 𝒖𝑛𝑢𝑚 )𝑑Ω (16) 𝑒𝑑 = ∫Ω (𝒖𝑒𝑥 )𝑑Ω where the superscript “ex” and “num” represent the analytic solution and numerical solution, respectively. The relative errors in the displacement norm with different layer number of refinement are plotted in Fig. 11. We can see that as the layer number increases, the relative errors in the displacement norm can be reduced greatly. The deformed surface of the half space and a detail view are shown in Fig. 12. It is observed again that local refinements can improve the accuracy of displacements, especially for that near the loading point. To verify the accuracy and efficiency of the proposed refinement method, we also solve this problem using uniform refinement meshes. Fig. 13 displays the relative error in displacement norm along with the number of PPs. It is found that the proposed refinement method can obtained more accurate results with less number of PPs compared to uniform refinement method. The number of PPs increases much faster as the layer number of refinement increases using uniform refinement

Fig. 10. The meshes with two layer refinement based on the initial 20 × 10 mesh.

this problem can be found in [49], 𝑢𝑟 = − 𝑢𝜃 =

(1 − 𝑣)𝐹 2𝐹 cos 𝜃 ln 𝑟 − 𝜃 sin 𝜃 + 𝐼 cos 𝜃 𝜋𝐸 𝜋𝐸

(1 − 𝑣)𝐹 (1 + 𝑣)𝐹 2𝐹 sin 𝜃 ln 𝑟 − 𝜃 cos 𝜃 + sin 𝜃 − 𝐼 sin 𝜃 𝜋𝐸 𝜋𝐸 𝜋𝐸

(15a)

(15b)

where constant I can’t be determined due to rigid body displacement if no vertical constraint is imposed. The computational region is taken as a rectangular (10 × 5) and plane stress condition is assumed. The computational parameters are Young’s modulus E = 107 , Poisson’s ratio 𝜈 = 0.2. In the computation, the condition ur |r = 5, 𝜃 = 0 = 0 is used to determine constant I. The left, right and

Fig. 12. The deformed surface with different layer number of refinement: (a) global view; (b) local view near the loading point.

47

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

Fig. 13. The relative errors in the displacement norm with the number of PPs.

Fig. 15. Normalized KⅠ as a function of the layer number of refinement with different schemes of representation of discontinuity.

Fig. 16. A star-shaped crack in a square plate under biaxial tension.

Fig. 14. Single edge-cracked plate under uniaxial tension: (a) geometry and loading; (b) the mesh with three layer refinement.

An initial uniform mesh of 10 × 20 pixels is used. Subsequently, the mesh is locally refined for different times (mesh with three layer refinement is plotted in Fig. 14(b)). The effect of representation of discontinuity is also considered. In addition to the aforementioned schemes (denoted by “visibility criterion” and “partition criterion”), a third method is used for comparison, where the partially cut PPs are viewed as ordinary patches without special treatment (denoted by “without treatment”). The normalized value of KⅠ is plotted in Fig. 15 as a function of the layer number of refinement with different schemes of representation of discontinuity. The error decreases rapidly as the layer number increases, indicating that the refinement significantly improves the precision. It is obvious that “visibility criterion” works best, followed by “partition criterion”, and “without treatment” leads greater errors. The explanation for “visibility criterion” better than “partition criterion” may be that, in order to cut open the PP, false discontinuous line is introduced in “partition criterion” method. Next, the variations of SIFs with different crack length are tested. The layer number of refinement is fixed to 3 with 10 × 20 initial mesh, and “visibility criterion” is adopted. The normalized KⅠ of different crack

method. This illustrates the accuracy and efficiency of the proposed refinement method for stress concentration problems. 5.2. Edge crack Shown in Fig. 14(a) is a rectangular plate with an edge crack. The plate is loaded in tension at the top and bottom with 𝜎=1.The parameters include: L = 2, W = 1, the crack length a = 0.5, Young’s modulus E = 3 × 105 , Poisson’s ratio 𝜈 = 0.25, and the plane strain condition is assumed. The reference solution for this problem is given by Tada et al. [50] √ 𝐾I = 𝜎 𝜋𝑎𝐹 (𝑎∕𝑊 ) (17) where F(a/W) is a geometry correction faction: √ 2𝑊 𝜋𝑎 0.752 + 2.02(𝑎∕𝑊 ) + 0.37(1 − sin 𝐹 (𝑎∕𝑊 ) = tan ⋅ 𝜋𝑎 2𝑊 cos 2𝜋𝑎 𝑊

𝜋𝑎 3 ) 2𝑊

(18)

where a/W ∈ (0, 1). 48

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51 Table 1 Normalized KⅠ for edge crack with different crack length. a

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Numerical KⅠ Reference KⅠ Normalized KⅠ

0.6700 0.6702 0.9997

1.0890 1.0833 1.0053

1.6199 1.6068 1.0081

2.3798 2.3630 1.0071

3.5601 3.5426 1.0050

5.5704 5.5511 1.0035

9.4827 9.4545 1.0030

19.1802 19.0123 1.0088

Table 2 Number of MPs and PPs for different a/W. a/W

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

MP PP

2110 2299

2422 2671

2582 2857

2662 2977

2708 3059

2716 3071

2716 3078

2716 3113

2684 3123

Fig. 17. The mesh for star-shaped crack with four layer refinement for a/W = 0.2: (a) global view; (b) local view near crack tip B.

length are listed in Table 1. It can be observed that the accuracy of all results is excellent and all the relative errors are less than one percent. 5.3. A star-shaped crack in a square plate To show the ability of proposed method for more complex crack problems, we analyze a star-shaped crack subjected to bi-axial tension in plane strain condition, as shown in Fig. 16. The normalized stress intensity factors at tips A and B are defined as ( √ ) ( √ ) ( √ ) 𝐹𝐼𝐴 = 𝐾𝐼𝐴 ∕ 𝜎 𝜋𝑎 , 𝐹𝐼𝐵 = 𝐾𝐼𝐵 ∕ 𝜎 𝜋𝑎 , 𝐹II𝐵 = 𝐾II𝐵 ∕ 𝜎 𝜋𝑎 (19) Fig. 18. Sketch of the beam under four-point bending.

In this example, the plate dimension is fixed to be W = 5, the material constants are Young’s modulus E = 200,000 and Poisson’s ratio v = 0.3, the bi-axial tension 𝜎 = 1. We investigate the influence of finiteness of plate on SIFs for ratios of a/W from 0.1 to 0.9. A uniform mesh of 25 × 25 with 4 layer refinement is used for each case. A typical mesh for a/W = 0.2 and the detail view of crack tip B are shown in Fig. 17. Number of MPs and PPs for different ratios of a/W are listed in Table 2. The computed SIF results are summarized and compared with three reference solutions, respectively, from Cheung et al. [51], Daux et al. [52], and Ma et al. [31] in Table 3. Results are in good agreement with the existing solutions. It can be seen that the proposed method can handle complex crack problems in a natural way, no extra freedom degree is introduced.

tions of the beam with two pre-existing cracks are illustrated in Fig. 18. The material properties are 3 × 105 for Young’s modulus, 0.25 for Poisson’s ratio, and the test is modeled in plane strain. A uniform mesh of 20 × 5 is used and 3 layer refinement is adopted. The crack increment is set to be 0.01. As the crack propagates, we refine the mesh in the vicinity of the new crack tip based on previous mesh. The paths of crack growth are plotted in Fig. 19, where four cyan dots represent the force bearing points. It is found that two cracks symmetrically propagate in the direction of maximal load on the opposite side of the beam. The predicated propagating trajectory is quite smooth and agrees well with results by four-node quadrilateral element NMM with continuous nodal stress (Quad4-CNS (NMM)) [32] (see Fig. 20(a)) and by virtual node method for polygonal element (VNMPE) [54] (see Fig. 20(b)) . This example shows that the proposed method is promising for crack propagation problems.

5.4. Crack growth of a beam under four-point bending This bending test, initially proposed by Bocca et al. [53], has been studied widely [32,54]. The initial configuration and boundary condi49

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51 Table 3 Normalized SIFs comparison for a star-shaped crack. a/W

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

𝐹𝐼𝐴

𝐹𝐼𝐵

#

∗

∗∗

∗∗∗

0.7443 0.7631 0.7918 0.8331 0.8905 0.9736 1.1091 1.3524 1.9080

0.7408 0.7570 0.7846 0.8255 0.8815 0.9758 1.1142 — —

0.7511 0.7670 0.7931 0.8287 0.8864 0.9673 1.0971 1.3423 1.9037

0.758 0.771 0.789 0.821 0.887 0.971 1.107 1.340 1.930

𝐵 𝐹𝐼𝐼

#

∗

∗∗

∗∗∗

0.7430 0.7630 0.7960 0.8464 0.9243 1.0439 1.2393 1.5766 2.2336

0.7408 0.757 0.7884 0.8365 0.9087 1.0448 1.1936 — —

0.7690 0.7683 0.7983 0.8466 0.9255 1.0445 1.2367 1.5624 2.1927

0.767 0.771 0.798 0.854 0.924 1.040 1.234 1.559 2.166

#

∗

∗∗

∗∗∗

0.0001 0.0004 0.0022 0.0085 0.0184 0.0372 0.0633 0.0871 0.0903

0.0000 0.0004 0.0022 0.0070 0.0168 0.0338 0.0529 — —

0.0001 0.0005 0.0021 0.0080 0.0184 0.0364 0.0593 0.0864 0.0868

0.000 0.000 0.002 0.007 0.016 0.036 0.061 0.082 0.089

(‘–’ means no corresponding solution). # Current results using MLS-based NMM; ∗ Cheung et al. [51]; ∗ ∗ Daux et al. [52]; ∗ ∗ ∗ Ma et al. [31].

Fig. 20. Crack growth paths of a beam under four-point bending by (a). four-node quadrilateral element NMM with continuous nodal stress [32], and (b). virtual node method for polygonal element [54].

Acknowledgments This work has been supported by the Natural Science Foundation of China (under #11602165 and #51479131), and Supported by Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Grant NO. Z015010. K.X.’s research was partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grant #72031326.

Fig. 19. Crack growth paths of beam under four-point bending at different steps.

6. Conclusions References In this paper, a structured mesh refinement technique has been developed in MLS-based numerical manifold method. Hierarchical meshes are used around each crack tip. An eccentric domain of influence is suggested near the transition region to ensure integral accuracy, as well as to reduce the quantity of calculation. The visibility criterion is introduced to cope with the discontinuity near the crack tip, and is proven to be simple and adequate. To demonstrate the ability of the proposed method, four examples are tested and the results are in good agreement with literatures. It is shown that 3 or 4 layer local refinement is sufficient enough around each crack tip. The proposed procedure can handle complex crack problems in a natural way without any extra degree of freedoms, and the crack tip can stop anywhere. Further investigations of proposed method are worthy and attractive, especially for problems of dynamic crack propagation and three-dimensional crack propagation.

[1] Wawrzynek PA, Ingraffea AR. An interactive approach to local remeshing around a propagating crack. Finite Elem Anal Des 1989;5(1):87–96. [2] Habraken A, Cescotto S. An automatic remeshing technique for finite element simulation of forming processes. Int J Numer Methods Eng 1990;30(8):1503–25. [3] Lo S, Lee C. Solving crack problems by an adaptive refinement procedure. Eng Fract Mech 1992;43(2):147–63. [4] Bouchard P, Bay F, Chastel Y, Tovena I. Crack propagation modelling using an advanced remeshing technique. Comput Methods Appl Mech Eng 2000;189(3):723–42. [5] Khoei AR, Azadi H, Moslemi H. Modeling of crack propagation via an automatic adaptive mesh refinement based on modified superconvergent patch recovery technique. Eng Fract Mech 2008;75(10):2921–45. [6] Belytschko T, Lu YY, Gu L. Crack propagation by element-free Galerkin methods. Eng Fract Mech 1995;51(2):295–315. [7] Belytschko T, Lu YY, Gu L. Element‐free Galerkin methods. Int J Numer Methods Eng 1994;37(2):229–56. [8] Liu GR, Tu ZH. An adaptive procedure based on background cells for meshless methods. Comput Methods Appl Mech Eng 2002;191(s17–18):1923–43. 50

F. Liu, K. Xia

Engineering Analysis with Boundary Elements 84 (2017) 42–51

[9] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. Int J Numer Methods Eng 2005;63(11):1559–82. [10] Schweitzer MA. An adaptive HP-version of the multilevel particle–partition of unity method. Comput Methods Appl Mech Eng 2009;198(13–14):1260–72. [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] Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Meth Engng 2000;48:1741–60. [13] Ventura G, Xu JX, Belytschko T. A vector level set method and new discontinuity approximations for crack growth by EFG. Int J Numer Methods Eng 2002;54(6):923–44. [14] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69(7):813–33. [15] Stazi FL, Budyn E, Chessa J, Belytschko T. An extended finite element method with higher-order elements for curved cracks. Comput Mech 2003;31(1–2):38–48. [16] Shi GH. Manifold Method of Material Analysis. In: Proceedings of the 9th army conference on applied mathematics and computing. Minneapolis, Minnesota, F: Army research office; 1991. [17] Shi GH, Goodman RE. Two dimensional discontinuous deformation analysis. Int J Numer Anal Methods Geomech 1985;9(6):541–56. [18] Wu Z, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech 2012;39:38–53. [19] Wang SL, Ge XR. Application of manifold method in simulating crack propagation. Chin J Rock Mech Eng 1997;16(5):405–10. [20] Yang YT, Tang XH, Zheng H, Liu QS, He L. Three-dimensional fracture propagation with numerical manifold method. Eng Anal Bound Elem 2016;72:65–77. [21] Tsay RJ, Chiou YJ, Chuang WL. Crack growth prediction by manifold method. J Eng Mech 1999;125(8):884–90. [22] Chiou YJ, Lee YM, Tsay RJ. Mixed mode fracture propagation by manifold method. Int J Fract 2002;114(4):327–47. [23] 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. [24] Liu ZJ, Zheng H. Two-dimensional numerical manifold method with multilayer covers. Sci China Technol Sci 2016;59(4):515–30. [25] 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. [26] 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. [27] 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. [28] 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. [29] Li SC, Cheng YM. Enriched meshless manifold method for two-dimensional crack modeling. Theor Appl Fract Mech 2005;44(3):234–48. [30] Wong LNY, Wu Z. Application of the numerical manifold method to model progressive failure in rock slopes. Eng Fract Mech 2014;119:1–20. [31] 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.

[32] Yang YT, Sun GH, Zheng H, Fu XD. A four-node quadrilateral element fitted to numerical manifold method with continuous nodal stress for crack analysis. Comput Struct 2016;177:69–82. [33] Fries TP. A corrected XFEM approximation without problems in blending elements. Int J Numer Methods Eng 2008;75(5):503–32. [34] Tarancon JE, Vercher A, Giner E, Fuenmayor FJ. Enhanced blending elements for XFEM applied to linear elastic fracture mechanics. Int J Numer Methods Eng 2009;77(1):126–48. [35] Zhu QZ. On enrichment functions in the extended finite element method. Int J Numer Methods Eng 2012;91(2):186–217. [36] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Methods Eng 1987;24(2):337–57. [37] 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. [38] Dolbow J, Belytschko T. Numerical integration of the Galerkin weak form in meshfree methods. Comput Mech 1999;23(3):219–30. [39] Ma GW, An XM, He L. The numerical manifold method: a review. Int J Comput Methods 2010;07(01):1–32. [40] Samet H. The quadtree and related hierarchical data structures. ACM Comput Surv (CSUR) 1984;16(2):187–260. [41] Yerry M, Shephard M. A modified quadtree approach to finite element mesh generation. IEEE Comput Graph Appl 1983;3(1):39–46. [42] Fries TP, Byfut A, Alizada A, Cheng KW, Schröder A. Hanging nodes and XFEM. Int J Numer Methods Eng 2011;86(4–5):404–30. [43] Tabarraei A, Sukumar N. Adaptive computations on conforming quadtree meshes. Finite Elem Anal Des 2005;41(7–8):686–702. [44] Zheng C, Wu SC, Tang XH, Zhang JH. A meshfree poly-cell Galerkin (MPG) approach for problems of elasticity and fracture. Comput Model Eng 2008:149–78. [45] Belytschko T, Gu L, Lu Y. Fracture and crack growth by element free Galerkin methods. Model Simul Mater Sci Eng 1994;2(3A):519. [46] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Eng Asme 1963;85(4):527. [47] Hellen TK. On the method of virtual crack extensions. Int J Numer Methods Eng 1975;9(1):187–207. [48] Stern M, Becker EB, Dunham RS. A contour integral computation of mixed-mode stress intensity factors. Int. J Fract 1976;12(3):359–68. [49] Timoshenko S, Goodier JN. Theory of elasticity. 3d ed. New York: McGraw-Hill; 1970. [50] Tada H, Paris PC, Irwin GR. The stress analysis of cracks handbook. 3rd ed. New-York: ASME Press; 2000. [51] Cheung YK, Woo CW, Wang YH. A general method for multiple crack problems in a finite plate. Comput Mech 1992;10(5):335–43. [52] Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Methods Eng 2000;48(12):1741–60. [53] Bocca P, Carpinteri A, Valente S. Size effects in the mixed mode crack propagation: softening and snap-back analysis. Eng Fract Mech 1990;35(1–3):159–70. [54] Tang XH, Wu SC, Zheng C, Zhang JH. A novel virtual node method for polygonal elements. Appl Math Mech 2009;30(10):1233–46.

51