# Hybrid Analytical and MLS-Based NMM for the Determination of Generalized Stress Intensity Factors

Authors: Feng Liu; Hong Zheng; Xiuli Du
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The numerical manifold method (NMM) is characterized by its two cover systems, the mathematical cover and the physical cover. In the standard NMM, the mathematical cover is required to cover the whole problem domain. In this study, however, around each crack tip we specify a small domain on which the displacement is taken as the truncated Williams’ displacement series. And accordingly all such small domains are not covered by the mathematical cover that only covers the rest of the problem domain. Meanwhile, the mathematical cover is constructed by designating all supports of the scattered nodes arising in the moving least squares interpolation as the mathematical patches. In this way, any physical patch contains no crack tip and can be approximated by polynomials. As a result, no blending element issue exists as in the extended finite element method and NMM. In addition to high precision, the proposed procedure is especially suitable for the situation where a crack tip is very close to other cracks, a case difficult to treat by the interaction integral procedure that is commonly used in the extraction of the stress intensity factors of mixed mode cracks.


---

## Body

## 1. Introduction

The extended finite element method (XFEM [1, 2]) was first introduced by Belytschko and Black and Moës et al. in 1999. The main idea of XFEM is to superpose enriched terms to the classical finite element approximation. The crack is represented independently of the mesh, allowing for simulation of the crack propagation without remeshing. It is known that the convergence rate of XFEM is influenced by the transition layer between the singular enrichment area and the rest of the domain. Some strategies for improving the convergence rate of XFEM can be found in [3–5]. Another effective approach to model crack problems is the numerical manifold method (NMM [6–14]). Besides XFEM and NMM, there are some other related methods, such as mesh-free method [15–18], cracking particles method [19], extended isogeometric analysis (XIGA [20]), and immersed particles method [21].Usually, neither XFEM nor NMM directly estimates the stress intensity factors (SIFs). Instead, they utilize a postprocessing procedure, called the interaction integral procedure [2], to extract SIFs of a mixed mode crack. If the crack is very close to another crack, however, the interaction integral procedure cannot work very well because a local domain containing the crack tip, selected for the calculation of the J integral, might be cut by the neighboring crack unless the local domain is selected very small. But if the local domain is too small, the precision of the resulting SIFs is compromised considerably.To solve the above issue, some special techniques have been developed to conduct the crack analysis, where the SIFs are selected as the basic unknowns and can be solved directly. Liu et al. [22] enriched the approximation surrounding the crack tip with the truncated Williams’ series and force the enriched DOFs to be equal. The enriched DOFs are actually the generalized SIFs. The numerical experiences indicate that this approach does not always reach high precision since the blending elements are still present.Réthoré et al. [23] developed a hybrid analytical and extended finite element method, abbreviated as HAX-FEM. By HAX-FEM, the problem domain Ω is subdivided into two overlapping subdomains ΩW and ΩX. Here, ΩW is a small domain around the crack tip and, for the simplicity in presentation, Ω is assumed to contain only one crack tip. On ΩW the displacement is expressed by a truncated Williams’ series, while the displacement on ΩX is approximated by the conventional XFEM. Due to the presence of the overlapping subdomain, the energy coupling procedure was utilized in the weak formulation so as not to commit the variational crime. In order not to avoid the overlapping subdomain, another artificial thing, Passieux et al. [24] proposed that ΩW and ΩX do not overlap, and the displacement continuity across the interface between ΩW and ΩX is enforced by the Lagrange multiplier method, which will lead to an indefinite system of linear equations. And the approximation to the Lagrange multiplier is not so easy such that the LBB condition is satisfied. To avoid the introduction of the Lagrange multiplier, recently Zhuang et al. [25] suggested that in ΩW the stress is taken as the basic field variable instead of the displacement, where the weak formulation is derived from the mixed energy functional proposed by Long and Zhao [26]. This procedure is able to achieve an accurate evaluation of SIFs, but the displacement in ΩW is hard to calculate accurately. The open degree of cracks is important for the situations such as the exploitation of shale gas by the hydraulic fracture technique.In this study, the MLS-based NMM is used to evaluate directly the SIFs of mixed mode cracks. The MLS-based NMM enjoys the merits owned by the element-free Galerkin (EFG [15]) method and the NMM [6], which has been verified in the analysis of free boundary value problems such as the unconfined seepage analysis [27], and will be further confirmed by the challenging examples on crack analysis in the sequel.


