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

Authors: Hui Sun; Qin Yan; Bing Han; Shuchen Li; Xianda Feng
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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 large-deformation 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.


---

## Body

## 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 high-geostress 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 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. 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 particle-in-cell method [14]. Bai et al. and Meng et al. [15, 16] used the FEMLIP to simulate mantle plume-lithosphere interactions. The 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. 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 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–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]. 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 large-deformation problems in geotechnical engineering are still in the test-and-trial 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 Interpolation-Based 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 Figure1, the initial configuration of the continuum at the initial moment, Ω0, 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 x=ϕX,t. 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, x=ϕX,t (initial coordinate and time).Figure 1
Reference coordinate system.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 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. This 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]. 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 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. 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 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.


### 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 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. Figure2 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 Lagrangian-mesh FEM for tracking internal variables.Figure 2
Dual discretization of Eulerian meshes and material particles.
(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 problem-solving 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. Figures3–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.Figure 3
The cover system of the numerical manifold method.Figure 4
The cover system of the FEM.Figure 5
The cover system of the meshless method.


## 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 Figure1, the initial configuration of the continuum at the initial moment, Ω0, 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 x=ϕX,t. 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, x=ϕX,t (initial coordinate and time).Figure 1
Reference coordinate system.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 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. This 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]. 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 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. 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 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.


