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

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.


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][4][5].Another effective approach to model crack problems is the numerical manifold method (NMM [6][7][8][9][10][11][12][13][14]).Besides XFEM and NMM, there are some other related methods, such as mesh-free method [15][16][17][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  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 Ω  and Ω  .Here, Ω  is a small domain around the crack tip and, for the simplicity in presentation, Ω is assumed to contain only one crack tip.On Ω  the displacement is expressed by a truncated Williams' series, while the displacement on Ω  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 Ω  and Ω  do not overlap, and the displacement continuity across the interface between Ω  and Ω  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 Ω  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 Ω  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.

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 -axis and the polar coordinates centered at the crack tip are designated  and  (counterclockwise from the positive -axis), Williams' series (truncated to  terms) for the displacement can be written as where  = /(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;  and V are Young's modulus and Poisson's ratio.The terms containing   and   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  0 and  0 and the rigid body rotation with respect to the crack tip depends on  2 .The second term in the mode I expansion is often referred to as the elastic stress component.The coefficients  1 and  1 are related to the mode I and mode II stress intensity factors (SIFs) as Traditionally, the SIFs are important for the determination of the initiation and propagation of cracks in brittle materials.However, some recent studies show that higherorder 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.

Problem Formulation
Let us consider the elastic body with a traction-free crack surface Γ  shown in Figure 1 (reproduced from [24]).The prescribed displacement u is imposed on Γ  , while traction  0 is acted on Γ  .The domain is divided into two nonoverlapping subdomains Ω  and Ω  with unknown displacements u  and u  , respectively.In domain Ω  , the approximation of displacement u  is described by the MLS-based NMM, which will be discussed in the next section.In domain Ω  , we approximate the displacement by the truncated Williams' series given by (1).
In [24], the displacement continuity across the interface Γ  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: where  is the infinitesimal strain tensor, D is elastic matrix, Γ  is the interface between Ω  and Ω  ,   is the penalty factor, and  is the Lagrange multiplier associated with the essential boundary condition.We assume that the boundary conditions are held by Ω  only since Ω  is a quite small domain around the crack tip.The equation can be extended to multiple crack tips problem straightly.

Standard Subdomain.
In domain Ω  , 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 Ω  some nodes.Unlike EFG, all these nodes, called mathematical nodes, can be allowed to be outside Ω  , as long as their influence domains have intersection with Ω  .All the influence domains of these mathematical nodes constitute the MC covering Ω  .Let us cut the mathematical patch of node  (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  2 , and associate  2 with triangle EBF.At this moment,  2 has the same coordinate as .The degrees of freedom associated with the two physical nodes  and  2 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].
For zero-order NMM, which is used in this paper, the displacement interpolation can be expressed as where   is the set of all physical nodes,   is the degree of freedom, and   () is the corresponding shape function with both associated with node-.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].

Analytical Patch.
On domain Ω  , the displacement is described by the truncated Williams' series (1) only.So no mesh is needed theoretically.Rewrite (1) as where  is the maximum order of the truncated Williams' series;  1 ,  2 ,  3 , and  4 are extracted from (1).Since the SIFs and higher-order terms are related to special   and   above, generalized stress intensity factors can be evaluated directly.

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.

System of Equations.
Though no mesh is needed for MLS interpolation, a background mesh is essential for numerical integration.In practice, a circle of radius   = ℎ  centered at the crack tip is defined first.Here, ℎ is the average element size of the mesh and   is a dimensionless size controlling the region of Ω  .All the elements in or in intersection with the circle constitute subdomain Ω  , while the rest constitute subdomain Ω  .The outer boundary of Ω  , which equals the inner boundary of Ω  , is the common interface Γ  .
Putting (u  , u  ) = 0 and replacing u  () and u  () by ( 4) and ( 5), we have the equations system as follows: where u is the degrees of freedom in domain Ω  , a is the degrees of freedom in domain Ω  , and Λ is the Lagrange multipliers associated with Γ  .
To obtain B  and B  , we should utilize the compatibility equation: Substituting ( 4) into (8) and comparing the result with equation   () = B  U  , matrix B  is obtained.We can derive B  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.

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.

Edge-Cracked Plate under Remote Tension.
The edgecracked plate considered here is shown in Figure 3(a).To compare with [24], the parameters are as follows:  = 17 mm,  = 7 mm,  = 3.5 mm, Poisson's ratio V = 0.3, and Young's modulus  = 200 GPa.A state of plane stress condition is assumed.
A uniform mesh of 49 × 119 elements is used.The mesh, the patch (with   = 2), and the crack (the red line) are plotted in Figure 3(b).The reference solution [24] of  I for this problem is  I0 = 2.9637 MPa√.
The normalized value of  I is plotted as a function of the dimensionless size   in Figure 4 with different orders  of truncated Williams' expansion.A similar tendency is observed as in [24].The dependence on the dimensionless size   vanishes when  is higher than 2. For the same , the precision decreases as the size   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,   = 1 leads to a wrong solution.For general case,   = 2 or   = 3 is recommended and used in the following examples.To consider the convergence of proposed method, the relative error of  I is plotted as a function of the mesh size ℎ in Figure 5 with different orders  of truncated Williams' series.The relative error of  I is defined as The convergence rates of relative error of  I 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( err )/ log(ℎ).Again we can see that when  is above 2, accurate SIFs can be obtained.For the range of penalty factor in 10 3 ∼10 5 , almost the same results are generated.

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 (Figure 6(a)).In order to compare the results with the reference, the following geometry and loading parameters are used:  = ℎ = 4,  = 1, and crack length  = 2.A 80×80 mesh is used.In this case,  = 14 together with   = 3 are used to calculate the coefficients of Williams' expansion.The penalty factor is assumed to be 1000.Figure 6(b) shows the   -stress nephogram, where the red lines represent the common interface of Ω  and Ω  .Table 1 shows the coefficients  1 - 5 ,  1 , and  3 - 5 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  1 - 5 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 halfcrack model with mixed mode leads to wrong results for the coefficients of higher-order parameters.Now consider the centre-cracked plate under pure shearing (CCPS).It is a typical example of mode II.The geometrical parameters of the CCPS (Figure 7(a)), the mesh, and material properties are the same as CCPT.Here, a uniformly distributed shear force  = 1 is imposed.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.

A Star-Shaped
Crack in a Square Plate.At last a starshaped 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 Figure 8.
In this example, the plate dimension is fixed to be  = 5, and material constants are Young's modulus  = 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 / = 0.1.The cracks are moved with a small distance to avoid special case caused by symmetry.Here  = 5 and   = 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.

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.

Figure 1 :
Figure 1: Decomposition of the domain into a standard domain and an analytical patch.

Figure 2 :
Figure 2: A crack cuts one mathematical patch into two physical patches.

Figure 4 :Figure 5 :
Figure 4: Normalized  I as function of the size of the patch for different values of .

Figure 7 (
b) gives the   -stress nephogram and again the red lines represent the interfaces.The coefficients  1 - 5 ,  1 , and  3 - 5 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.
The normalized stress intensity factors at tips  and  are defined as

Figure 8 :
Figure 8: A star-shaped crack in a square plate under biaxial tension.

Figure 9 :
Figure 9: Star-shaped crack and the hierarchical meshes: (a) global view; (b) local view with analytical patches.