## 2. Elastic Crack Tip Asymptotic Fields

The crack tip asymptotic fields, including displacement and stress fields, of a plane crack in a linear elastic brittle material with traction-free faces can be expressed in the so-called Williams’ expansions [28]. If the crack lies on the negative x-axis and the polar coordinates centered at the crack tip are designated r and θ (counterclockwise from the positive x-axis), Williams’ series (truncated to N terms) for the displacement can be written as(1)ux=∑n=0Nrn/22μAnκ+n2+-1ncos⁡n2θhhhhhhhhhhhhh-n2cos⁡n2-2θhhhhhhhhh-Bnκ+n2--1nsinn2θhhhhhhhhhhhhhhh-n2sinn2-2θ,uy=∑n=0Nrn/22μAnκ-n2--1nsinn2θhhhhhhhhhhhhh+n2sinn2-2θhhhhlhhhh+Bnκ-n2+-1ncos⁡n2θhhhhhhhhhhhhhh+n2cos⁡n2-2θ,where μ=E/(2(1+v)) is the shear modulus; the Kolosov constant κ is equal to (3-v)/(1+v) for plane stress and 3-4v for plane strain conditions; E and v are Young’s modulus and Poisson’s ratio. The terms containing An and Bn are related to the mode I and mode II parts of deformation, respectively. Here is a brief review of some notable properties of this expansion. The rigid body translation of the crack depends on the coefficients A0 and B0 and the rigid body rotation with respect to the crack tip depends on B2. The second term in the mode I expansion is often referred to as the elastic T-stress component. The coefficients A1 and B1 are related to the mode I and mode II stress intensity factors (SIFs) as(2)A1=KI2π,B1=-KII2π.Traditionally, the SIFs are important for the determination of the initiation and propagation of cracks in brittle materials. However, some recent studies show that higher-order terms of the asymptotic field are of great relevance for predicting the constraint of elastoplastic crack tip fields [29] and for interpreting the size effect of quasi-brittle materials. So evaluating the higher-order terms of Williams’ expansions, called generalized stress intensity factors together with the conventional SIFs, is also necessary.


## 3. Problem Formulation

Let us consider the elastic body with a traction-free crack surfaceΓc shown in Figure 1 (reproduced from [24]). The prescribed displacement u- is imposed on Γu, while traction t0 is acted on Γt. The domain is divided into two nonoverlapping subdomains ΩX and ΩW with unknown displacements uX and uW, respectively. In domain ΩX, the approximation of displacement uX is described by the MLS-based NMM, which will be discussed in the next section. In domain ΩW, we approximate the displacement by the truncated Williams’ series given by (1).Figure 1
Decomposition of the domain into a standard domain and an analytical patch.In [24], the displacement continuity across the interface ΓW is enforced by the Lagrange multiplier method. As the number of crack tips increases, the number of Lagrange multipliers will increase rapidly. So the penalty method is adopted in this study. In fact, the results created by the two methods are very small due to the numerical excellence of NMM; see [30], for example. In this way, we construct the conventional potential as follows:(3)πuX,uW=∫ΩX12εuXTDε(uX)dΩ+∫ΩW12εuWTDε(uW)dΩ+∫ΓW12kpuX-uW2dΓ-∫ΓtuXTt0dΓ+∫ΓuλTuX-u-dΓ,where ε is the infinitesimal strain tensor, D is elastic matrix, ΓW is the interface between ΩX and ΩW, kp is the penalty factor, and λ is the Lagrange multiplier associated with the essential boundary condition. We assume that the boundary conditions are held by ΩX only since ΩW is a quite small domain around the crack tip. The equation can be extended to multiple crack tips problem straightly.


### 3.1. Standard Subdomain