## 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 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. Figure2 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 Lagrangian-mesh FEM for tracking internal variables.Figure 2
Dual discretization of Eulerian meshes and material particles.
(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 problem-solving 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. Figures3–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.Figure 3
The cover system of the numerical manifold method.Figure 4
The cover system of the FEM.Figure 5
The cover system of the meshless method.


## 3. Barycentric Interpolation-Based Lagrangian Integration Point Method

### 3.1. Construction of Barycentric Interpolation Trial Functions

The subdivision of an arbitrary polygonΩ can be expressed as follows: Ω=∪e=1nΩ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 P1,P2,...,Pn with an arbitrary point inside the polygonal element, P∈Ωe. Designating the interior angle of the triangles as αi and with PPi=X−Xi, where X is the coordinate of P and X=x,y, the following shape function can be constructed for the polygonal unit [50]:(1)Nix=wi∑j=1nwj,i=1,2,…,n,where wi is the weight function for node Pi.Figure 6
Polygonal element.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,Nixi=1,2,…,n, for a polygonal computational domain, Ωe, should satisfy the following three basic conditions [45]:(1)
Nonnegativity and interpolation:(2)0≤NiP≤1,NiPj=δij,whereδij is the Kronecker delta.(2)
Unit partition:(3)∑i=1nNiP=1.(3)
Linear completeness:(4)P=∑i=1nNiPPi.Assume a polygonal domain,Ωe, and a unit circle with a point in the polygon, P, as its center. The circle and line segment PPi then intersect at Hi. A polygon circumscribing the circle and crossing Hi can then be derived, as shown in Figure 7.Figure 7
Geometric structure of the barycentric interpolation function.According to the Gauss divergence theorem, a polygon-bounded areaQ1, Q2, ... Qn satisfies the following equation:(5)∫V∇AdV=∫∂SAndS,where ∂S is the domain boundary and n is the outward normal vector of the boundary unit.WhenA = 1, (5) can be rewritten as follows:(6)∫∂VndS=0.The unit outward normal vector of a polygonQ1, Q2, ... Qn can be expressed as follows: n=1/P−PiPPi⟶.The boundary of the polygon can be integrated as follows:(7)∑i=1nQi−Qi−1P−PiPPi⟶=∑i=1nQi−Qi−1P−PiPi−P=0.Equation (7) can be rewritten through transposition as follows:(8)P∑i=1nQi−Qi−1P−Pi=∑i=1nQi−Qi−1P−PiPi.Thus, the coordinate of pointP can be expressed as follows:(9)P=∑i=1nQi−Qi−1/P−PiPi∑i=1nQi−Qi−1/P−Pi.The weight function is defined as the ratio of the side length,Qi, Qi−1, of the polygon Q1, Q2, ... Qn to the length of PPi:(10)wiP=Qi−Qi−1P−Pi.The following equation can be derived from Figure7:(11)Qi−Qi−1=tanαi−12+tanαi2.Substituting (11)into (10), the following equation can then be derived:(12)wiP=tanαi−1/2+tanαi/2P−Pi.The barycentric interpolation function for the polygonal unit can then be derived from (1).


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

#### 3.2.1. Balance Equation

Assume a problem in a planar domain. Its boundary, Г, can then be expressed as follows:Γ=Γv+Γs, where Γv is the velocity boundary and Γs is the stress boundary. The control equations and boundary conditions for the planar domain can then be expressed as follows:(1)
Balance equation:(13)∇σ+fv=0,in domainΩ.(2)
Velocity and stress boundaries:(14)σn=fs,on stress boundaryΓs,v=v¯,on velocity boundaryΓv,where σ is the stress tensor, ∇ is the differential operator, fv is the volumetric force vector, fs is the face force vector, n is the outward unit normal vector of boundary Γs, and v is the velocity.The above balance equation is essentially a Navier–Stokes equation with inertia not considered:(15)ρDvDt=−gradp+div2ηe˙+f,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.


#### 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:(16)x=Nxn.Similarly, the relationship between the displacements (velocities) of any point in the unit and the node can be expressed as follows:(17)u=Nδ,υ=Nυn.The relationship between the strain and displacement of a unit can be expressed as follows:(18)ε=εxεyγxy=Lu=LNδ=Bδ,where L is the differential operator and L=∂/∂x00∂/∂y∂/∂y∂/∂x, B=B1e…Bie…Bne, and Bie=∂Ni/∂x00∂Ni/∂y∂Ni/∂y∂Ni/∂x. 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:(19)ε˙=LNυn=Bυn.


#### 3.2.3. Constitutive Equations

The constitutive relationship for representing the stress-strain behavior of objects of viscous constitution can be expressed as follows:(20)σ=Dυε˙,where Dυ is the viscosity matrix, which can be expressed as follows:(21)Dυ=ξ+2ηξξ000ξξ+2ηξ000ξξξ+2η000000η000000η000000η,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 Kv=ξ+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.


### 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:(22)∫Ωε˙∗:σdV−∫Ωυ∗TfυdV−∫Γsυ∗TfsdS=0,where υ∗ is the virtual velocity and ε˙∗ is the virtual strain rate.According to (19) and (20), (22) can be rewritten as follows:(23)∫Ωυn∗TBTDυBυndV−∫Ωυ∗TfυdV−∫Γsυ∗TfsdS=0.The relationship between the velocity of any point in the unit and the nodal velocity can be expressed as the following trial function:(24)υ=Nυn.Thus, the following equation can be derived:(25)υn∗T∫ΩBTDυBdVυn−∫ΩNTfυdV−∫ΓsNTfsdS=0.According to the arbitrariness of virtual velocity, the following equation can be derived:(26)Kυn=fn,where K=∫ΩΒeTDvBdΩ and fn=∫ΩNTfυdV+∫ΓsNTfsdS.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 viscoelastic-plastic constitutive model was adopted, which consisted of viscous, elastic, and plastic components connected in series, as shown in Figure8.Figure 8
Viscoelastic-plastic constitutive model.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 viscoelastic-plastic models, the total strain rate can be expressed as the sum of viscous, elastic, and plastic strain rates:(27)ε˙tot=ε˙v+ε˙e+ε˙p.Dividing the strain rate into two components, partial and spherical strain rates, the following equation can be derived:(28)e˙tot=e˙υ+e˙p+e˙e,trε˙tot=trε˙v+trε˙p+trε˙e,where trε˙ 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:(29)s⌣2μ+s2η=e˙tot−e˙p=e˙υ+e˙e,p˙Ke+pKυ=trε˙tot−trε˙p=trε˙υ+trε˙e,where s⌣ is the Jaumann stress rate and satisfies the following equation: s⌣t+Δte=st+Δte−st/Δte−wtst+stwt.Designatingα=η/μ and β=Kυ/Ke and defining the equivalent shear viscosity ηeff=ηΔte/Δte+α and the equivalent bulk viscosity Kυeff=KυΔte/Δte+β, the first line of (29) can be rewritten as follows:(30)αst+Δte−stΔte−wtst+stwt+st+Δte=2ηe˙υ+e˙e.Thus, the following equation can be derived:(31)α+ΔteΔtest+Δte−αΔtest−αwtst−stwt=2ηe˙υ+e˙e.Designatinge˙ as e˙υ+e˙e and trε˙ as trε˙υ+trε˙e, the following equation can then be derived:(32)st+Δte=2ηΔteα+Δtee˙+αα+Δtest+αΔteα+Δtewtst−stwt.Designatingηeff=ηΔte/Δte+α, (32) can then be rewritten as follows:(33)st+Δte=ηeff2e˙t+Δte+stμΔte+wtst−stwtμ.Similarly, the following equation can be derived:(34)pt+Δte=Kυefftrε˙t+Δte+ptKeΔte.According to (15), with the inertial force not considered, the following equation can be derived:(35)2ηeffdive˙t+Δte+Kveffgradtrε˙t+Δte=f,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 f1, f2, andf3, respectively, and satisfy the following equation:(36)f1=fextt+Δte,f2=−ηeffdivstμΔte+divwtst−stwtμ−KυeffgradptKeΔte,f3=2ηeffdive˙pt+Δte+Kveffgradtrε˙pt+Δte.


### 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.Forn integration points, if appropriate weight coefficients and integration location are selected and the integrated function is a polynomial raised to a maximum power of 2n−1, the following equation stands:(37)∫ΩeψdΩe≈∑n=1Pωnψxn,where ψ is the given function, Ωe is the unit volume, P is the number of integration points (particles), ωn is the weight function, and xn 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):(38)ψx=α0+α1x+α2x2+,…,+αmxm.Integrating (38) in the range of [−1, 1], and then n + 1 equalities corresponding to ωn can be obtained according to the following equation:(39)∑n=1Pωn=2,∑n=1Pωnxn=0,∑n=1Pωnxn2=23,∑n=1pωnxn3=0.…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:(40)M0=∑nωn.The bilinear constraint is defined as follows:(41)Mξ=∑nωnξn,Mζ=∑nωnζn,Mξζ=∑nωnξnζn.The following definitions are then given:(42)Lξ=∑nξn2,Lζ=∑nζn2,Lξζ=∑nξn2ζn2.The value ofωn can then be obtained through the following iteration:(43)ωni+1=1M0i+1ωni−MξiLξi−MζiLζi−MξζiLξζi.The convergence condition for the iteration is defined as follows:(44)∑j=1npωjiξj+∑j=1npωjiζj+∑j=1npωjiξjζj≤ε.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:(45)xpt+Δt=xpt+Δt∑nNnxpυn.


## 3.1. Construction of Barycentric Interpolation Trial Functions

The subdivision of an arbitrary polygonΩ can be expressed as follows: Ω=∪e=1nΩ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 P1,P2,...,Pn with an arbitrary point inside the polygonal element, P∈Ωe. Designating the interior angle of the triangles as αi and with PPi=X−Xi, where X is the coordinate of P and X=x,y, the following shape function can be constructed for the polygonal unit [50]:(1)Nix=wi∑j=1nwj,i=1,2,…,n,where wi is the weight function for node Pi.Figure 6
Polygonal element.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,Nixi=1,2,…,n, for a polygonal computational domain, Ωe, should satisfy the following three basic conditions [45]:(1)
Nonnegativity and interpolation:(2)0≤NiP≤1,NiPj=δij,whereδij is the Kronecker delta.(2)
Unit partition:(3)∑i=1nNiP=1.(3)
Linear completeness:(4)P=∑i=1nNiPPi.Assume a polygonal domain,Ωe, and a unit circle with a point in the polygon, P, as its center. The circle and line segment PPi then intersect at Hi. A polygon circumscribing the circle and crossing Hi can then be derived, as shown in Figure 7.Figure 7
Geometric structure of the barycentric interpolation function.According to the Gauss divergence theorem, a polygon-bounded areaQ1, Q2, ... Qn satisfies the following equation:(5)∫V∇AdV=∫∂SAndS,where ∂S is the domain boundary and n is the outward normal vector of the boundary unit.WhenA = 1, (5) can be rewritten as follows:(6)∫∂VndS=0.The unit outward normal vector of a polygonQ1, Q2, ... Qn can be expressed as follows: n=1/P−PiPPi⟶.The boundary of the polygon can be integrated as follows:(7)∑i=1nQi−Qi−1P−PiPPi⟶=∑i=1nQi−Qi−1P−PiPi−P=0.Equation (7) can be rewritten through transposition as follows:(8)P∑i=1nQi−Qi−1P−Pi=∑i=1nQi−Qi−1P−PiPi.Thus, the coordinate of pointP can be expressed as follows:(9)P=∑i=1nQi−Qi−1/P−PiPi∑i=1nQi−Qi−1/P−Pi.The weight function is defined as the ratio of the side length,Qi, Qi−1, of the polygon Q1, Q2, ... Qn to the length of PPi:(10)wiP=Qi−Qi−1P−Pi.The following equation can be derived from Figure7:(11)Qi−Qi−1=tanαi−12+tanαi2.Substituting (11)into (10), the following equation can then be derived:(12)wiP=tanαi−1/2+tanαi/2P−Pi.The barycentric interpolation function for the polygonal unit can then be derived from (1).


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

### 3.2.1. Balance Equation

Assume a problem in a planar domain. Its boundary, Г, can then be expressed as follows:Γ=Γv+Γs, where Γv is the velocity boundary and Γs is the stress boundary. The control equations and boundary conditions for the planar domain can then be expressed as follows:(1)
Balance equation:(13)∇σ+fv=0,in domainΩ.(2)
Velocity and stress boundaries:(14)σn=fs,on stress boundaryΓs,v=v¯,on velocity boundaryΓv,where σ is the stress tensor, ∇ is the differential operator, fv is the volumetric force vector, fs is the face force vector, n is the outward unit normal vector of boundary Γs, and v is the velocity.The above balance equation is essentially a Navier–Stokes equation with inertia not considered:(15)ρDvDt=−gradp+div2ηe˙+f,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.


### 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:(16)x=Nxn.Similarly, the relationship between the displacements (velocities) of any point in the unit and the node can be expressed as follows:(17)u=Nδ,υ=Nυn.The relationship between the strain and displacement of a unit can be expressed as follows:(18)ε=εxεyγxy=Lu=LNδ=Bδ,where L is the differential operator and L=∂/∂x00∂/∂y∂/∂y∂/∂x, B=B1e…Bie…Bne, and Bie=∂Ni/∂x00∂Ni/∂y∂Ni/∂y∂Ni/∂x. 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:(19)ε˙=LNυn=Bυn.


### 3.2.3. Constitutive Equations

The constitutive relationship for representing the stress-strain behavior of objects of viscous constitution can be expressed as follows:(20)σ=Dυε˙,where Dυ is the viscosity matrix, which can be expressed as follows:(21)Dυ=ξ+2ηξξ000ξξ+2ηξ000ξξξ+2η000000η000000η000000η,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 Kv=ξ+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.


## 3.2.1. Balance Equation

Assume a problem in a planar domain. Its boundary, Г, can then be expressed as follows:Γ=Γv+Γs, where Γv is the velocity boundary and Γs is the stress boundary. The control equations and boundary conditions for the planar domain can then be expressed as follows:(1)
Balance equation:(13)∇σ+fv=0,in domainΩ.(2)
Velocity and stress boundaries:(14)σn=fs,on stress boundaryΓs,v=v¯,on velocity boundaryΓv,where σ is the stress tensor, ∇ is the differential operator, fv is the volumetric force vector, fs is the face force vector, n is the outward unit normal vector of boundary Γs, and v is the velocity.The above balance equation is essentially a Navier–Stokes equation with inertia not considered:(15)ρDvDt=−gradp+div2ηe˙+f,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.


## 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:(16)x=Nxn.Similarly, the relationship between the displacements (velocities) of any point in the unit and the node can be expressed as follows:(17)u=Nδ,υ=Nυn.The relationship between the strain and displacement of a unit can be expressed as follows:(18)ε=εxεyγxy=Lu=LNδ=Bδ,where L is the differential operator and L=∂/∂x00∂/∂y∂/∂y∂/∂x, B=B1e…Bie…Bne, and Bie=∂Ni/∂x00∂Ni/∂y∂Ni/∂y∂Ni/∂x. 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:(19)ε˙=LNυn=Bυn.


## 3.2.3. Constitutive Equations

The constitutive relationship for representing the stress-strain behavior of objects of viscous constitution can be expressed as follows:(20)σ=Dυε˙,where Dυ is the viscosity matrix, which can be expressed as follows:(21)Dυ=ξ+2ηξξ000ξξ+2ηξ000ξξξ+2η000000η000000η000000η,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 Kv=ξ+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.


## 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:(22)∫Ωε˙∗:σdV−∫Ωυ∗TfυdV−∫Γsυ∗TfsdS=0,where υ∗ is the virtual velocity and ε˙∗ is the virtual strain rate.According to (19) and (20), (22) can be rewritten as follows:(23)∫Ωυn∗TBTDυBυndV−∫Ωυ∗TfυdV−∫Γsυ∗TfsdS=0.The relationship between the velocity of any point in the unit and the nodal velocity can be expressed as the following trial function:(24)υ=Nυn.Thus, the following equation can be derived:(25)υn∗T∫ΩBTDυBdVυn−∫ΩNTfυdV−∫ΓsNTfsdS=0.According to the arbitrariness of virtual velocity, the following equation can be derived:(26)Kυn=fn,where K=∫ΩΒeTDvBdΩ and fn=∫ΩNTfυdV+∫ΓsNTfsdS.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 viscoelastic-plastic constitutive model was adopted, which consisted of viscous, elastic, and plastic components connected in series, as shown in Figure8.Figure 8
Viscoelastic-plastic constitutive model.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 viscoelastic-plastic models, the total strain rate can be expressed as the sum of viscous, elastic, and plastic strain rates:(27)ε˙tot=ε˙v+ε˙e+ε˙p.Dividing the strain rate into two components, partial and spherical strain rates, the following equation can be derived:(28)e˙tot=e˙υ+e˙p+e˙e,trε˙tot=trε˙v+trε˙p+trε˙e,where trε˙ 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:(29)s⌣2μ+s2η=e˙tot−e˙p=e˙υ+e˙e,p˙Ke+pKυ=trε˙tot−trε˙p=trε˙υ+trε˙e,where s⌣ is the Jaumann stress rate and satisfies the following equation: s⌣t+Δte=st+Δte−st/Δte−wtst+stwt.Designatingα=η/μ and β=Kυ/Ke and defining the equivalent shear viscosity ηeff=ηΔte/Δte+α and the equivalent bulk viscosity Kυeff=KυΔte/Δte+β, the first line of (29) can be rewritten as follows:(30)αst+Δte−stΔte−wtst+stwt+st+Δte=2ηe˙υ+e˙e.Thus, the following equation can be derived:(31)α+ΔteΔtest+Δte−αΔtest−αwtst−stwt=2ηe˙υ+e˙e.Designatinge˙ as e˙υ+e˙e and trε˙ as trε˙υ+trε˙e, the following equation can then be derived:(32)st+Δte=2ηΔteα+Δtee˙+αα+Δtest+αΔteα+Δtewtst−stwt.Designatingηeff=ηΔte/Δte+α, (32) can then be rewritten as follows:(33)st+Δte=ηeff2e˙t+Δte+stμΔte+wtst−stwtμ.Similarly, the following equation can be derived:(34)pt+Δte=Kυefftrε˙t+Δte+ptKeΔte.According to (15), with the inertial force not considered, the following equation can be derived:(35)2ηeffdive˙t+Δte+Kveffgradtrε˙t+Δte=f,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 f1, f2, andf3, respectively, and satisfy the following equation:(36)f1=fextt+Δte,f2=−ηeffdivstμΔte+divwtst−stwtμ−KυeffgradptKeΔte,f3=2ηeffdive˙pt+Δte+Kveffgradtrε˙pt+Δte.


## 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.Forn integration points, if appropriate weight coefficients and integration location are selected and the integrated function is a polynomial raised to a maximum power of 2n−1, the following equation stands:(37)∫ΩeψdΩe≈∑n=1Pωnψxn,where ψ is the given function, Ωe is the unit volume, P is the number of integration points (particles), ωn is the weight function, and xn 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):(38)ψx=α0+α1x+α2x2+,…,+αmxm.Integrating (38) in the range of [−1, 1], and then n + 1 equalities corresponding to ωn can be obtained according to the following equation:(39)∑n=1Pωn=2,∑n=1Pωnxn=0,∑n=1Pωnxn2=23,∑n=1pωnxn3=0.…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:(40)M0=∑nωn.The bilinear constraint is defined as follows:(41)Mξ=∑nωnξn,Mζ=∑nωnζn,Mξζ=∑nωnξnζn.The following definitions are then given:(42)Lξ=∑nξn2,Lζ=∑nζn2,Lξζ=∑nξn2ζn2.The value ofωn can then be obtained through the following iteration:(43)ωni+1=1M0i+1ωni−MξiLξi−MζiLζi−MξζiLξζi.The convergence condition for the iteration is defined as follows:(44)∑j=1npωjiξj+∑j=1npωjiζj+∑j=1npωjiξjζj≤ε.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:(45)xpt+Δt=xpt+Δt∑nNnxpυn.


## 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.Figure9 illustrates the computational model. The 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 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 9
Computational model.Table 1
Physical-mechanical parameters of the slope.
Density (kg·m−3)Bulk modulus (MPa)Shear modulus (MPa)Cohesion (kPa)Angle of friction (°)Maximum tensile strength (MPa)2,000100300.35200.0045Figure 10
Partial strain rate yielded by the method proposed in this study.Figure 11
Potential slip surfaces computed using the two methods.Figure11 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.Figure 12
Slope collapse process. (a) Step 8. (b) Step 12. (c) Step 28. (d) Step 50.
(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].



---

*Source: 6743204-2020-08-29.xml*