Independent cover meshless particle method for complex geotechnical engineering

Abstract

A new Independent Cover Meshless Particle (ICMP) method is proposed for the analysis of complex geotechnical engineering. In the ICMP method, the independent rectangular cover regardless of the shape of the analysis model is employed as the influence domain of each discrete node, the general polynomial is employed as the meshless interpolation function of the independent nodal cover, and the Cartesian Transformation Method (CTM) is used for the numerical integration of the nodal covers cut by material interfaces, joints, cracks and faults. The present method has a simple formulation and a low computational cost, and is easy for the numerical analysis and modeling of complex geotechnical engineering. Several typical numerical examples are presented to demonstrate the accuracy and robustness of the proposed method.

Keywords

meshless method / particle method / independent cover / CTM / geotechnical engineering

Introduction

Large-scale geotechnical engineering generally passes through different stratums, involves the developed faults, joints and cracks, and encounters engineering geological problems such as high-pressure fissure water, high ground pressure and rock burst. It also involves complex constructions of excavation, pouring and supporting. Taking into account the complex terrain conditions, geometric shapes and boundary conditions, the traditional numerical methods find great challenge and difficulties in the safety assessment and numerical analysis of large-scale geotechnical engineering, especially in mesh generation and mesh moving for crack growth.

In order to resolve these difficulties, in recent decades, significant development and progress have been achieved for the research on new numerical methods, such as the Meshless Method (MM) [¹–⁵], Numerical Manifold Method (NMM) [⁶–¹⁰], Extended Finite Element Method (XFEM) [¹¹–¹³], Generalized Finite Element Method (GFEM) [¹⁴–¹⁶], Peridynamics (PD) [¹⁷–¹⁹]. These methods have showed their potential and advantages, and have been developed and widely used in different areas. However, there still exist obvious deficiencies in these methods, which seriously hinder their further applications in practical engineering problems. For example, the meshless methods have advantages in avoiding of element limitation to interpolation functions based on arbitrarily distributed nodes. However, the meshless approaches are generally rational functions and confront with the difficulties and problems of high computational cost, complex numerical implementation, and the shortcomings of arbitrariness in the selection of nodal influence covers. The meshless particle method [²⁰,²¹] based on meshless interpolation function shows more flexible by using the special point integration, but it has yet to overcome the above-mentioned fundamental difficulties in meshless interpolations. The NMM employs the mathematical cover independent of analysis area to define the computation accuracy, and thus gets rid of the difficulty of mesh generation for material interface, faults, joints and cracks, and is convenient to simulate complex moving boundary problems. However, the absence of a systematic cover generation theory and algorithms has long become the bottle neck of the NMM and impeded the progress of the research and application of the NMM [²²]. The similar difficulties are found in XFEM and GFEM.

In this paper, a new meshless particle method based on independent nodal cover is developed for the analysis of complex geotechnical engineering. In the ICMP (Independent Cover Meshless Particle) method, the independent rectangular cover regardless of the shape of analysis model is defined to describe the influence domain of each discrete node, the general polynomial is employed as the meshless interpolation function of the independent nodal cover, and the Cartesian Transformation Method (CTM) is used for the numerical integration of the nodal covers cut by material interfaces, joints and faults. The new method has the advantages of simple formulation and implementation, and low computation burden, and it can be simple and convenient to build 2D/3D numerical model for complex geotechnical problems. The ICMP overcomes most of the difficulties encountered in the above mentioned numerical methods, and shows a potential prospect for the simulation of failure process of geotechnical structure and dynamic propagation of the joints and cracks.

The paper is organizes as follows. Section 2 introduces the basic concepts of the ICMP. Section 3 provides a detailed description of the implementation of the numerical integration for independent nodal covers with the usage of the CTM. Section 4 presents four typical numerical examples to demonstrate the accuracy and robustness of the proposed method. Finally, the key points and main conclusions are summarized in Section 5.

Introduction of the basic concepts of the ICMP

Meshless interpolation for each independent nodal cover

Consider an arbitrary two-dimensional analysis domain W as shown in Fig. 1. There are a joint and a material interface represented by lines denoted as 1 and 2, respectively. A rectangular mesh regardless of analysis model is used to cover the analysis domain W as shown in Fig. 2. The domain W is discretized by a set of uniformly distributed nodes i (i = 1‒20), and each rectangular cover C_i is defined as the influence domain (or cover) of the discrete node i. The coordinates of node i are taken as the central position of the rectangular cover C_i.

In the ICMP, a particle approach is used to simplify the topological operation. The basic idea is to use a group of uniformly distributed nodes to discretize the arbitrary analysis domain. Generally, a single particle/node represents an influence domain assigned with central position, rectangular cover C_i and material parameters. As seen, this approach would greatly reduce the complexity of the algorithm for mesh generation and simplify the analysis process of the complex geotechnical problems.

Define the material interface, joint and boundary lines as the physical lines. For the domains doesn’t cut by physical lines, such as cover C_{\mathrm{10}}, the approximation for cover C_{\mathrm{10}} is defined as

(1){\overline{\mathbf{u}}}^{\mathrm{10}}(\overline{x},\overline{y})={\mathbf{\Phi }}^{\mathrm{10}}(\overline{x},\overline{y})\cdot {\mathbf{a}}^{\mathrm{10}},

where {\overline{\mathbf{u}}}^{\mathrm{10}}={[\overline{u}(\overline{x},\overline{y}),\overline{v}(\overline{x},\overline{y})]}^{\mathrm{T}} is the displacement approximation of cover C_{\mathrm{10}} along the x and y axes, {\mathbf{a}}^{\mathrm{10}}={[a_{\mathrm{1}},a_{\mathrm{2}},\dots ,a_n]}^{\mathrm{T}} is the DOF (degree of freedom) of cover C_{\mathrm{10}}, {\mathbf{\Phi }}^{\mathrm{10}} is the polynomial basis function where

(2){\mathbf{\Phi }}^i(\overline{x},\overline{y})=\left[\begin{array}{l}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overline{x}}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overline{y}}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdots \\ \mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overline{x}}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overline{y}}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cdots \end{array}\right],

where {\overline{x}}_{\mathrm{0}}=\overline{x}-{\overline{x}}_{\mathrm{10}}, {\bar{y}}_{\mathrm{0}}=\bar{y}-{\bar{y}}_{\mathrm{10}}, ({\bar{x}}_{\mathrm{10}},{\bar{y}}_{\mathrm{10}}) is the coordinate of node 10. In this work, a two-order polynomial is employed as the basis function.