In domainΩX, MLS-based NMM is used to approximate the displacement. The salient feature of NMM is the introduction of two covers, namely, the mathematical cover (MC) and the physical cover (PC). The MC consists of mathematical patches (MP), each of which is a rectangle or hexagon usually. The MCs do not need to match the boundary of the domain, as long as they cover the whole domain. The construction of physical patches is accomplished by cutting the mathematical patches using discontinuity lines, such as cracks and material boundaries. All the physical patches form the PC covering exactly the problem domain. The main purpose to introduce two covers is to solve in a unified way the continuum and discontinuum problems. Moreover, the best approximation to the primal field variables can be guaranteed because we can always use a uniform mesh to construct the MC.Here, we give a brief illustration on how NMM simulates the discontinuity. First, we deploy overΩX some nodes. Unlike EFG, all these nodes, called mathematical nodes, can be allowed to be outside ΩX, as long as their influence domains have intersection with ΩX. All the influence domains of these mathematical nodes constitute the MC covering ΩX. Let us cut the mathematical patch of node P (denoted as quadrilateral ABCD) by a crack (the bold line) into 2 patches: pentagon AEFCD and triangle EBF, respectively, as shown in Figure 2. Both of the two patches are physical patches of the PC in construction. We add a new physical node, denoted by P2, and associate P2 with triangle EBF. At this moment, P2 has the same coordinate as P. The degrees of freedom associated with the two physical nodes P and P2 are separate, reproducing strong discontinuity along the crack. So, in NMM, there is no need to introduce the Heaviside enrichment function which is widely used in XFEM. Details can be found in [9, 12].Figure 2
A crack cuts one mathematical patch into two physical patches.For zero-order NMM, which is used in this paper, the displacement interpolation can be expressed as(4)uX(x)=∑i∈NPNi(x)ui(x)=NXUX,where NP is the set of all physical nodes, ui is the degree of freedom, and Ni(x) is the corresponding shape function with both associated with node-i. Here, we use the MLS shape function with linear basis and rectangular influence domain instead of the FEM shape function. The MLS interpolation is discussed in detail in [15, 31].


### 3.2. Analytical Patch

On domainΩW, the displacement is described by the truncated Williams’ series (1) only. So no mesh is needed theoretically. Rewrite (1) as(5)uW(x)=uxuy=∑n=0Nrn/22μf1f2f3f4AnBn=NWUW,where N is the maximum order of the truncated Williams’ series; f1, f2, f3, and f4 are extracted from (1). Since the SIFs and higher-order terms are related to special An and Bn above, generalized stress intensity factors can be evaluated directly.


### 3.3. Matching Interface

The third term on the right of (3) stands for the matching of standard subdomain and analytical patch. No more unknowns are introduced through penalty method.


### 3.4. System of Equations

Though no mesh is needed for MLS interpolation, a background mesh is essential for numerical integration. In practice, a circle of radiusRW=hrW centered at the crack tip is defined first. Here, h is the average element size of the mesh and rW is a dimensionless size controlling the region of ΩW. All the elements in or in intersection with the circle constitute subdomain ΩW, while the rest constitute subdomain ΩX. The outer boundary of ΩW, which equals the inner boundary of ΩX, is the common interface ΓW.PuttingδπuX,uW=0 and replacing uX(x) and uW(x) by (4) and (5), we have the equations system as follows:(6)KX+HXLTGTLKW+HW0G00uaΛ=P0Q,where(7)KX=∫ΩXBXTDBXdΩ,KW=∫ΩWBWTDBWdΩ,HX=∫ΓXNXTNXdΩ,HW=∫ΓWNWTNWdΩ,L=∫ΓXNWTNXdΩ,G=∫ΓuNλTNXdΩ,P=∫ΓtNXTt0dΩ,Q=∫ΓuNλTu¯dΩ.u is the degrees of freedom in domain ΩX, a is the degrees of freedom in domain ΩW, and Λ is the Lagrange multipliers associated with Γu.To obtainBX and BW, we should utilize the compatibility equation:(8)εij=12∂ui∂xj+∂uj∂xi.Substituting (4) into (8) and comparing the result with equation εX(x)=BXUX, matrix BX is obtained. We can derive BW in a similar way.The background mesh is used for integrating the matrices in (6). For the elements cut by no discontinuity, 4×4 Gauss’s points are used. For the elements cut by cracks or boundaries they are subdivided into a set of triangles and each triangle integrated by 13-point Hammer integration. For the element with a crack tip inside, a special technique [12] is used to eliminate the singularity.


## 3.1. Standard Subdomain

