Study on the Manifold Cover Lagrangian Integral Point Method Based on Barycentric Interpolation

Qilu Transportation Development Group Co., Ltd., Yinfeng Fortune Plaza D Block, No. 1 Lg’ao West Road, Jinan, Shandong Province, China College of Civil Engineering, Taiyuan University of Technology, No. 79 Yingze West Main Street, Taiyuan, Shanxi Province, China School of QiluTransportation, Shandong University, No. 12550 East Erhuan Road, Jinan, Shandong Province, China School of Civil Engineering and Architecture, University of Jinan, No. 336, West Road of Nan Xinzhuang, Jinan 250022, Shandong Province, China


Introduction
e 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 high-geostress soft rocks [2][3][4].
us, 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. e finite element method (FEM) and the finite difference method (FDM) are now the most common methods in geotechnical engineering for analyzing continuous deformations. e 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 side-slope collapses, many researchers have conducted in-depth investigation of the constitutive model of geomaterials from the perspective of their solid-liquid transformation [7,8], but the problem of computation termination caused by meshdistortion remains unresolved. e 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]. e 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]. e FEMLIP originates from the FEM and the particle-in-cell method [14]. Bai et al. and Meng et al. [15,16] used the FEMLIP to simulate mantle plume-lithosphere interactions. e FEMLIP has also been used to simulate mudflows [17,18]. Li et al. [19,20] simulated mudflows using the FEMLIP and hydro-elastoplastic models, which consider the rainfall-induced solid-liquid 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. e 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 high-order 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 second-order 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][29][30], and Gao et al. [31] developed an improved meshless numerical manifold method with higher computational efficiency by incorporating a complex variable moving least-square 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]. e 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. e 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].
However, numerical simulation methods for solving large-deformation problems in geotechnical engineering are still in the test-and-trial stage and are not ready for wide applications. us, 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. us, this method is not affected by mesh distortion and allows unlimited computation in the predefined domain. e method was applied to a typical side slope. e resultof slip surface of the slope was then compared with that obtained using the Bishop method of slices. e method proposed in this study was then applied to simulate the large deformation (collapse) of the side slope, thereby confirming its validity and feasibility.

Motion Description Methods in Continuum Mechanics.
In continuum mechanics, there are two classical methods for the description of motion: Lagrangian description and Eulerian description. ese 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, Ω 0 , is determined by coordinate X. e configuration at time point t, or the present configuration Ω, is determined by coordinate x.
e relationship between the initial and present configurations can be expressed as the equation x � ϕ(X, t). e 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, x � ϕ(X, t) (initial coordinate and time).
In the Lagrangian description, mesh nodes (together with the material) are fixed to and move with material points.
e Lagrangian description, focusing on material points and using the coordinate of material points for describing their motion, is mainly used for the stress-strain 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. is characteristic allows an accurate representation of the object's motion and boundaries. However, when used for solving large-deformation problems, the Lagrangian description very often suffers mesh distortion-caused computation termination.
In the Eulerian description, mesh nodes and materials are independent from each other, and meshes are usually fixed [48]. e Eulerian description focuses on the field, with all physical quantities of material points represented as functions of spatial coordinate and time. e 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 large-deformation problems, the Eulerian description does not have the problem of mesh distortion because the meshes do not change with the movement of material points. us, it is appropriate for solving fluid problems, but it cannot capture object boundaries because it focuses only on a certain material point. e 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 problem-solving according to defined parameters, which are not completely the same as the Eulerian description. Even so, when used for solving complex large-deformation problems, the ALE still cannot avoid the mesh distortion-caused problem of computation termination.