The interpolation function defined in Eq. (1) is regardless of the shape of the analysis domain and is actually a kind of meshless approximation based on discrete nodes.

As shown in Fig. 3, physical lines 1 and 2 partition cover C_{\mathrm{15}} into four independent sub-domains {\mathrm{15}}^{\mathrm{1}}, {\mathrm{15}}^{\mathrm{2}}, {\mathrm{15}}^{\mathrm{3}} and {\mathrm{15}}^{\mathrm{4}}. Suppose that physical line 1 is described by function F_{\mathrm{1}}(\overline{\mathbf{x}})=\mathrm{0}, and physical line 2 is describe by function F_{\mathrm{2}}(\overline{\mathbf{x}})=\mathrm{0}, denoted as F_{\mathrm{10}} and F_{\mathrm{20}} respectively [²³]. Physical line 1 partitions cover C_{\mathrm{15}} into two parts F_{\mathrm{1}}(\overline{\mathbf{x}})>0 and F_{\mathrm{1}}(\overline{\mathbf{x}})<0, denoted as F_{\mathrm{1A}} and F_{\mathrm{1}\mathrm{B}}, respectively. Physical line 2 partitions the cover C_{\mathrm{15}} into two parts F_{\mathrm{2}}(\overline{\mathbf{x}})>0 and F_{\mathrm{2}}(\overline{\mathbf{x}})<0, denoted as F_{\mathrm{2}\mathrm{A}} and F_{\mathrm{2}\mathrm{B}}, respectively. By using the symbol functions F_{\mathrm{1}} and F_{\mathrm{2}}, the sub-domain {\mathrm{15}}^{\mathrm{1}}, {\mathrm{15}}^{\mathrm{2}}, {\mathrm{15}}^{\mathrm{3}} and {\mathrm{15}}^{\mathrm{4}} can be record as

(3)\left[\begin{array}{l}{\mathrm{15}}^{\mathrm{1}}\\ {\mathrm{15}}^{\mathrm{2}}\\ {\mathrm{15}}^{\mathrm{3}}\\ {\mathrm{15}}^{\mathrm{4}}\end{array}\right]=\left[\begin{array}{l}{F}_{\mathrm{1A}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}_{\mathrm{2}\mathrm{A}}\\ F_{\mathrm{1}\mathrm{B}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}_{\mathrm{2}\mathrm{A}}\\ F_{\mathrm{1B}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}_{\mathrm{2}\mathrm{B}}\\ F_{\mathrm{1A}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{F}_{\mathrm{2}\mathrm{B}}\end{array}\right].

It can be seen that, the complex geometry algorithm partitioning the sub-domain and the time–consuming judgment of the point–polygon relations can be avoided with the definition of Eq. (3). For example, the sub-domain number {\mathrm{15}}^{\mathrm{1}} of integration point K in Fig. 3 can be easily identified by values (F_{\mathrm{1A}}, F_{\mathrm{2}\mathrm{A}}) of the symbol functions F_{\mathrm{1}} and F_{\mathrm{2}}.

Each of the sub-domains {\mathrm{15}}^{\mathrm{1}}, {\mathrm{15}}^{\mathrm{2}}, {\mathrm{15}}^{\mathrm{3}} and {\mathrm{15}}^{\mathrm{4}} has independent displacement and deformation, the displacement function of sub-domain {\mathrm{15}}^i can be defined as {\overline{\mathbf{u}}}^{{\mathrm{15}}^{\mathrm{1}}}, {\overline{\mathbf{u}}}^{{\mathrm{15}}^{\mathrm{2}}}, {\overline{\mathbf{u}}}^{{\mathrm{15}}^{\mathrm{3}}}and {\overline{\mathbf{u}}}^{{\mathrm{15}}^{\mathrm{4}}} by Eq. (1).

Substituting Eq. (1) into the strain–displacement relations of the linear elastic problem undergoing infinitesimal deformation, we can get the strain of cover C_{\mathrm{i}} as

(4)\mathbf{\epsilon }=\mathbf{L}{\overline{\mathbf{u}}}^{C_{\mathrm{i}}}=\mathbf{L}{\mathbf{\Phi }}^{C_{\mathrm{i}}}{\mathbf{a}}^{C_{\mathrm{i}}}=\mathbf{B}{\mathbf{a}}^{C_{\mathrm{i}}},

where \mathbf{\epsilon } is the strain vector and \mathbf{{\rm B}} is the strain matrix, where

(5)\mathbf{\epsilon }={\left[{\epsilon }_{\overline{x}}{\epsilon }_{\overline{y}}{\gamma }_{\overline{x}\hspace{0.17em}\overline{y}}\right]}^{\mathrm{T}},

(6)\mathbf{B}=\mathbf{L}{\mathbf{\Phi }}^{C_{\mathrm{i}}},

where \mathbf{L} is the differential operator of the plane problem and

(7)\mathbf{L}=\left[\begin{array}{ccc}\partial /\partial \overline{x}& \square & 0\\ 0& \square & \partial /\partial \overline{y}\\ \partial /\partial \overline{y}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& \square & \partial /\partial \overline{x}\end{array}\right].

For the isotropic linear elastic material, the stress of an arbitrary point in cover C_i can be expressed as

(8)\mathbf{\sigma }=\mathbf{D}\mathbf{B}{\mathbf{a}}^{C_i},

where \mathbf{D} is the elastic matrix, and

(9)\mathbf{D}={\displaystyle \frac{E_{\mathrm{0}}}{\mathrm{1}-v_{\mathrm{0}}^{\mathrm{2}}}}\left[\begin{array}{l}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{v}_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\\ v_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\\ \mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\mathrm{1}-v_{\mathrm{0}}}{\mathrm{2}}}\end{array}\right],

where E_{\mathrm{0}}=E is the Young’s modulus and v_{\mathrm{0}}=v is the Poisson’s ration for plane stress case.

Thus, the strain energy of cover C_{\mathrm{i}} can be expressed as