In domainΩX, MLS-based NMM is used to approximate the displacement. The salient feature of NMM is the introduction of two covers, namely, the mathematical cover (MC) and the physical cover (PC). The MC consists of mathematical patches (MP), each of which is a rectangle or hexagon usually. The MCs do not need to match the boundary of the domain, as long as they cover the whole domain. The construction of physical patches is accomplished by cutting the mathematical patches using discontinuity lines, such as cracks and material boundaries. All the physical patches form the PC covering exactly the problem domain. The main purpose to introduce two covers is to solve in a unified way the continuum and discontinuum problems. Moreover, the best approximation to the primal field variables can be guaranteed because we can always use a uniform mesh to construct the MC.Here, we give a brief illustration on how NMM simulates the discontinuity. First, we deploy overΩX some nodes. Unlike EFG, all these nodes, called mathematical nodes, can be allowed to be outside ΩX, as long as their influence domains have intersection with ΩX. All the influence domains of these mathematical nodes constitute the MC covering ΩX. Let us cut the mathematical patch of node P (denoted as quadrilateral ABCD) by a crack (the bold line) into 2 patches: pentagon AEFCD and triangle EBF, respectively, as shown in Figure 2. Both of the two patches are physical patches of the PC in construction. We add a new physical node, denoted by P2, and associate P2 with triangle EBF. At this moment, P2 has the same coordinate as P. The degrees of freedom associated with the two physical nodes P and P2 are separate, reproducing strong discontinuity along the crack. So, in NMM, there is no need to introduce the Heaviside enrichment function which is widely used in XFEM. Details can be found in [9, 12].Figure 2
A crack cuts one mathematical patch into two physical patches.For zero-order NMM, which is used in this paper, the displacement interpolation can be expressed as(4)uX(x)=∑i∈NPNi(x)ui(x)=NXUX,where NP is the set of all physical nodes, ui is the degree of freedom, and Ni(x) is the corresponding shape function with both associated with node-i. Here, we use the MLS shape function with linear basis and rectangular influence domain instead of the FEM shape function. The MLS interpolation is discussed in detail in [15, 31].


## 3.2. Analytical Patch

On domainΩW, the displacement is described by the truncated Williams’ series (1) only. So no mesh is needed theoretically. Rewrite (1) as(5)uW(x)=uxuy=∑n=0Nrn/22μf1f2f3f4AnBn=NWUW,where N is the maximum order of the truncated Williams’ series; f1, f2, f3, and f4 are extracted from (1). Since the SIFs and higher-order terms are related to special An and Bn above, generalized stress intensity factors can be evaluated directly.


## 3.3. Matching Interface

The third term on the right of (3) stands for the matching of standard subdomain and analytical patch. No more unknowns are introduced through penalty method.


## 3.4. System of Equations

Though no mesh is needed for MLS interpolation, a background mesh is essential for numerical integration. In practice, a circle of radiusRW=hrW centered at the crack tip is defined first. Here, h is the average element size of the mesh and rW is a dimensionless size controlling the region of ΩW. All the elements in or in intersection with the circle constitute subdomain ΩW, while the rest constitute subdomain ΩX. The outer boundary of ΩW, which equals the inner boundary of ΩX, is the common interface ΓW.PuttingδπuX,uW=0 and replacing uX(x) and uW(x) by (4) and (5), we have the equations system as follows:(6)KX+HXLTGTLKW+HW0G00uaΛ=P0Q,where(7)KX=∫ΩXBXTDBXdΩ,KW=∫ΩWBWTDBWdΩ,HX=∫ΓXNXTNXdΩ,HW=∫ΓWNWTNWdΩ,L=∫ΓXNWTNXdΩ,G=∫ΓuNλTNXdΩ,P=∫ΓtNXTt0dΩ,Q=∫ΓuNλTu¯dΩ.u is the degrees of freedom in domain ΩX, a is the degrees of freedom in domain ΩW, and Λ is the Lagrange multipliers associated with Γu.To obtainBX and BW, we should utilize the compatibility equation:(8)εij=12∂ui∂xj+∂uj∂xi.Substituting (4) into (8) and comparing the result with equation εX(x)=BXUX, matrix BX is obtained. We can derive BW in a similar way.The background mesh is used for integrating the matrices in (6). For the elements cut by no discontinuity, 4×4 Gauss’s points are used. For the elements cut by cracks or boundaries they are subdivided into a set of triangles and each triangle integrated by 13-point Hammer integration. For the element with a crack tip inside, a special technique [12] is used to eliminate the singularity.


## 4. Numerical Examples

Several examples are studied to show the efficiency of the proposed method. The first one is an edge-cracked plate under remote tension. The choice of truncated order and the convergence are investigated. Secondly the centre-cracked plate under uniaxial tension (CCPT) and pure shearing (CCPS) are used to evaluate the generalized stress intensity factors. At last a star crack is provided to show the excellent performance of proposed method for complex cracks. In all these examples, the rectangular weight function with double node spacing is adopted.


