Research Article  Open Access
Hui Sun, Qin Yan, Bing Han, Shuchen Li, Xianda Feng, "Study on the Manifold Cover Lagrangian Integral Point Method Based on Barycentric Interpolation", Mathematical Problems in Engineering, vol. 2020, Article ID 6743204, 11 pages, 2020. https://doi.org/10.1155/2020/6743204
Study on the Manifold Cover Lagrangian Integral Point Method Based on Barycentric Interpolation
Abstract
To achieve numerical simulation of large deformation evolution processes in underground engineering, the barycentric interpolation test function is established in this paper based on the manifold cover idea. A largedeformation numerical simulation method is proposed by the double discrete method with the fixed Euler background mesh and moving material points, with discontinuous damage processes implemented by continuous simulation. The material particles are also the integration points. This method is called the manifold cover Lagrangian integral point method based on barycentric interpolation. The method uses the Euler mesh as the background integral mesh and describes the deformation behavior of macroscopic objects through the motion of particles between meshes. Therefore, this method can avoid the problem of computation termination caused by the distortion of the mesh in the calculation process. In addition, this method can keep material particles moving without limits in the set region, which makes it suitable for simulating large deformation and collapse problems in geotechnical engineering. Taking a typical slope as an example, the results of a slope slip surface obtained using the manifold cover Lagrangian integral point method based on barycentric interpolation proposed in this paper were basically consistent with the theoretical analytical method. Hence, the correctness of the method was verified. The method was then applied for simulating the collapse process of the side slope, thereby confirming the feasibility of the method for computing large deformations.
1. Introduction
The surrounding rocks of tunnels and other underground engineering structures mainly have three deformation failure mechanisms: rock burst, collapse, and large deformation [1]. Rock burst is the brittle failure of hard rocks under high geostress. Collapse and spalling are local deformation failures of surrounding rocks with a certain level of structural constraint. Large deformation of surrounding rocks is very common, especially in tunnels built in loess strata and highgeostress soft rocks [2–4]. Thus, understanding the evolution of large deformation of tunnels has extremely great importance for evaluating tunnel stability, investigating the root cause of collapses, and controlling engineering hazards.
Numerical computation methods, because of their proven effectiveness in tunnel stability analysis, have been developing rapidly and are applied widely in geotechnical engineering. The finite element method (FEM) and the finite difference method (FDM) are now the most common methods in geotechnical engineering for analyzing continuous deformations. The FEM, developed in the 1950s, is the earliest numerical method [5]. With the rapid development of computer technology, the FEM has been applied widely in many fields [6] and is now one of the most widely applied and established numerical methods. However, when used for computing large deformations in geotechnical engineering, the FEM may suffer mesh distortion, which may result in termination of computation. To simulate large sideslope collapses, many researchers have conducted indepth investigation of the constitutive model of geomaterials from the perspective of their solidliquid transformation [7, 8], but the problem of computation termination caused by meshdistortion remains unresolved. The FDM solves differential equations by reducing them to approximate difference equations, so it can, in some sense, account for the effect of mesh node movement on subsequent computations [9]. The FDM has been applied widely in geotechnical engineering [10, 11], but the magnitude of deformation is very limited.
To resolve this problem, improved FEMs, numerical manifold methods, discrete element methods, and meshless methods have been applied widely.
In the 2000s, the FEM with Lagrangian integration points (FEMLIP) was proposed in geophysics to represent mantle convection [12, 13]. The FEMLIP originates from the FEM and the particleincell method [14]. Bai et al. and Meng et al. [15, 16] used the FEMLIP to simulate mantle plumelithosphere interactions. The FEMLIP has also been used to simulate mudflows [17, 18]. Li et al. [19, 20] simulated mudflows using the FEMLIP and hydroelastoplastic models, which consider the rainfallinduced solidliquid transformation of geomaterials.
Shi [21] developed a numerical manifold method by discretizing materials using mathematical and physical manifolds and defining the relationship between manifold functions using weight functions. The method simulates the motion of material media by considering the interaction and relative motion between blocks and thus has good adaptability in large deformation computation. And since then, many improvements have been made to the method, expanding its scope of application. For example, Lu [22] developed a numerical manifold method with higher computational accuracy and larger scope of application by using highorder displacement manifold functions. Additionally, Wang et al. [23] developed a numerical manifold method based on the variational principles of parameters and applied it in an elastoplastic analysis of materials, and Zhang et al. [24] developed a secondorder numerical manifold method (including its computational codes) and applied it for simulating the entire process of structural failures. Dong et al. [25] developed an improved numerical manifold method that accounts for the reinforcing effect of surrounding rocks. Cao and Su [26] developed a numerical manifold method for simulating anchorage bolts, thereby improving the flexibility and accuracy of the numerical manifold method for simulating anchorage devices. Li and Cheng [27] proposed meshless numerical manifold methods for simulating crack propagation by establishing trial functions using unit partition and finite cover techniques [28–30], and Gao et al. [31] developed an improved meshless numerical manifold method with higher computational efficiency by incorporating a complex variable moving leastsquare method.
Discrete element methods (DEMs) are effective numerical computation methods, especially for solving noncontinuum problems. By the type of basic unit simulated, DEMs can be classified into two types: particle DEMs, which use discs or spheres as the basic unit, and block DEMs, which use blocky masses as the basic unit. Early DEMs include the DEM proposed by Peter Cundall and the discontinuous deformation analysis (DDA) proposed by Shi [32]. The DDA simulates a structure by segmenting it into blocks, accounts for the interaction between and the deformation of the blocks, and then solves the problem using the principle of minimum potential energy. The DDA is appropriate for solving large deformation and failure problems in geotechnical engineering. Aside from the DEM and DDA, the UDEC and PFC3D software packages developed by the Itasca engineering and software company have been used widely in geotechnical engineering [33].
Research of meshless methods dates back to the 1970s. Gingold and Monaghan and Lucy [34, 35] proposed smoothed particle hydrodynamics (SPH), which was first applied in astrophysics. A meshless method establishes shape functions without gridding; instead, it constructs and solves approximate functions using discrete points. Many meshless methods have been developed using different approximation functions [36–41] or different equivalent differential equations (e.g., Galerkin, least square, and collocation point methods) [42–47].
However, numerical simulation methods for solving largedeformation problems in geotechnical engineering are still in the testandtrial stage and are not ready for wide applications. Thus, this study proposes a numerical simulation method for solving large deformation and collapse problems in geotechnical engineering. Particularly, the method constructs barycentric interpolation trial functions based on the idea of manifold cover, realizes double discretization of Eulerian background meshes and moving material points, uses material particles as mesh integration points, and represents the macroscopic deformation behavior of objects using the motion of particles between meshes. Thus, this method is not affected by mesh distortion and allows unlimited computation in the predefined domain. The method was applied to a typical side slope. The resultof slip surface of the slope was then compared with that obtained using the Bishop method of slices. The method proposed in this study was then applied to simulate the large deformation (collapse) of the side slope, thereby confirming its validity and feasibility.
2. Barycentric InterpolationBased Manifold Cover FEMLIP
2.1. Motion Description Methods in Continuum Mechanics
In continuum mechanics, there are two classical methods for the description of motion: Lagrangian description and Eulerian description. These two methods are also commonly used for describing the motion of finite elements.
In the reference coordinate system shown in Figure 1, the initial configuration of the continuum at the initial moment, , is determined by coordinate X. The configuration at time point t, or the present configuration , is determined by coordinate x. The relationship between the initial and present configurations can be expressed as the equation . The Lagrangian description describes the motion of the continuum as a function of the initial coordinate, X, whereas the Eulerian description describes the motion as a function of the present coordinate, (initial coordinate and time).
In the Lagrangian description, mesh nodes (together with the material) are fixed to and move with material points. The Lagrangian description, focusing on material points and using the coordinate of material points for describing their motion, is mainly used for the stressstrain analysis of solid structures. When using the Lagrangian description for object analysis, the change in the object’s shape is completely consistent with that in the discrete meshes. This characteristic allows an accurate representation of the object’s motion and boundaries. However, when used for solving largedeformation problems, the Lagrangian description very often suffers mesh distortioncaused computation termination.
In the Eulerian description, mesh nodes and materials are independent from each other, and meshes are usually fixed [48]. The Eulerian description focuses on the field, with all physical quantities of material points represented as functions of spatial coordinate and time. The spatial coordinate is a purely mathematical, independent variable, and the motion of the material point at a fixed location is described temporally. When used for solving largedeformation problems, the Eulerian description does not have the problem of mesh distortion because the meshes do not change with the movement of material points. Thus, it is appropriate for solving fluid problems, but it cannot capture object boundaries because it focuses only on a certain material point.
The arbitrary Lagrangian–Eulerian (ALE) description has the advantages of both Lagrangian and Eulerian descriptions. It uses the Lagrangian description for representing the boundaries of mass and the Eulerian description for representing the interior of objects [49]. In the ALE, the meshes exist independently of the material points. However, they can be adjusted appropriately during problemsolving according to defined parameters, which are not completely the same as the Eulerian description. Even so, when used for solving complex largedeformation problems, the ALE still cannot avoid the mesh distortioncaused problem of computation termination.
2.2. Geometric Description of the Lagrangian Integration Point Method
The FEMLIP uses fixed Eulerian background meshes and discretizes materials using particles, with the properties of materials stored in the particles. It originated from the particleincell (PIC) method used in plasma research. During computation, the background meshes are fixed, and particles are used as the principal medium for transmitting the physical information of materials. Figure 2 illustrates the dual discretization of integration units and material particles and the motion of particles between Eulerian background meshes. Because of the discretization of both mesh nodes and material points, this method has both the capability of the Eulerianmesh FEM for computing large deformations and the capability of the Lagrangianmesh FEM for tracking internal variables.
(a)
(b)
(c)
(d)
2.3. Geometric Description of the Manifold Cover
The manifold cover system is a concept originating from numerical manifold methods. A manifold cover system consists of mathematical and physical manifolds. The problemsolving domain is defined using physical manifolds, and mathematical manifolds are used to define the integration domain and construct trial functions. The accuracy of the computational solutions usually depends on the density of mathematical manifolds. The FEM is a particular form of the manifold cover method, that is, a manifold cover method with the mathematical and physical manifolds completely overlapping. The displacement functions of the cover usually use constant functions, and the weight functions usually use linear functions. The meshless methods developed in recent years treat the effect of nodes or particles as mathematical covers and construct approximation functions through mathematical approximation. Figures 3–5 illustrate the cover system of the three methods. Here, thick lines represent physical covers, and thin lines represent mathematical covers. In the FEM, each node has its domain of influence, as shown in Figure 4. The influence domain of a node consists of all of its adjacent units (those highlighted with red lines in the figure). For planar rectangular meshes, the domain of influence of each node covers another eight nodes, and each node is covered by the influence domains of another eight nodes.
3. Barycentric InterpolationBased Lagrangian Integration Point Method
3.1. Construction of Barycentric Interpolation Trial Functions
The subdivision of an arbitrary polygon can be expressed as follows: . Each of the resulting polygonal unit can be illustrated as in Figure 6. Assume a polygonal unit with n sides, . Subdivide the polygon into n triangles by connecting the vertexes of the polygon with an arbitrary point inside the polygonal element, . Designating the interior angle of the triangles as and with , where X is the coordinate of P and , the following shape function can be constructed for the polygonal unit [50]:where is the weight function for node .
In finite element computation, to define boundary conditions of the model and represent the rigid object displacement mode and linear displacement field accurately, the shape function, , for a polygonal computational domain, , should satisfy the following three basic conditions [45]:(1)Nonnegativity and interpolation: where is the Kronecker delta.(2)Unit partition:(3)Linear completeness:
Assume a polygonal domain, , and a unit circle with a point in the polygon, P, as its center. The circle and line segment PP_{i} then intersect at H_{i}. A polygon circumscribing the circle and crossing H_{i} can then be derived, as shown in Figure 7.
According to the Gauss divergence theorem, a polygonbounded area Q_{1}, Q_{2}, ... Q_{n} satisfies the following equation:where is the domain boundary and n is the outward normal vector of the boundary unit.
When A = 1, (5) can be rewritten as follows:
The unit outward normal vector of a polygon Q_{1}, Q_{2}, ... Q_{n} can be expressed as follows: .
The boundary of the polygon can be integrated as follows:
Equation (7) can be rewritten through transposition as follows:
Thus, the coordinate of point P can be expressed as follows:
The weight function is defined as the ratio of the side length, Q_{i}, Q_{i−1}, of the polygon Q_{1}, Q_{2}, ... Q_{n} to the length of PP_{i}:
The following equation can be derived from Figure 7:
Substituting (11)into (10), the following equation can then be derived:
The barycentric interpolation function for the polygonal unit can then be derived from (1).
3.2. Basic Equations for the Barycentric InterpolationBased Lagrangian Integration Point Method
3.2.1. Balance Equation
Assume a problem in a planar domain. Its boundary, Г, can then be expressed as follows: , where is the velocity boundary and is the stress boundary. The control equations and boundary conditions for the planar domain can then be expressed as follows:(1)Balance equation:(2)Velocity and stress boundaries:where is the stress tensor, is the differential operator, is the volumetric force vector, is the face force vector, is the outward unit normal vector of boundary , and is the velocity.
The above balance equation is essentially a Navier–Stokes equation with inertia not considered:where is the density of the material, is the pressure, is the viscosity, is the deviation part of the strain rate, and is the volumetric force.
3.2.2. Geometric Equations
According to the unit trial functions, the relationship between the coordinates of any point in the unit and the node in the local coordinate system can be expressed as follows:
Similarly, the relationship between the displacements (velocities) of any point in the unit and the node can be expressed as follows:
The relationship between the strain and displacement of a unit can be expressed as follows:where L is the differential operator and , , and . Deriving both sides of (18) with respect to time, the relationship between the strain rate and the nodal velocity can then be obtained as follows:
3.2.3. Constitutive Equations
The constitutive relationship for representing the stressstrain behavior of objects of viscous constitution can be expressed as follows:where is the viscosity matrix, which can be expressed as follows:where η is the shear viscosity and ξ is the second viscosity parameter. These two parameters are similar with the Lamé elastic constants in the theory of elasticity, and the volume viscosity . For viscoelasticplastic constitutive models, the shear and volume viscosities are replaced with equivalent shear and volume viscosities, the computation methods for which will be detailed in the next section.
3.3. Solving Control Equations
According to the principle of virtual work, the weak form of integration of the discretization equation can be expressed as follows:where is the virtual velocity and is the virtual strain rate.
According to (19) and (20), (22) can be rewritten as follows:
The relationship between the velocity of any point in the unit and the nodal velocity can be expressed as the following trial function:
Thus, the following equation can be derived:
According to the arbitrariness of virtual velocity, the following equation can be derived:where and .
The mesh node velocity was resolved and then substituted into (19) to obtain the unit strain rate.
Considering the temporal variation of large deformations, a viscoelasticplastic constitutive model was adopted, which consisted of viscous, elastic, and plastic components connected in series, as shown in Figure 8.
The major parameters of the constitutive model include viscous shear modulus, viscous bulk modulus, elastic shear modulus, elastic bulk modulus, cohesion, and angle of friction.
For viscoelasticplastic models, the total strain rate can be expressed as the sum of viscous, elastic, and plastic strain rates:
Dividing the strain rate into two components, partial and spherical strain rates, the following equation can be derived:where is the trace of the strain rate matrix.
Stress rate is affected by the rigid object rotation. The component of the stress rate that is not affected by the rigid object rotation is referred to as the Jaumann stress rate, , which is obtained by deducting the effect of the rigid object rotation. Thus, the relationship between stress and strain rate can be expressed as follows:where is the Jaumann stress rate and satisfies the following equation: .
Designating and and defining the equivalent shear viscosity and the equivalent bulk viscosity , the first line of (29) can be rewritten as follows:
Thus, the following equation can be derived:
Designating as and as , the following equation can then be derived:
Designating , (32) can then be rewritten as follows:
Similarly, the following equation can be derived:
According to (15), with the inertial force not considered, the following equation can be derived:where f consists of three components: (1) external force, (2) the force resulting from the previous time step, and (3) the force resulting from the plastic deformation. The three components are designated as , , and, respectively, and satisfy the following equation:
3.4. Numerical Integration
Traditional FEMs usually select integration points according to predefined schemes so that the most accurate integration with respect to a minimum number of points can be realized (for example, Gauss integration). However,when mesh nodes move as the material deforms, traditional FEMs are prone to stop computation because of mesh distortion when used for computing large deformations. To resolve this problem, fixed Eulerian background meshes were used, and the particles moving inside the meshes were adopted to describe material deformation and object motion. With the moving particles as integration points, the instantaneous relationship between the material points and meshes can be established for each computational increment, thereby enabling accurate tracking of material deformation and object motion and ensuring mesh adaptability. Because the location of particles changes continuously, the integration weight needs to be updated to realize correct integration.
For n integration points, if appropriate weight coefficients and integration location are selected and the integrated function is a polynomial raised to a maximum power of 2_{n−1}, the following equation stands:where ψ is the given function, is the unit volume, P is the number of integration points (particles), ω_{n} is the weight function, and x_{n} is the location of integration points.
For onedimensional problems (linear units), the integrated function can be expressed as the following polynomial (integrated in the range of −1 to +1):
Integrating (38) in the range of [−1, 1], and then n + 1 equalities corresponding to ω_{n} can be obtained according to the following equation:
In principle, the value of ω_{n} can be estimated by defining an appropriate number of constraints. However, in real operations, this requires a large amount of computation.,Simultaneously, ωn calculated maybe negative values, and reasonable convergence cannot be ensured. For real applications, the relationship between the value of ωn of a particle and the mass or volume of the material can be established only when the value is positive.
For viscous flows using bilinear units, if only steady and linear constraints are used, then the optimum convergence rate under the unlimited mesh constraint can be obtained. The steady constraint is defined as follows:
The bilinear constraint is defined as follows:
The following definitions are then given:
The value of ω_{n} can then be obtained through the following iteration:
The convergence condition for the iteration is defined as follows:
The right side of equation (44) is the convergence accuracy that can be defined case by case.
After the mesh node velocity is obtained, the constantly changing location of particles can be updated as follows:
4. An Example Application
The stability of the side slope presented by Dawson et al., which has been cited by numerous studies, is analyzed using the method proposed above. The slip surface of the side slope obtained using the method was then compared with the Bishop method of slices in the literature. On this basis, the large deformation (collapse) of the side slope was simulated using the method proposed in this study, thereby confirming its accuracy and feasibility.
Figure 9 illustrates the computational model. The model adopted normal constraints, except for the upper surface. Table 1 shows the physicalmechanical parameters of the model. Figure 10 shows the shear strain rate yielded by the method proposed in this study. The potential slip surface of the side slope was then identified based on the strain rate and was then compared with that obtained using the Bishop method of slices, as shown in Figure 11.