(10){\mathrm{\Pi }}^C={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{({\mathbf{a}}^{C_i})}^{\mathrm{T}}{\displaystyle {\iint }_{C_i}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{B}t\mathrm{d}\overline{x}\mathrm{d}\overline{y}{\mathbf{a}}^{C_i}},

where t is the thickness of the analysis domain.

Connection between the adjacent independent nodal covers

With the definition of the meshless interpolation in Eq. (1), the influence domains of the adjacent nodes have their own independent freedoms and deformations, for example, covers C_{\mathrm{10}} and C_{\mathrm{11}}^{\mathrm{1}}, and covers C_{\mathrm{11}}^{\mathrm{2}} and C_{\mathrm{11}}^{\mathrm{3}}, as shown in Figs. 4 and 5, respectively. Considering the two adjacent covers C_{\mathrm{11}}^{\mathrm{2}} and C_{\mathrm{11}}^{\mathrm{3}}, the displacements of the share boundary between covers C_{\mathrm{11}}^{\mathrm{2}} and C_{\mathrm{11}}^{\mathrm{3}} by Eq. (1) are totally discontinuous, but the displacements along the share boundary should actually be continuous.

Suppose that C_{\mathrm{11}}^{\mathrm{2}} and C_{\mathrm{11}}^{\mathrm{3}} are connected by a fictitious thin-layer (or link element) L_C with width b and length l as shown in Fig. 5, where b\ll l [²²]. Because b\ll l in link element L_C, the deformation along the x direction can be ignored in comparison with the deformation along the y direction, the strain and displacement can then be approximated as

(11){\epsilon }_y\approx {\displaystyle \frac{v^{{\mathrm{11}}^{\mathrm{2}}}-v^{{\mathrm{11}}^{\mathrm{3}}}}{b}},\hspace{0.17em}{\gamma }_{yx}\approx {\displaystyle \frac{u^{{\mathrm{11}}^{\mathrm{2}}}-u^{{\mathrm{11}}^{\mathrm{3}}}}{b}},\hspace{0.17em}{\epsilon }_x\approx \mathrm{0.}

By substituting Eq. (1) into Eq. (11), we can get the strain of L_C as

(12){\mathbf{\epsilon }}^l={\displaystyle \frac{\mathrm{1}}{b}}{\mathbf{N}}^l{\mathbf{a}}^l,

where {\epsilon }^l={[{\gamma }_{yx},\hspace{0.17em}{\epsilon }_y]}^{\mathrm{T}}, {\mathbf{N}}^lis the shape function of link element L_C, and

(13){\mathbf{N}}^l={\mathbf{\lambda }}^l\left[{\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{2}}}\left({\overline{\mathbf{x}}}^p\right),\hspace{0.17em}{\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{3}}}\left({\overline{\mathbf{x}}}^p\right)\right],

where {\overline{\mathbf{x}}}^p is the coordinate of an arbitrary point P in link element L_C, {\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{2}}} is the shape function of point P in cover C_{\mathrm{11}}^{\mathrm{2}}, {\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{3}}} is the shape function of point P in cover C_{\mathrm{11}}^{\mathrm{3}}, {\mathbf{a}}^l is the freedom vector of link element L_C, and

(14){\mathbf{a}}^l=\left[\begin{array}{l}{\mathbf{a}}^{{\mathrm{11}}^{\mathrm{2}}}\\ {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{3}}}\end{array}\right],

where {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{3}}} is the DOF of cover C_{\mathrm{11}}^{\mathrm{3}}, {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{2}}} is the DOF of cover C_{\mathrm{11}}^{\mathrm{2}}, {\mathbf{\lambda }}^l is the coordinate transformation matrix, and

(15){\mathbf{\lambda }}^l=\left[\begin{array}{cc}\cos \hspace{0.17em}\theta & \sin \hspace{0.17em}\theta \\ -\sin \hspace{0.17em}\theta & \cos \hspace{0.17em}\theta \end{array}\right].

The stress of L_C can then be expressed as

(16){\mathbf{\sigma }}^l={\mathbf{D}}^l{\mathbf{\epsilon }}^l,

where

(17){\mathbf{\sigma }}^l={[{\tau }_{yx},{\sigma }_y]}^{\mathrm{T}},

(18){\mathbf{D}}^l=\left[\begin{array}{cc}{\mathit{G}}_{\mathrm{0}}\hspace{0.17em}& \mathrm{0}\\ \mathrm{0}& {\mathit{E}}_{\mathrm{0}}\end{array}\right],

where G_{\mathrm{0}}={\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+v)}} and E_{\mathrm{0}}={\displaystyle \frac{E}{\mathrm{1}-v^{\mathrm{2}}}} for plane stress case.

The width b should be taken the value much less than that of the length l for link element L_{\mathrm{C}} to achieve a higher solution accuracy. Numerical results indicate that the accuracy of the ICMP is satisfactory while b\le {\mathrm{10}}^{-\mathrm{3}}l. In this work, width b is taken as b={\mathrm{10}}^{-\mathrm{3}}l for all numerical examples.

Thus, the strain energy of link element L_C can be expressed as

(19){\mathrm{\Pi }}^l={\displaystyle \frac{\mathrm{1}}{\mathrm{2}b}}{({\mathbf{a}}^l)}^{\mathrm{T}}{\displaystyle {\int }_{-l/\mathrm{2}}^{l/\mathrm{2}}{({\mathbf{N}}^l)}^{\mathrm{T}}{\mathbf{D}}^l{\mathbf{N}}^{l}t\mathrm{d}s{\mathbf{a}}^l}.

Joint element

There is a joint between cover C_{\mathrm{11}}^{\mathrm{1}} and cover C_{\mathrm{11}}^{\mathrm{2}} as shown in Fig. 6. A Goodman-type element L_S is used to simulate the joint. With a similar derivation process, the strain energy of element L_S is expressed as

(20){\mathrm{\Pi }}^S={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{({\mathbf{a}}^S)}^{\mathrm{T}}{\displaystyle {\int }_{-l/\mathrm{2}}^{l/\mathrm{2}}{({\mathbf{N}}^S)}^{\mathrm{T}}{\mathbf{D}}^S{\mathbf{N}}^{S}t\mathrm{d}\overline{x}{\mathbf{a}}^S},

where {\mathbf{N}}^S is the shape function of element L_S, and

(21){\mathbf{N}}^S={\mathbf{\lambda }}^S\left[{\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{1}}}\left({\overline{\mathbf{x}}}^p\right),\hspace{0.17em}-{\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{2}}}\left({\overline{\mathbf{x}}}^p\right)\right],