### 4.1. Edge-Cracked Plate under Remote Tension

The edge-cracked plate considered here is shown in Figure3(a). To compare with [24], the parameters are as follows: L=17 mm, W=7 mm, a=3.5 mm, Poisson’s ratio v=0.3, and Young’s modulus E=200 GPa. A state of plane stress condition is assumed.Figure 3
Edge-cracked plate under remote tension: (a) geometry and loading; (b) mesh and patch forrW=2.
(a)
(b)A uniform mesh of49×119 elements is used. The mesh, the patch (with rW=2), and the crack (the red line) are plotted in Figure 3(b). The reference solution [24] of KI for this problem is KI0=2.9637 MPam.The normalized value ofKI is plotted as a function of the dimensionless size rW in Figure 4 with different orders N of truncated Williams’ expansion. A similar tendency is observed as in [24]. The dependence on the dimensionless size rW vanishes when N is higher than 2. For the same N, the precision decreases as the size rW increases, indicating that the analytical patch only influences near the crack tip. Since the influence domain of nodes in MLS is larger than that of XFEM, rW=1 leads to a wrong solution. For general case, rW=2 or rW=3 is recommended and used in the following examples.Figure 4
NormalizedKI as function of the size of the patch for different values of N.To consider the convergence of proposed method, the relative error ofKI is plotted as a function of the mesh size h in Figure 5 with different orders N of truncated Williams’ series. The relative error of KI is defined as(9)Kerr⁡=KI-KI0KI0.The convergence rates of relative error of KI are around 2, larger than the value 1.25 to 1.5 in [24]. By the convergence rate we mean the average slope of log⁡(Kerr⁡)/log⁡(h). Again we can see that when N is above 2, accurate SIFs can be obtained.Figure 5
Relative error ofKI as function of node spacing for different values of N.For the range of penalty factor in 103~105E, almost the same results are generated.


### 4.2. A Centre-Cracked Plate under Uniaxial Tension or Pure Shearing

The centre-cracked plate under uniaxial tension (CCPT) is selected to investigate the pure mode I crack parameters (Figure6(a)). In order to compare the results with the reference, the following geometry and loading parameters are used: W=h=4, σ=1, and crack length a=2. A 80×80 mesh is used. In this case, N=14 together with rW=3 are used to calculate the coefficients of Williams’ expansion. The penalty factor is assumed to be 1000E. Figure 6(b) shows the σy-stress nephogram, where the red lines represent the common interface of ΩX and ΩW. Table 1 shows the coefficients A1–A5, B1,and B3–B5 for CCPT specimens in comparison with the results provided by other researchers [32, 33]. It seems that the results here are in good agreement with the results of finite element overdeterministic (FEOD [21]) method. Only A1–A5 are given for hybrid crack element (HCE) and boundary collocation method (BCM) in [32]. In fact, the results of HCE in Table 1 are calculated by half-crack model with only mode I considered. HCE considering half-crack model with mixed mode leads to wrong results for the coefficients of higher-order parameters.Table 1
CoefficientsA1–A5,  B1, and B3–B5 for CCPT specimens (a/W=0.25).
A
1
B
1
A
2
A
3
B
3
A
4
B
4
A
5
B
5
Present
0.7674
−0.0002
−0.2761
0.1830
0.0001
0.0026
0.0001
−0.0235
−0.0003
FEOD
0.7679
0.0000
−0.2778
0.1866
0.0000
0.0030
−0.0001
−0.0275
−0.0000
HCE
0.7665
—
−0.2779
0.1915
—
−0.0018
—
−0.0235
—
BCM
0.7680
—
−0.2777
0.1866
—
0.0030
—
−0.0279
—Figure 6
The CCPT specimen: (a) geometry and loading; (b)σy-stress nephogram for rW=3.
(a)
(b)Now consider the centre-cracked plate under pure shearing (CCPS). It is a typical example of mode II. The geometrical parameters of the CCPS (Figure7(a)), the mesh, and material properties are the same as CCPT. Here, a uniformly distributed shear force τ=1 is imposed. Figure 7(b) gives the τxy-stress nephogram and again the red lines represent the interfaces. The coefficients A1–A5, B1, and B3–B5 are shown in Table 2. Good agreement is observed between the present results and those of FEOD, HCE, and BCM. Again results of HCE in Table 2 are calculated by half-crack model with only mode II considered.Table 2
CoefficientsA1–A5, B1, and B3–B5 for CCPS specimens (a/W=0.25).
A
1
B
1
A
2
A
3
B
3
A
4
B
4
A
5
B
5
Present
−0.0002
−0.7422
0.0000
0.0005
−0.1848
−0.0009
−0.0047
−0.0002
0.0264
FEOD
0.0010
−0.7403
−0.0011
0.0011
−0.1840
−0.0003
−0.0024
−0.0000
0.0234
HCE
—
−0.7414
—
—
−0.1835
—
−0.0045
—
0.0270
BCM
—
−0.7422
—
—
−0.1841
—
−0.0034
—
0.0242Figure 7
The CCPS specimen: (a) geometry and loading; (b)τxy-stress nephogram for rW=3.
(a)
(b)Both of these two examples illustrate that the present method can predict the generalized stress intensity factors precisely. No special treatment is needed throughout the whole process.