Figure 11 shows that the slip surfaces yielded by two methods are basically consistent, thereby confirming that the method proposed in this study is valid and feasible for simulating large deformations in engineering. The evolution of the large deformation (collapse) of the slope was further observed, with the viscosity of the soil considered, as shown in Figure 12.
(a)
(b)
(c)
(d)
5. Conclusions
This study proposed a barycentric interpolation manifold method with Lagrangian integration points, a continuous numerical simulation method for simulating the discontinuous failure process of large deformations, by constructing barycentric interpolation trial functions based on the idea of the manifold cover, adopting double discretization of Eulerian background meshes and moving material points, and using material particles as mesh integration points. Our conclusions can be summarized as follows:(1)The method adopts Eulerian background integration meshes and represents the macroscopic deformation behavior of objects using the motion of particles between meshes, thereby avoiding the computation termination caused by mesh distortion that usually occurs with large deformation computation. Moreover, the method allows unlimited computation in the predefined domain and is suitable for simulating large deformations and collapses in geotechnical engineering.(2)The method proposed in this study was applied to a typical side slope. The slip surface yielded by the method was basically consistent with that yielded by a theoretical analysis method, thereby confirming the validity of the method proposed in this study. The method was then applied to simulate the collapse process of the side slope, confirming the feasibility of the method for computing large deformations [51].
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was financially supported by the National Natural Science Foundation of China (nos. 51879150 and 51809115), Qilu Construction Projects of Science and Technology in 2016 (no. 2016B20), and the Shandong Provincial Natural Science Foundation, China (nos. ZR2019QEE003 and 2018GGX104024). These financial supports are gratefully acknowledged.
References
 Y. Jiang, Y. Li, T. Li, and L. Wang, “Study of the classified system of types and mechanism of great distortion in tunnel and underground engineering,” Journal of Geological Hazards and Environment Preservation, vol. 15, no. 4, pp. 46–51, 2004. View at: Google Scholar
 B. Kou, L. Zhang, and C. Han, “Control study on large deformation of loess tunnel in inner Mongolia,” Subgrade Engineering, vol. 4, pp. 162–166, 2017. View at: Google Scholar
 Z. Liu, M. Wang, and J. Fang, “Insitu study of reinforcement systems for tunnels under high geostress and large deformation,” China Civil Engineering Journal, vol. 43, no. 5, pp. 111–116, 2010. View at: Google Scholar
 H. Wang and Z. Zhong, “Research on the large deformation of ultrashallow and largebias separated tunnel in loess Plateau and its treatment,” Chinese Journal of Underground Space and Engineering, vol. 8, no. 2, pp. 1841–1845, 2012. View at: Google Scholar
 C. Xu and B. Hua, Finite Element Theory, Method and Program of Solid Mechanics, Water Resources and Electric Power Press, Beijing, China, 1983.
 E. E. Alonso, A. Gens, and C. H. Delahaye, “Influence of rainfall on the deformation and stability of a slope in overconsolidated clays: a case study,” Hydrogeology Journal, vol. 11, no. 1, pp. 174–192, 2003. View at: Google Scholar
 M. Cai, Rock Mechanics and Engineering, Science Press, Beijing, China, 2013.
 Y. Nakata, D. Liu, M. Hyodo, N. Yoshimoto, and Y. Kato, Unsaturated Soil: Theoretical and Numerical Advances in Unsaturated Soil Mechanics, Newcastle, Australia, 2010.
 Y. Chen and D. Xu, FLAC/FLAC3D Foundation and Engineering Examples, China Water & Power Press, Beijing, China, 2009.
 S. Sun, Application of FLAC3D in Geotechnical Engineering, China Water & Power Press, Beijing, China, 2011.
 B. Zhu, The Finite Element Method Theory and Applications, China Water & Power Press, Beijing, China, 1998.
 L. Moresi, F. Dufour, and H.B. Mühlhaus, “Mantle convection modeling with viscoelastic/brittle lithosphere: numerical methodology and plate tectonic modeling,” Pure and Applied Geophysics, vol. 159, pp. 2335–2356, 2002. View at: Google Scholar
 L. Moresi, F. Dufour, and H.B. Mühlhaus, “A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials,” Journal of Computational Physics, vol. 184, pp. 476–497, 2003. View at: Google Scholar
 F. H. Harlow, “The particleincell computing method for fluid dynamics,” Computer Methods in Physics, vol. 3, pp. 319–343, 1964. View at: Google Scholar
 F. Bai, Z. Chen, and W. Bai, “Numerical simulation on the melting process of the mantle plumelithosphere interaction,” Chinese Journal of Geophysics, vol. 61, no. 4, pp. 1341–1351, 2018. View at: Google Scholar
 W. Meng, Z. Chen, and W. Bai, “Numerical simulation on process of the plumelithosphere interaction,” Chinese Journal of Geophysics, vol. 58, no. 2, pp. 495–503, 2015. View at: Google Scholar
 N. Prime, F. Dufour, and F. Darve, “Unified model for geomaterial solid/fluid states and the transition in between,” Journal of Engineering Mechanics, vol. 140, no. 6, pp. 682–694, 2013. View at: Google Scholar
 N. Prime, F. Dufour, and F. Darve, “Solidfluid transition modelling in geomaterials and application to a mudflow interacting with an obstacle,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 13, pp. 1341–1361, 2014. View at: Google Scholar
 Z. H. Li, F. Dufour, and F. Darve, “Hydroelastoplastic modelling with a solid/fluid transition,” Computers and Geotechnics, vol. 75, pp. 69–79, 2016. View at: Google Scholar
 Z. H. Li, F. Dufour, and F. Darve, “Modelling rainfallinduced mudflows using FEMLIP and a unified hydroelastoplastic model with solidfluid transition,” European Journal of Environmental & Civil Engineering, vol. 22, no. 4, pp. 1–31, 2016. View at: Google Scholar
 G.H. Shi, Manifold Method of Material Analysis, Applied Mathmatics and Computing, Minneapolish, MI, USA, 1992.
 M. Lu, HighOrder Manifold Method with Simplex Integration, Analysis of Discontinuous deformation, Wuhan, China, 2002.
 S. Wang, W. Zhu, S. Li, and S. Chen, “Numerical manifold method of elastoplastic analysis for rockmass,” Chinese Journal of Rock Mechanics and Engineering, vol. 21, no. 6, pp. 900–904, 2002. View at: Google Scholar
 X. Zhang, K. Song, and M. Lu, “Weighted residual method with compactly supported trial functions,” Acta Mechanica Sinica, vol. 35, no. 1, pp. 43–49, 2003. View at: Google Scholar
 Z. Dong, A. Wu, and X. Ding, “Rockbolts simulation by numerical manifold method and its preliminary application,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 20, pp. 3754–3760, 2005. View at: Google Scholar
 W. Cao and B. Su, “Simulation of numerical manifold method for supports by reinforced bolt and its application,” Chinese Journal of Geotechnical Engineering, vol. 23, no. 5, pp. 581–583, 2001. View at: Google Scholar
 S. Li and Y. Cheng, “Meshless numerical manifold method based on unit partition,” Chinese Journal of Theoretical and Applied Mechanics, vol. 36, no. 4, pp. 496–500, 2004. View at: Google Scholar
 Q. Jiang, S. Deng, and C. Zhou, “Study of threedimensional highorder numerical manifold method,” Rock and Soil Mechanics, vol. 27, no. 9, pp. 1472–1474, 2006. View at: Google Scholar
 S. Li, Y. Cheng, and S. Li, “Enriched meshless manifold method,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 12, pp. 2065–2073, 2005. View at: Google Scholar
 S. Li and Z. Wang, Meshless Barycentric Interpolation Collocation Method with High Precision, Science Press, Beijing, China, 2012.
 H. Gao, Y. Cheng, and H. Jiang, “A meshless manifold method with complex variables for elasticity,” Chinese Journal of Applied Mechanics, vol. 27, no. 1, pp. 15–19, 2010. View at: Google Scholar
 G.h. Shi, Discontinuous Deformation Analysis: A New Numerical Model for the Statics and Dynamics of Block Systems, University of California, Berkeley, CA, USA, 1988.
 I. Itasca, PFC3DParticle Flow Code in 3 Dimensions, Science Press, Minneapolis, MI, USA, 2003.
 R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamicsTheory and application to nonspherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, no. 3, pp. 375–389, 1977. View at: Google Scholar
 L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” The Astronomical Journal, vol. 8, no. 12, pp. 1013–1024, 1977. View at: Google Scholar
 I. Babuška and J. M. Melenk, “The partition of unity method,” International Journal for Numerical Methods in Engineering, vol. 40, no. 4, pp. 727–758, 1997. View at: Google Scholar
 I. Babuška and Z. Zhang, “The partition of unity method for the elastically supported beam,” Computer Methods in Applied Mechanics & Engineering, vol. 152, no. 1–2, pp. 1–18, 1998. View at: Google Scholar
 P. Lancaster and K. kauskas, “Surfaces generated by moving least square methods,” Mathematics of Computation, vol. 37, no. 155, pp. 141–158, 1981. View at: Google Scholar
 F. B. Liu and Y. M. Cheng, “The improved elementFree Galerkin method based on the nonsingular weight functions for inhomogeneous swelling of polymer gels,” International Journal of Applied Mechanics, vol. 10, no. 4, 2018. View at: Google Scholar
 F. B. Liu and Y. M. Cheng, “The improved elementfree Galerkin method based on the nonsingular weight functions for elastic large deformation problems,” International Journal of Computational Materials Science and Engineering, vol. 7, 2018. View at: Google Scholar
 F. B. Liu, Q. Wu, and Y. M. Cheng, “A meshless method based on the nonsingular weight functions for elastoplastic large deformation problems,” International Journal of Applied Mechanics, vol. 11, no. 1, 2019. View at: Google Scholar
 Y. J. Deng, X. Q. He, L. G. Sun, S. H. Yi, and Y. Dai, “An improved interpolating complex variable element free Galerkin method for the pattern transformation of hydrogel,” Engineering Analysis with Boundary Elements, vol. 113, pp. 99–109, 2020. View at: Google Scholar
 Y. J. Deng and X. Q. He, “An improved interpolating complex variable meshless method for bending problem of Kirchhoff plates,” International Journal of Applied Mechanics, vol. 9, no. 6, 2017. View at: Google Scholar
 Y. J. Deng, C. Liu, M. J. Peng, and Y. M. Cheng, “The interpolating complex variable element free Galerkin method for temperature field problems,” International Journal of Applied Mechanics, vol. 7, no. 2, 2015. View at: Google Scholar
 A. Tabarraei and N. Sukumar, “Application of polygonal finite elements in linear elasticity,” International Journal of Computational Methods, vol. 3, no. 4, pp. 503–520, 2006. View at: Google Scholar
 X. Zhang, W. Hu, X. Pan, and W. Lu, “Meshless weighted leastsquare method,” Acta Mechanica Sinica, vol. 35, no. 4, pp. 425–431, 2003. View at: Google Scholar
 X. Zhang, X. Liu, K. Song, and M. Lu, “Leastsquare collocation meshless method,” International Journal for Numerical Methods in Engineering, vol. 51, no. 9, pp. 1089–1100, 2001. View at: Google Scholar
 Z. Ding, Large Strain Consolidation Theory and Finite Element Method of Eulerian Description, Chang’an University, Chang’an, China, 2002.
 X. Zhang, M. Lu, and J. Wang, “Research progress in arbitrary LagrangianEulerian method,” Chinese Journal of Computational Mechanics, vol. 14, no. 1, pp. 91–102, 1997. View at: Google Scholar
 Z. Wang, Study on Rational Element Method, Shanghai University, Shanghai, China, 2004.
 G. Zhang, “Structure analysis by secondorder manifold method,” Journal of China Institute of Water Resources and Hydropower Research, vol. 1, no. 1, pp. 63–69, 2003. View at: Google Scholar
Copyright
Copyright © 2020 Hui Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.