{\overline{\mathbf{x}}}^p is the coordinate of an arbitrary point P in link element L_S, {\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{2}}} is the shape function of point P in cover C_{\mathrm{11}}^{\mathrm{2}}, {\mathbf{\Phi }}^{{\mathrm{11}}^{\mathrm{1}}} is the shape function of point P in cover C_{\mathrm{11}}^{\mathrm{1}}, {\mathbf{\lambda }}^S is the coordinate trasformation matrix similar to Eq. (15), and

(22){\mathbf{a}}^S=\left[\begin{array}{l}{\mathbf{a}}^{{\mathrm{11}}^{\mathrm{3}}}\\ {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{2}}}\end{array}\right],

where {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{1}}} is the DOF of cover C_{\mathrm{11}}^{\mathrm{1}}, {\mathbf{a}}^{{\mathrm{11}}^{\mathrm{2}}} is the DOF of cover C_{\mathrm{11}}^{\mathrm{2}}, and

(23){\mathbf{D}}^S=\left[\begin{array}{cc}{k}_S& \mathrm{0}\\ \mathrm{0}& k_n\end{array}\right],

where k_S and k_n are the normal stiffness and shear stiffness of element L_S, respectively. In general, the coefficients k_S and k_n of the joint should be analyzed according to the actual measurement or field experiment. Please refer to [²⁴] for the detailed calculation method of the joint stiffness.

Imposition of displacement boundary condition

Suppose that the displacement of boundary line e-f of cover C_{\mathrm{13}} is fixed along the x and y directions as shown in Fig. 7. The displacement boundary condition can also be imposed by using a fictitious thin-layer element L_b similar to section 2.2. The strain energy of fictitious element L_b is expressed as

(24){\mathrm{\Pi }}^b={\displaystyle \frac{\mathrm{1}}{\mathrm{2}b}}{({\mathbf{a}}^b)}^{\mathrm{T}}{\displaystyle {\int }_{-l/\mathrm{2}}^{l/\mathrm{2}}{({\mathbf{N}}^b)}^{\mathrm{T}}{\mathbf{D}}^b{\mathbf{N}}^{b}t\mathrm{d}\overline{x}{\mathbf{a}}^b},

where

(25){\mathbf{N}}^b={\mathbf{\lambda }}^b{\mathbf{\Phi }}^{\mathrm{13}}\left({\overline{\mathbf{x}}}^p\right),

{\overline{\mathbf{x}}}^p is the coordinate of an arbitrary point P in element L_b, {\mathbf{\Phi }}^{\mathrm{13}} is the shape function of point P in cover C_{\mathrm{13}}, {\mathbf{\lambda }}^b is the coordinate transformation matrix similar to Eq. (15), {\mathbf{a}}^{\mathrm{13}} is the DOF of cover C_{\mathrm{13}}, {\mathbf{D}}^b is the elastic matrix similar to Eq. (18).

Imposition of load boundary condition

There are uniformly distributed load {\mathbf{f}}_{\mathrm{0}}={[f\overline{x}(x)\hspace{0.17em}f\overline{y}(x)]}^{\mathrm{T}} along the boundary of cover C_{\mathrm{2}} as shown in Fig. 8. The potential energy of force {\mathbf{f}}_{\mathrm{0}} can be expressed as

(26){\mathrm{\Pi }}^d=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{({\mathbf{a}}^d)}^{\mathrm{T}}{\displaystyle {\int }_{\overline{x}}{({\mathbf{N}}^d)}^{\mathrm{T}}{\mathbf{f}}_{\mathrm{0}}td\overline{x}},

where {\mathbf{N}}^d={\mathbf{\lambda }}^d{\mathbf{\Phi }}^{\mathrm{2}}, {\mathbf{\Phi }}^{\mathrm{2}} is the shape function of cover C_2, {\mathbf{\lambda }}^d is the coordinate transformation matrix similar to Eq. (15), and {\mathbf{a}}^d is the DOF of cover C_2.

Other types of displacement and force boundary conditions can be derived through the same way.

Equilibrium equation

By using the Eqs. (10), (19), (20), (24) and (26), the total potential energy of domain \mathrm{\Omega } is written as

(27)\mathrm{\Pi }=\mathrm{\Sigma }({\mathrm{\Pi }}^c+{\mathrm{\Pi }}^l+{\mathrm{\Pi }}^S+{\mathrm{\Pi }}^b+{\mathrm{\Pi }}^d),

invoking \delta \mathrm{\Pi }=\mathrm{0} results in the following discrete equation

(28){\displaystyle \frac{\partial \mathrm{\Pi }}{\partial a}}=\mathrm{\Sigma }({\mathbf{K}}^C{\mathbf{a}}^C+{\mathbf{K}}^l{\mathbf{a}}^l+{\mathbf{K}}^S{\mathbf{a}}^S+{\mathbf{K}}^b{\mathbf{a}}^b-{\mathbf{F}}^d),

where,

(29){\mathbf{K}}^C={\displaystyle {\iint }_{C_i}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{B}t\mathrm{d}\overline{x}\mathrm{d}\overline{y}},

(30){\mathbf{K}}^l={\displaystyle \frac{\mathrm{1}}{b}}{\displaystyle {\int }_L{({\mathbf{N}}^l)}^{\mathrm{T}}{\mathbf{D}}^l{\mathbf{N}}^{l}t\mathrm{d}s},

(31){\mathbf{K}}^S={\displaystyle {\int }_L{({\mathbf{N}}^S)}^{\mathrm{T}}{\mathbf{D}}^S{\mathbf{N}}^{S}t\mathrm{d}\overline{x}},

(32){\mathbf{K}}^b={\displaystyle \frac{\mathrm{1}}{b}}{\displaystyle {\int }_L{({\mathbf{N}}^b)}^{\mathrm{T}}{\mathbf{D}}^b{\mathbf{N}}^{b}t\mathrm{d}\overline{x}},

(33){\mathbf{F}}^d={\displaystyle \int {({\mathbf{N}}^d)}^{\mathrm{T}}{\mathbf{f}}_{\mathrm{0}}t\mathrm{d}\overline{x}}.

Assembling all of the above stiffness matrices and force vectors, the equilibrium equation for the domain \mathrm{\Omega } is then expressed as

(34)\mathbf{K}\cdot \mathbf{a}=\mathbf{F},