### 4.3. A Star-Shaped Crack in a Square Plate

At last a star-shaped crack is analyzed which is often viewed as a representation of complex cracks. The star-shaped crack in a square plate subjected to biaxial tension is shown in Figure8. The normalized stress intensity factors at tips A and B are defined as(10)FIA=KIAσπa,FIB=KIBσπa,FIIB=KIIBσπa.Figure 8
A star-shaped crack in a square plate under biaxial tension.In this example, the plate dimension is fixed to beW=5, and material constants are Young’s modulus E=200000 and Poisson’s ratio v=0.3. The biaxial tension σ is taken to be unity. The influence of the finiteness of plate on the SIFs is also investigated. Since the mesh refinement with no limit is one of the major advantages for MLS interpolation, hierarchical meshes are used around every crack tip here. The mesh (with 4 layers of hierarchical meshes near the tip) and detailed view with analytical patches are plotted in Figure 9 with a/W=0.1. The cracks are moved with a small distance to avoid special case caused by symmetry. Here N=5 and rW=3 are selected. The computed SIF results are summarized and compared with three reference solutions, respectively, from Cheung et al. [34], Daux et al. [35], and Ma et al. [7] in Table 3. From this table, good agreement can be easily seen. This example also shows that the proposed method can implement local refinement and handle complex crack straightly.Table 3
Normalized SIFs comparison for a star-shaped crack.
a
/
W
F
I
A
F
I
B
F
II
B
T
*
**
***
T
*
**
***
T
*
**
***
0.1
0.7447
0.7408
0.7511
0.758
0.7463
0.7408
0.7690
0.767
−0.0003
0.0000
0.0001
0.000
0.2
0.7629
0.7570
0.7670
0.771
0.7649
0.757
0.7683
0.771
0.0005
0.0004
0.0005
0.000
0.3
0.7913
0.7846
0.7931
0.789
0.7961
0.7884
0.7983
0.798
0.0018
0.0022
0.0021
0.002
0.4
0.8319
0.8255
0.8287
0.821
0.8470
0.8365
0.8466
0.854
0.0074
0.0070
0.0080
0.007
0.5
0.8887
0.8815
0.8864
0.887
0.9229
0.9087
0.9255
0.924
0.0186
0.0168
0.0184
0.016
0.6
0.9712
0.9758
0.9673
0.971
1.0404
1.0182
1.0445
1.040
0.0355
0.0338
0.0364
0.036
0.7
1.1035
1.1142
1.0971
1.107
1.2334
1.1936
1.2367
1.234
0.0600
0.0529
0.0593
0.061
0.8
1.3398
—
1.3423
1.340
1.5631
—
1.5624
1.559
0.0847
—
0.0864
0.082
0.9
1.8801
—
1.9037
1.930
2.1968
—
2.1927
2.166
0.0870
—
0.0868
0.089
“—” means no corresponding solution.
*Cheung et al. (1992) [34];  **Daux et al. (2000) [35]; ***Ma et al. (2009) [7].Figure 9
Star-shaped crack and the hierarchical meshes: (a) global view; (b) local view with analytical patches.
(a)
(b)


## 4.1. Edge-Cracked Plate under Remote Tension