Geometric Description of the Lagrangian Integration Point
Method.
e FEMLIP uses fixed Eulerian background meshes and discretizes materials using particles, with the properties of materials stored in the particles. It originated from the particle-in-cell (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 Eulerian-mesh FEM for computing large deformations and the capability of the Lagrangianmesh FEM for tracking internal variables.

Geometric Description of the Manifold
Cover. e manifold cover system is a concept originating from numerical manifold methods. A manifold cover system consists of mathematical and physical manifolds. e problemsolving domain is defined using physical manifolds, and mathematical manifolds are used to define the integration domain and construct trial functions. e accuracy of the computational solutions usually depends on the density of mathematical manifolds. e FEM is a particular form of the manifold cover method, that is, a manifold cover method with the mathematical and physical manifolds completely overlapping. e displacement functions of the cover usually use constant functions, and the weight functions usually use linear functions. e 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. e 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.

Construction of Barycentric Interpolation Trial Functions.
e subdivision of an arbitrary polygon Ω can be expressed as follows: Ω � ∪ n e�1 Ω e . Each of the resulting polygonal unit Ω e can be illustrated as in Figure 6. Assume a polygonal unit with n sides, Ω e . Subdivide the polygon into n triangles by connecting the vertexes of the polygon (P 1 , P 2 , ..., P n ) with an arbitrary point inside the polygonal element, P ∈ Ω e . Designating the interior angle of the triangles as α i and with PP i � ‖X − X i ‖, where X is the coordinate of P and X � (x, y), the following shape function can be constructed for the polygonal unit [50]: where w i is the weight function for node P i .
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, N i (x) i � 1, 2, . . . , n, for a polygonal computational domain, Ω e , should satisfy the following three basic conditions [45]: (1) Nonnegativity and interpolation: where δ ij is the Kronecker delta.
(2) Unit partition: (3) Linear completeness: Assume a polygonal domain, Ω e , and a unit circle with a point in the polygon, P, as its center. e 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 zS is the domain boundary and n is the outward normal vector of the boundary unit. When A � 1, (5) can be rewritten as follows: e unit outward normal vector of a polygon Q 1 , Q 2 , ... Q n can be expressed as follows: n � (1/‖P − P i ‖)PP i ��→ . e boundary of the polygon can be integrated as follows: Equation (7) can be rewritten through transposition as follows: us, the coordinate of point P can be expressed as follows: e 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 : e following equation can be derived from Figure 7: Substituting (11)into (10), the following equation can then be derived: e barycentric interpolation function for the polygonal unit can then be derived from (1).

Basic Equations for the Barycentric Interpolation-Based
Lagrangian Integration Point Method

Balance Equation.
Assume a problem in a planar domain. Its boundary, Γ, can then be expressed as follows: where Γ v is the velocity boundary and Γ s is the stress boundary. e 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, f v is the volumetric force vector, f s is the face force vector, n is the outward unit normal vector of boundary Γ s , and v is the velocity. e above balance equation is essentially a Navier-Stokes equation with inertia not considered: where ρ is the density of the material, p is the pressure, η is the viscosity, _ e is the deviation part of the strain rate, and f is the volumetric force.

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: x � Nx n . (16) Similarly, the relationship between the displacements (velocities) of any point in the unit and the node can be expressed as follows: e relationship between the strain and displacement of a unit can be expressed as follows: where L is the differential operator 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:

Constitutive Equations.
e constitutive relationship for representing the stress-strain behavior of objects of viscous constitution can be expressed as follows: where D υ is the viscosity matrix, which can be expressed as follows: where η is the shear viscosity and ξ is the second viscosity parameter. ese two parameters are similar with the Lamé elastic constants in the theory of elasticity, and the volume viscosity K v � ξ + (2/3)η. For viscoelastic-plastic 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.

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: e relationship between the velocity of any point in the unit and the nodal velocity can be expressed as the following trial function: υ � Nυ n . (24) us, the following equation can be derived: According to the arbitrariness of virtual velocity, the following equation can be derived: where K � Ω Β eT D v BdΩ and f n � Ω N T f υ dV + Γ s N T f s dS. e mesh node velocity was resolved and then substituted into (19) to obtain the unit strain rate.
Considering the temporal variation of large deformations, a viscoelastic-plastic constitutive model was adopted, which consisted of viscous, elastic, and plastic components connected in series, as shown in Figure 8. e 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 viscoelastic-plastic 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 tr(_ ε) is the trace of the strain rate matrix. Stress rate is affected by the rigid object rotation. e 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. us, the relationship between stress and strain rate can be expressed as follows: where s ⌣ is the Jaumann stress rate and satisfies the following equation: s ⌣t+Δte � (s t+Δt e − s t )/Δt e − w t s t + s t w t .
Designating α � (η/μ) and β � (K υ /K e ) and defining the equivalent shear viscosity η eff � η(Δt e /Δt e + α) and the equivalent bulk viscosity (K υ ) eff � K υ (Δt e /Δt e + β), the first line of (29) can be rewritten as follows: 6 Mathematical Problems in Engineering us, the following equation can be derived: Designating _ e as _ e υ + _ e e and tr(_ ε) as tr(_ ε υ ) + tr(_ ε e ), the following equation can then be derived: Designating η eff � η(Δt e /Δt e + α), (32) can then be rewritten as follows: s t+Δt e � η eff 2_ e t+Δt e + s t μΔt e + w t s t − s t w t μ . (33) Similarly, the following equation can be derived: According to (15), with the inertial force not considered, the following equation can be derived: (35) 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. e three components are designated as f 1 , f 2 , andf 3 , respectively, and satisfy the following equation:

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, Ω e 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 one-dimensional 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 e steady constraint is defined as follows: e bilinear constraint is defined as follows: e following definitions are then given: e value of ω n can then be obtained through the following iteration: e convergence condition for the iteration is defined as follows: e 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:

An Example Application
e stability of the side slope presented by Dawson et al., which has been cited by numerous studies, is analyzed using the method proposed above. e 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. e model adopted normal constraints, except for the upper surface. Table 1 shows the physical-mechanical parameters of the model. Figure 10 shows the shear strain rate yielded by the   Mathematical Problems in Engineering method proposed in this study. e 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. e evolution of the large deformation (collapse) of the slope was further observed, with the viscosity of the soil considered, as shown in Figure 12.

Conclusions
is 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) e 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) e method proposed in this study was applied to a typical side slope. e 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. e 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
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.