where \mathbf{K} and \mathbf{F} are the equivalent stiffness matrix and force vector, respectively. \mathbf{a} is the vector of DOFs to be solved.

In the ICMP method, an arbitrary analysis domain is discretized by a set of uniformly distributed nodes, and an independent rectangular cover is employed as the influence domain for each discrete node. Although the fictitious thin layer element has to be used to impose the continuous condition between the adjacent independent nodal covers, the present method still keeps the characteristics of the meshless methods and has the advantages of simple formulation and implementation, and low computation burden, and it can be simple and convenient to build 2D/3D numerical model for complex geotechnical problems.

Integration of the nodal covers

Cartesian transformation method for evaluation of domain integrals

Consider a cover that isn’t intersected by physical lines, such as cover C_{\mathrm{10}} shown in Fig. 2. It’s easy to compute the numerical integration of cover C_{\mathrm{10}} by using conventional Gaussian quadrature method. But in general, cover C_i is partitioned into sub-domains by physical lines, for example, cover C_{\mathrm{15}} is partitioned into 4 sub-domains {\mathrm{15}}^{\mathrm{1}}, {\mathrm{15}}^{\mathrm{2}}, {\mathrm{15}}^{\mathrm{3}} and {\mathrm{15}}^{\mathrm{4}} as shown in Fig. 3. Each sub-domain of cover C_{\mathrm{15}} has independent DOF, and the Gauss quadrature is not directly applicable any more for the sub domains with the shape of arbitrary polygon. Similar problem occurs in the XFEM and NMM. Common method solving the difficulty is to further subdivide the sub-domains into sub-triangles as shown in Fig. 9, in which the hammer quadrature rule can be implemented directly. The principle of the method is simple, but its disadvantages are also obvious. The subdividing of the triangles not only increases the mesh-dependency of the method, but also leads to the complexity of geometric subdivision algorithm, especially for three-dimensional integral domain with arbitrary geometrical shapes.

In this study, the CTM, which is proposed by Khosravifard and Hematiyan [²⁵–²⁸], is used to integrate the stiffness matrix of cover C_{\mathrm{15}}. The CTM avoids the complex geometric subdivision of the integral domain, and is simple and convenient in operation. Next, we will take cover C_{\mathrm{15}} in Fig. 10 as an example to introduce the implementation of the CTM.

Consider a general integral over cover C_{\mathrm{15}}

(35)\mathrm{I}={\displaystyle {\int }_{\mathrm{\Omega }}{f}_i(\overline{x},\overline{y})\mathrm{d}\mathrm{\Omega }},

where f_i(\bar{x},\bar{y}) is the function to be integrated in cover C_{\mathrm{15}}, \mathrm{\Omega } represent the integration domain. According to the CTM, this domain integral can be transformed into two 1D integrals as

(36)\mathrm{I}={\displaystyle {\int }_{{\Gamma }_{\mathrm{2}}}g(\overline{y})\mathrm{d}\overline{y},}

and

(37)g(\overline{y})={\displaystyle {\int }_a^b{f}_i(\overline{x},\overline{y})\mathrm{d}\overline{x},}

where {\Gamma }_{\mathrm{2}}:{\bar{x}}_{\mathrm{4}}=b is the right boundary of cover C_{\mathrm{15}}, a and b are the coordinates of the vertices of cover C_{\mathrm{15}} along {\bar{x}}_{\mathrm{4}} direction.

Due to the transformation, we can use the Gauss quadrature (GQ) to compute the integration of cover C_{\mathrm{15}}, the procedure is as following:

1) Recording the location of the vertices of rectangular domain, intersection points between physical lines and domain boundary, and the intersection points between physical lines as \mathbf{X}=\{{\overline{x}}_{\mathrm{1}},{\overline{x}}_{\mathrm{2}},\dots \}, \mathbf{Y}=\{{\overline{y}}_{\mathrm{1}},{\overline{y}}_{\mathrm{2}},\dots \}, the boundary {\Gamma }_{\mathrm{2}} in Fig. 10 is partitioned into 3 portions as follows:

(38)L_{\mathrm{1}}:{\bar{y}}_{\mathrm{2}}\le \bar{y}\le {\bar{y}}_{\mathrm{1}},\hspace{0.17em}{L}_{\mathrm{2}}:{\bar{y}}_{\mathrm{3}}\le \bar{y}\le {\bar{y}}_{\mathrm{2}},\hspace{0.17em}{L}_{\mathrm{3}}:{\bar{y}}_{\mathrm{4}}\le \bar{y}\le {\bar{y}}_{\mathrm{3}}.

2) Define the portions as integration intervals, and an m-point Gaussian quadrature can be used for each individual interval. For example, we set three integration points on integration interval L_2 along the \bar{y} direction, as shown in Fig. 10. The value of the integral can be computed as follows:

(39)I={\displaystyle {\int }_{{\overline{y}}_{\mathrm{4}}}^{{\overline{y}}_{\mathrm{1}}}g(\overline{y})\mathrm{d}\overline{y}}={\mathrm{\Sigma }}_i^n{\mathrm{\Sigma }}_j^m{J}_i{W}_{j}G({\eta }_i),

where G({\eta }_i)=g(\bar{y}), {\eta }_i are the Gaussian points, n=\mathrm{3} is the number of intervals, m=\mathrm{3} is the number of integration points, J_i (i = 1, 2, 3) are the Jacobian values with J_{\mathrm{1}}={\displaystyle \frac{{\bar{y}}_{\mathrm{1}}-{\bar{y}}_{\mathrm{2}}}{\mathrm{2}}}, J_{\mathrm{2}}={\displaystyle \frac{{\bar{y}}_{\mathrm{2}}-{\bar{y}}_{\mathrm{3}}}{\mathrm{2}}} and J_{\mathrm{3}}={\displaystyle \frac{{\bar{y}}_{\mathrm{3}}-{\bar{y}}_{\mathrm{4}}}{\mathrm{2}}}, and W_j denotes the integration weight of the GQ method.

3) Elicit horizontal lines from the integration points along boundary {\Gamma }_{\mathrm{2}}, and define the horizontal lines as integration rays shown in Fig. 10. Each integration ray is divided into several integration intervals by boundary lines and physical lines. For example, there are three integration rays associated with interval L_2 along the y direction as shown in Fig. 10, a typical integration ray y=y({\eta }_i) is divided into three integration intervals, we can record the intervals as:

(40)S_{\mathrm{1}}:{\bar{x}}_{\mathrm{1}}\le \bar{x}\le {\bar{x}}_{\mathrm{2}},\hspace{0.17em}{S}_{\mathrm{2}}:{\bar{x}}_{\mathrm{2}}\le \bar{x}\le {\bar{x}}_{\mathrm{3}},\hspace{0.17em}{S}_{\mathrm{3}}:{\bar{x}}_{\mathrm{3}}\le \bar{x}\le {\bar{x}}_{\mathrm{4}}.

4) To compute the value of the integral in Eq. (35), the value of f_i(\bar{x},\bar{y}) must be calculated for each ray. Symbol function values in Eq. (3) is employed to identify the cover numbers of the integration points, and the m-point Gaussian quadrature is used for each individual interval. For example, integration intervals S_{\mathrm{1}}, S_{\mathrm{2}} and S_{\mathrm{3}} along the \bar{x} direction belonging to cover C_{\mathrm{15}}^{\mathrm{3}}, C_{\mathrm{15}}^{\mathrm{2}} and C_{\mathrm{15}}^{\mathrm{1}}, respectively, and m integration points are set along each interval S_i. The integration of f_i(\bar{x},\bar{y}) along \bar{x} direction can be calculated as follows:

(41)I={\displaystyle {\int }_a^b{f}_i(\overline{x},\overline{y})}\mathrm{d}\overline{x}={\mathrm{\Sigma }}_{r=\mathrm{1}}^n\left({\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{1}}{f}_i(\overline{x},\overline{y})}{J}_r\mathrm{d}\overline{x}\right)={\mathrm{\Sigma }}_r^n{\mathrm{\Sigma }}_S^m{J}_r{W}_S{F}_i{\mathrm{(\zeta }}_S),

where F_i{\mathrm{(\zeta }}_S)=f_i(\bar{x},\bar{y}), {\zeta }_S are the Gaussian points, n=\mathrm{3} is the number of intervals, m=\mathrm{3} is the number of integration points in the individual integration interval, J_r={\displaystyle \frac{{\bar{x}}_{i+\mathrm{1}}-{\bar{x}}_i}{\mathrm{2}}} is the Jacobian, and W_S denotes the integration weight of the GQ method.

Numerical implementation of the CTM

As can be observed from Eqs. (39) and (41), there have two Jacobians and two integration weights corresponding to each integration point. For the domain integral with different integrands in a same domain, the position of the integration points and the product of the associated Jacobians and weight of each point are stored in individual arrays. Then the interval is evaluated by the scalar product of the values of the integrand at integration points. The following relations clarify the procedure.

(1) Location of integration points are found and stored in an array.

(42)X=\left\{\begin{array}{l}{X}_{\mathrm{1}}\\ X_{\mathrm{2}}\\ ⋮\\ X_n\end{array}\right\},

where {\mathbf{X}}_i=\left[{\overline{x}}_i,{\overline{y}}_i\right], n is the number of the integration points.

(2) Jacobians and weights corresponding to each integration point are computed and the products of them are stored in an array.

(43){\mathbf{W}}^{\mathrm{2}D}=\{W_i^{\mathrm{2}D}\},\hspace{0.17em}{W}_i^{\mathrm{2}D}=W^{\mathrm{2}D}(X_i)=J_i^x\times J_i^y\times w_i^x\times w_i^y,

where J_i^x and J_i^y are the corresponding Jacobians for intervals along \bar{x} direction and \bar{y} direction, and w_i^x and w_i^y are the integration weights associated with the integration points of the intervals.

(3) The values of the integrand at integration points are computed

(44)\mathbf{F}=\left\{\begin{array}{l}{f}_{\mathrm{1}}\\ f_{\mathrm{2}}\\ ⋮\\ f_n\end{array}\right\},\hspace{0.17em}\mathrm{where}\hspace{0.17em}\hspace{0.17em}{f}_i=f({\mathbf{X}}_i).

(4) The value of the integral is computed by the scalar product of the two vectors \mathbf{F} and {\mathbf{W}}^{\mathrm{2}D}

(45)I={\mathbf{W}}^{\mathrm{2}D}\cdot \mathbf{F}={\mathrm{\Sigma }}_{i=\mathrm{1}}^N{W}^{\mathrm{2}D}({\mathbf{X}}_i)\times f({\mathbf{X}}_i).

Numerical examples

Cook beam

Cook beam subject to a uniform shear load F=\mathrm{1}/\mathrm{16} along the \bar{y} direction is considered as shown in Fig. 11. Plane stress condition is assumed. The material constants are E=\mathrm{1.0} and \mu =\mathrm{1}/\mathrm{3}. The exact solution for the displacement of point A is given by Cook et al. [²⁹].

The beam is discretized by 320 uniform nodes, as shown in Fig. 12. The vertical displacement V_A by the ICMP is 23.828, which is 0.551% away from the theoretical prediction of 23.96.

The discrete model with 20, 80, 320, 720 and 1280 distributed nodes are used to test the convergence. The errors for the vertical displacement V_Aversus node numbers are plotted in Fig. 13. As seen, the errors of the displacement solution rapidly convergence to the theoretical solution with the increasing of the node number.

Linear beam with a weak interlayer

Consider a linear beam with a weak interlayer, as shown in Fig. 14. The dimensions of the beam are L= 8.0 m and W= 1.0 m, the interlayer is located at \bar{x}= 3.1 m. a is the distance from the interlayer to the left side, b is the distance from the interlayer to the right side. The beam is subjected to a uniform traction \sigma =\mathrm{10}\hspace{0.17em}\mathrm{MPa} on the right side along the \bar{x} direction. The elastic material parameters of the beam are E= 10.0 MPa and \mu =\mathrm{0.25}. The problem is solved under a plane stress condition.

The problem is modeled using 128 uniformly distributed nodes, as shown in Fig. 14(b). The normal stiffness coefficient k_n is varied from 1.0 to 1000.0 MPa, and the tangential stiffness coefficient k_S is 1.0 MPa. The horizontal displacements of the midpoint C at the right end by the ICMP are compared with that of the theoretical solutions in Table 1. As seen, the ICMP shows a high accuracy for the this problem.

One-dimensional biomaterial bar