The edge-cracked plate considered here is shown in Figure3(a). To compare with [24], the parameters are as follows: L=17 mm, W=7 mm, a=3.5 mm, Poisson’s ratio v=0.3, and Young’s modulus E=200 GPa. A state of plane stress condition is assumed.Figure 3
Edge-cracked plate under remote tension: (a) geometry and loading; (b) mesh and patch forrW=2.
(a)
(b)A uniform mesh of49×119 elements is used. The mesh, the patch (with rW=2), and the crack (the red line) are plotted in Figure 3(b). The reference solution [24] of KI for this problem is KI0=2.9637 MPam.The normalized value ofKI is plotted as a function of the dimensionless size rW in Figure 4 with different orders N of truncated Williams’ expansion. A similar tendency is observed as in [24]. The dependence on the dimensionless size rW vanishes when N is higher than 2. For the same N, the precision decreases as the size rW increases, indicating that the analytical patch only influences near the crack tip. Since the influence domain of nodes in MLS is larger than that of XFEM, rW=1 leads to a wrong solution. For general case, rW=2 or rW=3 is recommended and used in the following examples.Figure 4
NormalizedKI as function of the size of the patch for different values of N.To consider the convergence of proposed method, the relative error ofKI is plotted as a function of the mesh size h in Figure 5 with different orders N of truncated Williams’ series. The relative error of KI is defined as(9)Kerr⁡=KI-KI0KI0.The convergence rates of relative error of KI are around 2, larger than the value 1.25 to 1.5 in [24]. By the convergence rate we mean the average slope of log⁡(Kerr⁡)/log⁡(h). Again we can see that when N is above 2, accurate SIFs can be obtained.Figure 5
Relative error ofKI as function of node spacing for different values of N.For the range of penalty factor in 103~105E, almost the same results are generated.


## 4.2. A Centre-Cracked Plate under Uniaxial Tension or Pure Shearing

The centre-cracked plate under uniaxial tension (CCPT) is selected to investigate the pure mode I crack parameters (Figure6(a)). In order to compare the results with the reference, the following geometry and loading parameters are used: W=h=4, σ=1, and crack length a=2. A 80×80 mesh is used. In this case, N=14 together with rW=3 are used to calculate the coefficients of Williams’ expansion. The penalty factor is assumed to be 1000E. Figure 6(b) shows the σy-stress nephogram, where the red lines represent the common interface of ΩX and ΩW. Table 1 shows the coefficients A1–A5, B1,and B3–B5 for CCPT specimens in comparison with the results provided by other researchers [32, 33]. It seems that the results here are in good agreement with the results of finite element overdeterministic (FEOD [21]) method. Only A1–A5 are given for hybrid crack element (HCE) and boundary collocation method (BCM) in [32]. In fact, the results of HCE in Table 1 are calculated by half-crack model with only mode I considered. HCE considering half-crack model with mixed mode leads to wrong results for the coefficients of higher-order parameters.Table 1
CoefficientsA1–A5,  B1, and B3–B5 for CCPT specimens (a/W=0.25).
A
1
B
1
A
2
A
3
B
3
A
4
B
4
A
5
B
5
Present
0.7674
−0.0002
−0.2761
0.1830
0.0001
0.0026
0.0001
−0.0235
−0.0003
FEOD
0.7679
0.0000
−0.2778
0.1866
0.0000
0.0030
−0.0001
−0.0275
−0.0000
HCE
0.7665
—
−0.2779
0.1915
—
−0.0018
—
−0.0235
—
BCM
0.7680
—
−0.2777
0.1866
—
0.0030
—
−0.0279
—Figure 6
The CCPT specimen: (a) geometry and loading; (b)σy-stress nephogram for rW=3.
(a)
(b)Now consider the centre-cracked plate under pure shearing (CCPS). It is a typical example of mode II. The geometrical parameters of the CCPS (Figure7(a)), the mesh, and material properties are the same as CCPT. Here, a uniformly distributed shear force τ=1 is imposed. Figure 7(b) gives the τxy-stress nephogram and again the red lines represent the interfaces. The coefficients A1–A5, B1, and B3–B5 are shown in Table 2. Good agreement is observed between the present results and those of FEOD, HCE, and BCM. Again results of HCE in Table 2 are calculated by half-crack model with only mode II considered.Table 2
CoefficientsA1–A5, B1, and B3–B5 for CCPS specimens (a/W=0.25).
A
1
B
1
A
2
A
3
B
3
A
4
B
4
A
5
B
5
Present
−0.0002
−0.7422
0.0000
0.0005
−0.1848
−0.0009
−0.0047
−0.0002
0.0264
FEOD
0.0010
−0.7403
−0.0011
0.0011
−0.1840
−0.0003
−0.0024
−0.0000
0.0234
HCE
—
−0.7414
—
—
−0.1835
—
−0.0045
—
0.0270
BCM
—
−0.7422
—
—
−0.1841
—
−0.0034
—
0.0242Figure 7
The CCPS specimen: (a) geometry and loading; (b)τxy-stress nephogram for rW=3.
(a)
(b)Both of these two examples illustrate that the present method can predict the generalized stress intensity factors precisely. No special treatment is needed throughout the whole process.


## 4.3. A Star-Shaped Crack in a Square Plate

At last a star-shaped crack is analyzed which is often viewed as a representation of complex cracks. The star-shaped crack in a square plate subjected to biaxial tension is shown in Figure8. The normalized stress intensity factors at tips A and B are defined as(10)FIA=KIAσπa,FIB=KIBσπa,FIIB=KIIBσπa.Figure 8
A star-shaped crack in a square plate under biaxial tension.In this example, the plate dimension is fixed to beW=5, and material constants are Young’s modulus E=200000 and Poisson’s ratio v=0.3. The biaxial tension σ is taken to be unity. The influence of the finiteness of plate on the SIFs is also investigated. Since the mesh refinement with no limit is one of the major advantages for MLS interpolation, hierarchical meshes are used around every crack tip here. The mesh (with 4 layers of hierarchical meshes near the tip) and detailed view with analytical patches are plotted in Figure 9 with a/W=0.1. The cracks are moved with a small distance to avoid special case caused by symmetry. Here N=5 and rW=3 are selected. The computed SIF results are summarized and compared with three reference solutions, respectively, from Cheung et al. [34], Daux et al. [35], and Ma et al. [7] in Table 3. From this table, good agreement can be easily seen. This example also shows that the proposed method can implement local refinement and handle complex crack straightly.Table 3
Normalized SIFs comparison for a star-shaped crack.
a
/
W
F
I
A
F
I
B
F
II
B
T
*
**
***
T
*
**
***
T
*
**
***
0.1
0.7447
0.7408
0.7511
0.758
0.7463
0.7408
0.7690
0.767
−0.0003
0.0000
0.0001
0.000
0.2
0.7629
0.7570
0.7670
0.771
0.7649
0.757
0.7683
0.771
0.0005
0.0004
0.0005
0.000
0.3
0.7913
0.7846
0.7931
0.789
0.7961
0.7884
0.7983
0.798
0.0018
0.0022
0.0021
0.002
0.4
0.8319
0.8255
0.8287
0.821
0.8470
0.8365
0.8466
0.854
0.0074
0.0070
0.0080
0.007
0.5
0.8887
0.8815
0.8864
0.887
0.9229
0.9087
0.9255
0.924
0.0186
0.0168
0.0184
0.016
0.6
0.9712
0.9758
0.9673
0.971
1.0404
1.0182
1.0445
1.040
0.0355
0.0338
0.0364
0.036
0.7
1.1035
1.1142
1.0971
1.107
1.2334
1.1936
1.2367
1.234
0.0600
0.0529
0.0593
0.061
0.8
1.3398
—
1.3423
1.340
1.5631
—
1.5624
1.559
0.0847
—
0.0864
0.082
0.9
1.8801
—
1.9037
1.930
2.1968
—
2.1927
2.166
0.0870
—
0.0868
0.089
“—” means no corresponding solution.
*Cheung et al. (1992) [34];  **Daux et al. (2000) [35]; ***Ma et al. (2009) [7].Figure 9
Star-shaped crack and the hierarchical meshes: (a) global view; (b) local view with analytical patches.
(a)
(b)


## 5. Conclusions

A hybrid analytical and numerical manifold method is proposed to analyze crack problems. The analytical subdomain near the crack tip is described by the truncated Williams’ expansions and the rest is modeled by the MLS-based NMM. Then, the two parts are coupled by the penalty method. The variational formulation and discrete equations are carefully discussed. The introduction of the analytical patch is to directly determine generalized stress intensity factors, while the advantage of MLS-based NMM is the simplicity in dealing with complex cracks and local refinement. The numerical experiments show that the proposed method is a promising approach for crack problems.



---

*Source: 219657-2015-05-19.xml*