A one-dimensional (1-D) biomaterial bar as shown in Fig. 15 is studied. The dimension of the bar is (0, L). The material discontinuity is located at \bar{x}=\zeta. The elastic moduli are E_{\mathrm{1}} and E_{\mathrm{2}} for materials 1 and 2, respectively. The cross-section of the bar is assumed as A. The bar is fixed at the left end and subjected to a tensile loading F at the right end.

The tensile loading results a constant stress in the bar as

(46)\mathrm{\sigma }=F/A,\hspace{0.17em}\overline{x}\in (\mathrm{0},L).

According to Hook’s law, the strain field in the bar is obtained as

(47)\mathrm{\epsilon }=\{\begin{array}{l}{\displaystyle \frac{F}{E_{\mathrm{1}}A}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\overline{x}\in (0,\mathrm{\zeta }]\\ {\displaystyle \frac{F}{E_{\mathrm{2}}A}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\overline{x}\in [\mathrm{\zeta },L)\end{array}.

The displacement field in the bar is finally calculated as

(48)\mathrm{u}(x)=\{\begin{array}{l}{\displaystyle \frac{F}{E_{\mathrm{1}}A}}\overline{x},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\overline{x}\in (0,\zeta ]\\ {\displaystyle \frac{F}{E_{\mathrm{1}}A}}\mathrm{\zeta +}{\displaystyle \frac{F}{E_{\mathrm{2}}A}}(\overline{x}-\zeta ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\overline{x}\in [\zeta ,L)\end{array}.

The 1-D biomaterial bar problem can be simulated in a 2-D domain by simply letting the Poisson’s ratio be 0. Considering a 2-D domain \mathrm{\Omega }=(\mathrm{0},L)\times (\mathrm{0},L) with the material interface {\Gamma }^{\mathrm{1}-\mathrm{2}} located at \zeta \times (\mathrm{0},L). Denote the Young’s modulus and Passion’s ratio in {\mathrm{\Omega }}_{\mathrm{1}}=\zeta \times (\mathrm{0},L) by E_{\mathrm{1}}, {\mu }_1 and that in {\mathrm{\Omega }}_{\mathrm{2}}=(\mathrm{\zeta ,}L)\times (\mathrm{0},L) by E_{\mathrm{2}}, {\mu }_2. The ICMP model for a 2-D bimaterial problem is depicted in Fig. 16(a). The problem domain is discretized by 100 discrete nodes, as shown in Fig. 16(b). In the ICMP calculation, L= 1, \zeta = 0.25, A= 1, F = 1, E_{\mathrm{1}}= 1, E_{\mathrm{2}}=4, {\mu }_1={\mu }_{\mathrm{2}}=\mathrm{0}. The displacements along the \bar{x} direction computed by the ICMP are compared with that of the theoretical solutions given by Eq. (48), as shown in Fig. 17. The variation of the displacement across the material interface is clearly seen. The good agreement between the ICMP results and the theoretical solutions demonstrates the correctness and the accuracy of the proposed method.

Circular inclusion in an infinite plate under uniaxial tension

A circular inclusion in an infinite plate subjected to a uniaxial tension in the horizontal direction is studied. Due to the symmetry of the plate, a quarter of the plate is modeled as shown in Fig. 18(a). The dimension of the model is L=\mathrm{1}\hspace{0.17em}\mathrm{m}, the radius of the inclusion is a=\mathrm{0.1}\hspace{0.17em}\mathrm{m}, and the uniaxial tension is \sigma =\mathrm{1}\hspace{0.17em}\mathrm{kN/m}. Denote the Young’s modulus and Poisson’s ratio of the plate by E_{\mathrm{1}}, {\mu }_1 and that of the inclusion by E_{\mathrm{2}}, {\mu }_2, respectively. Lames elastic constants of the plate by {\lambda }_{\mathrm{1}}, {\gamma }_{\mathrm{1}} and that of the inclusion by {\lambda }_{\mathrm{2}}, {\gamma }_{\mathrm{2}}, respectively.

Two different cases regarding to {\lambda }_{\mathrm{2}}/{\lambda }_{\mathrm{1}} are considered. For soft inclusion problem {\lambda }_{\mathrm{2}}/{\lambda }_{\mathrm{1}}=\mathrm{0.5} and for hard inclusion problem {\lambda }_{\mathrm{2}}/{\lambda }_{\mathrm{1}}=\mathrm{2}. The theoretical solutions for the stress of the problem are given by An et al. [³⁰,³¹].

Case 1: {\lambda }_{\mathrm{2}}/{\lambda }_{\mathrm{1}}=\mathrm{0.5} (soft inclusion problem). The problem domain is discretized by 3600 uniform discrete nodes, as shown in Fig. 18(b). The results of the stress {\sigma }_{\bar{x}} by the ICMP are compared with that of the theoretical solutions in Fig. 19(a). It can be seen that the ICMP results match well with the theoretical solutions. The ICMP result of the stress {\sigma }_{\bar{x}} within the inclusion is averaged at 0.747, which is 0.400% away from the theoretical prediction of 0.75. The maximum stress {\sigma }_{\bar{x}}=\mathrm{1.490} appears at \bar{y} = 0.1 m, which is 0.655% away from the theoretical prediction of 1.500.

The discrete model with \mathrm{20}\times \mathrm{20}, \mathrm{30}\times \mathrm{30}, \mathrm{45}\times \mathrm{45} and \mathrm{60}\times \mathrm{60} distributed nodes are used for convergence studies. The errors for the stress {\sigma }_{\bar{x}} of the vertex C versus node numbers are plotted in Fig. 19(b). As seen, the errors of the stress solution rapidly convergence to the theoretical solution with the increasing of the node number.

Case 2: {\lambda }_{\mathrm{2}}/{\lambda }_{\mathrm{1}}=\mathrm{2} (hard inclusion problem). The problem domain is discretized by 8100 uniformly distributed nodes. The results of the stress {\sigma }_{\bar{x}} by the ICMP are compared with that of the theoretical solutions in Fig. 20(a). It can be seen that the ICMP results match well with the theoretical solutions. The ICMP result of the stress {\sigma }_{\bar{x}} within the inclusion is averaged at 1.203, which is 0.267% away from the theoretical prediction of 1.200. The Minimum stress {\sigma }_{\bar{x}}=\mathrm{0.605} appears at \bar{y} = 0.1 m, which is 0.833% away from the theoretical prediction of 0.600.

Again, the convergence test is carried out to validate the proposed method, four different meshes with discrete nodes numbers as \mathrm{30}\times \mathrm{30}, \mathrm{45}\times \mathrm{45}, \mathrm{60}\times \mathrm{60}, \mathrm{15}\times \mathrm{15} are investigated. The errors for the stress {\sigma }_{\bar{x}} of the vertex C versus node numbers are plotted in Fig. 20(b). Figure 20(b) shows that the errors of the stress solution rapidly convergence to the theoretical solution with the increasing of the node number.

As seen, the ICMP shows high accuracy for this material discontinuities problem.

Conclusions

A new ICMP method based on independent nodal covers is proposed for the analysis of complex geotechnical engineering. Five numerical examples are presented to demonstrate the robustness of the proposed method and evaluate its accuracy. The key points and main conclusions are summarized below:

(1)In the ICMP, the independent rectangular covers are defined to describe the influence domains of discrete nodes, and the general polynomial is employed as the meshless interpolation function of the independent nodal cover. The main advantage of the present ICMP is to get rid of the mesh dependency in interpolation, and is simple and convenient for 2D/3D modeling of complex geotechnical problem.

(2)Symbol function is utilized to identify the sub-domain numbers of the corresponding integration points, for avoiding the usage of the complex geometrical subdivision algorithm and the time-consuming judgment of the point-polygon relations. A fictitious link element is employed to impose the continuous between two adjacent independent covers and enforce the essential boundary conditions. Furthermore, it is very convenient for the ICMP to compute the integral value of general integral domain by using the CTM, whose formula derivation and numerical implementation are simple.

(3)Although the ICMP is only used for two-dimensional elastic analyses in the paper, the extension of the present work to three-dimensional analyses, fracture analyses and progressive failure analyses is straightforward and feasible, and they will be the further work of the present authors

References

[1] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: An overview and recent developments. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 3–47

[2] Li S, Liu W K. Meshfree and particle methods and their applications. Applied Mechanics Reviews, 2002, 55(1): 1–34

[3] Liu G R, Nourbakhshnia N, Zhang Y W. A novel singular ES-FEM method for simulating singular stress fields near the crack tips for linear fracture problems. Engineering Fracture Mechanics, 2011, 78(6): 863–876

[4] Cai Y C, Han L, Tian L G, Zhang L. Meshless method based on shepard function and partition of unity for two-dimensional crack problems. Engineering Analysis with Boundary Elements, 2016, 65: 126–135

[5] Asadpour A. Thermo-elastic extended meshfree method for fracture without crack tip enrichment. Frontiers of Structural and Civil Engineering, 2015, 9(4): 441–447

[6] Shi G H. Manifold method of material analysis. In: Transactions of the 9th Army Conference on Applied Mathematics and Computing. Minneapolis, MN: 1991. 57–76

[7] Shi G H. Simplex integration for manifold method, FEM, DDA and analytical analysis. In: Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media. Albuquerque, NM: 1996, 205–262

[8] Ma G W, An X M. The numerical manifold method: A review. International Journal of Clinical Medicine, 2010, 7: 1–32

[9] Sun L, Zhao G F, Zhao J. Particle manifold method (PMM): A new continuum-discontinuum numerical model for geomechanics. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(12): 1711–1736

[10] Wu J, Cai Y C. A partition of unity formulation referring to the NMM for multiple intersecting crack analysis. Theoretical and Applied Fracture Mechanics, 2014, 72: 28–36

[11] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[12] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131–150

[13] Daux C, Moes N, Dolbow J. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 2000, 48: 1741–1760

[14] Duarte C A, Oden J T. An h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 237–262

[15] Duarte C A, Hamzeh O N, Liszka T J, Tworzydlo W W. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Computer Methods in Applied Mechanics and Engineering, 2001, 190(15-17): 2227–2262

[16] Peng L J, Li X K, Ge X R. Numerical simulation of dynamic crack propagation with generalized finite element method. Chinese Journal of Rock Mechanics and Engineering, 2004, 23: 3132–3137

[17] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[18] Kilic B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials. Dissertation for the Doctoral Degree, the University of Arizona, 2008

[19] Kilic B, Madenci E. Structural stability and failure analysis using peridynamic theory. International Journal of Non-linear Mechanics, 2009, 44(8): 845–854

[20] Rabczuk T, Areias P M A, Belytschko T. A simplified mesh-free method for shear bands with cohesive surfaces. International Journal of Non-linear Mechanics, 2007, 69: 993–1021

[21] Rabczuk T, Song J H, Belytschko T. Simuliations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics, 2009, 76(6): 730–741

[22] Cai Y C, Zhu H H, Zhuang X Y. A continuous/discontinuous deformation analysis (CDDA) method based on deformable blocks for fracture modeling. Frontiers of Structural and Civil Engineering, 2013, 7(4): 369–378

[23] Cai Y C, Wu J, Atluri S N. A new implementation of the numerical manifold method (NMM) for the modeling of non-collinear and intersecting cracks. Computer Modeling in Engineering & Sciences, 2013, 92(1): 63–85

[24] Zhu B F. The Finite Element Method Theory and Applications. Beijing: China Water & Power Press, 2009

[25] Hematiyan M R. A general method for evaluation of 2D and 3D domain integrals without domain discretization and its application in BEM. Computational Mechanics, 2007, 39(4): 509–520

[26] Hematiyan M R. Exact transformation of a wide variety of domain integrals into boundary integrals in boundary element method. Communications in Numerical Methods in Engineering, 2008, 24: 497–521

[27] Khosravifard A, Hematiyan M R. A new method for meshless integration in 2D and 3D Galerkin meshfree methods. Engineering Analysis with Boundary Elements, 2010, 34(1): 30–40

[28] Khosravifard A, Hematiyan M R, Marin L. Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method. Applied Mathematical Modelling, 2011, 35(9): 4157–4174

[29] Cook R D, Malkus D S, Plesha M E. Concepts and Applications of Finite Element Analysis. 3rd ed. Hoboken, NJ: Wiley, 1989

[30] An X, Ma G, Cai Y, Zhu H, An X M, Ma G W, Cai Y C. A new way to treat material discontinuities in the numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47-48): 3296–3308

[31] An X M, Ma G W, Cai Y C, Zhu H. Arbitrary discontinuities in the numerical manifold method. International Journal of Computational Methods, 2011, 8(2): 315–347

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg