A three-dimensional two-level gradient smoothing meshfree method for rainfall induced landslide simulations

Abstract

A three-dimensional two-level gradient smoothing meshfree method is presented for rainfall induced landslide simulations. The two-level gradient smoothing for meshfree shape function is elaborated in the three-dimensional Lagrangian setting with detailed implementation procedure. It is shown that due to the successive gradient smoothing operation without the requirement of derivative computation in the present formulation, the two-level smoothed gradient of meshfree shape function is capable of achieving a given influence domain more efficiently than the standard gradient of meshfree shape function. Subsequently, the two-level smoothed gradient of meshfree shape function is employed to discretize the weak form of coupled rainfall seepage and soil motion equations in a nodal integration format, as provides an efficient three-dimensional regularized meshfree formulation for large deformation rainfall induced landslide simulations. The exponential damage and pressure dependent plasticity relationships are utilized to describe the failure evolution in landslides. The plastic response of soil is characterized by the true effective stress measure, which is updated according to the rotationally neutralized objective integration algorithm. The effectiveness of the present three-dimensional two-level gradient smoothing meshfree method is demonstrated through numerical examples.

Keywords

meshfree method / landslide / rainfall / three-dimensional two-level gradient smoothing / nodal integration

Introduction

Computational simulation has been playing an increasingly important role on the modeling, prediction and mitigation of the rainfall induced slope failure, a typical and frequent natural disaster possibly with severe damage and impact to the society. The slope failure induced by rainfall encompasses the large deformation complex failure evolution under the circumstance of soil-water interaction and thus poses considerable difficulties for the conventional finite element analysis [¹]. On the other hand, meshfree methods based upon unstructured node approximation have experienced very rapid developments in the past more than two decades. According to the approximation on either strong or weak forms, meshfree methods can be classified into the strong form-based collocation type of methods and the weak form-based Galerkin type of methods. Smoothed particle hydrodynamics (SPH) method [²–⁴] is a typical approach among the many meshfree collocation methods, while the element free Galerkin method [⁵] and the reproducing kernel particle method [⁶] are two representatives of the widely used Galerkin type of meshfree methods. Meshfree methods have showed noticeable advantages on large deformation modeling, moving boundary simulation, and damage and failure analysis, etc. [⁷–¹⁷]. Comprehensive summary and classification of meshfree methods and their applications can be found in Refs. [¹⁸–²³].

As for the slope failure or landslide simulations in the context of meshfree collocation methods, Bui et al. [²⁴] performed the slope stability analysis and failure simulation with SPH. Pastor et al. [²⁵] carried out a landslide run-out analysis using a SPH depth-integrated model, and Hu et al. [²⁶] investigated the landslide run-out behavior through a three-dimensional (3D) SPH computation. A SPH simulation particularly focusing on the municipal solid waste landfills was also reported by Dai and Huang [²⁷]. In contrast to the collocation formulation, Galerkin meshfree methods are more preferred from the stability and accuracy points of view. Rabczuk and Areias [²⁸] proposed an interesting element free Galerkin analysis of geomaterial failure where the slip lines are modelled as slipped particles. The numerical manifold method, in combination with the vector sum method and the graph theory, was presented for slope stability assessment by Zhang et al. [²⁹] and Liu et al. [³⁰]. Nonetheless, the efficiency is a major concern for Galerkin meshfree methods such as the element free Galerkin method with moving least square (MLS) approximation [⁵] and the reproducing kernel particle method with reproducing kernel (RK) approximation [⁶]. The support overlapping and non-polynomial nature of MLS/RK shape functions necessitate higher order quadrature for the domain integration [³¹] and consequently the development of efficient schemes for meshfree domain integration has been an important topic of current interests [³²–³⁹]. A natural and efficient nodally integrated meshfree formulation suffers rank deficiency issue [⁴⁰]. Chen et al. [⁴¹,⁴²] proposed a stabilized conforming nodal integration (SCNI) algorithm which simultaneously achieve the efficiency as well as stability for Galerkin meshfree methods. This approach has been further developed by Wang et al. [¹³] for 3D simulation of large deformation failure evolution in soils. Kwok et al. [⁴³] presented a slope stability analysis and post-failure simulation by employing the semi-Lagrangian reproducing kernel particle method with non-conforming stabilized nodal integration [⁴⁴]. Moreover, the coupled finite element and material point method and the particle finite element method have also been adopted by Lian et al. [⁴⁵] and Zhang et al. [⁴⁶] for the landslide simulations, among others.

Another important aspect for the slope failure or landslide analysis is to properly deal with the material instability and related discretization sensitivity issue [⁴⁷]. In meshfree methods, the meshfree shape functions have an inherent non-local feature and their adjustable support size may serve as an intrinsic length scale and regularize the numerical solutions. This characteristic is systematically explored by Chen et al. [⁴⁸,⁴⁹] to develop a general reproducing kernel strain smoothing regularized framework for strain localization analysis, which is closely related to the gradient enhanced models [⁵⁰]. In order to improve the computational efficiency, Wang and Li [⁵¹] introduced the stabilized conforming nodal integration into the general strain smoothing operation and developed a two-level strain smoothing meshfree method for small deformation elastic damage analysis without the problem of discretization sensitivity. Subsequently, based upon the two-level strain or gradient smoothing technique a two-dimensional (2D) meshfree formulation was developed by Wang et al. [⁵²] for rainfall-induced slope failure analysis.

In this work, a 3D two-level gradient smoothing meshfree method for rainfall-induced landslide simulations is further developed. The two-level gradient smoothing operation is carried out for the Lagrangian meshfree shape functions referring to the initial configuration, which shows obvious advantages regarding the efficiency and kernel stability considerations [⁵,¹³,⁵³]. In particular, the smoothed gradient computation for meshfree shape function is discussed in detail. It is shown that through the two-level gradient smoothing operation, the resulting smoothed gradient of meshfree shape function turns out to have a relatively larger influence domain with less computational effort compared with the standard meshfree shape function gradient, as is beneficial for the discretization insensitive modeling of strain softening problems arising in landslides. Thereafter, the two-level smoothed gradient of meshfree shape function is employed to discretize the weak form of coupled soil motion and seepage equations. The damage and plastic responses of soil are described by the exponential damage evolution law [⁵⁴–⁵⁶] and the Drucker-Prager model with the true effective stress measure. The rotationally neutralized algorithm [⁵⁷] is employed for the objective incremental stress update. The effectiveness of the present 3D two-level gradient smoothing meshfree formulation for landslide simulations are illustrated through numerical examples.

The organization of this paper is as follows. In Section 2, the basic equations of soil motion and rainfall infiltration are outlined, which are followed by the damage and plasticity relationships. In Section 3, the Lagrangian meshfree approximation is briefly described. Then the two-level gradient smoothing of meshfree shape function and its numerical implementation are elaborated in the 3D setting. Subsequently, the discrete two-level gradient smoothing meshfree equations and stress update procedure are given for landslide simulations. In Section 4, numerical demonstration of the proposed method is given. Finally, conclusions are drawn in Section 5.

Basic equations

Kinematics and balance equations

In the landslide modeling, the soil-water mixture is considered. For a soil-water mixture, a material point \mathit{X} in the initial configuration {\mathrm{\Omega }}_{\mathrm{0}} with boundary {\mathrm{\Gamma }}_{\mathrm{0}} moves to the current position \mathit{x} after a motion defined by \mathit{x}=\mathit{\text{φ}}(\mathit{X},t), and the current configuration is denoted by \mathrm{\Omega } with boundary \mathrm{\Gamma }. Referring to the current configuration, the soil-water mixture density \rho can be expressed as [⁵⁸]

(1)\rho =(\mathrm{1}-n){{\rho }_{\mathrm{s}}}_{}+n{S_{\mathrm{r}}}_{}{{\rho }_{\mathrm{w}}}_{},

where n is the soil porosity, {S_{}}_{\mathrm{r}} is the degree of saturation, {{\rho }_{}}_{\mathrm{s}} and {{\rho }_{}}_{\mathrm{w}} are the densities of soil skeleton and water, respectively. Equation (1) also implies that the air phase is neglected herein.

According to the motion of \mathit{x}=\phi (\mathit{X},t), the displacement gradient \mathit{F}, rate of deformation tensor \dot{\mathrm{\epsilon }}, and spin tensor \dot{\omega } are given by:

(2)\mathit{F}=\mathit{1}+{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}},

(3)\dot{\epsilon }\text{=}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\nabla \dot{\mathit{u}}+{\left(\nabla \dot{\mathit{u}}\right)}^{}\right],

(4)\dot{\omega }\text{=}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\nabla \dot{\mathit{u}}-{\left(\nabla \dot{\mathit{u}}\right)}^{}\right],

in which \mathit{1} is the second order identity tensor, \nabla is the spatial gradient operator, \mathit{u} and \dot{\mathit{u}} represent the displacement and velocity, respectively. During the deformation the soil porosity n satisfies the following relationship resulting from the soil phase mass balance [⁵⁹]:

(5)n=\mathrm{1}-{\displaystyle \frac{\mathrm{1}-n_{\mathrm{0}}}{J}},J=\mathrm{\text{det}}(\mathit{F}),

with n_{\mathrm{0}} being the initial soil porosity.

The equation of motion and the continuity equation for the soil-water mixture considered herein are [⁶⁰]

(6)\nabla \cdot \mathit{\text{σ}}+\mathit{b}=\rho \ddot{\mathit{u}}\mathrm{\text{in}}\mathrm{\Omega },

(7){({{\rho }_{\mathrm{w}}}_{}n{S_{\mathrm{r}}}_{})}_{,t}+\nabla \cdot ({{\rho }_{\mathrm{w}}}_{}{{\mathit{v}}_{\mathrm{w}}}_{})=\mathrm{0}\mathrm{\text{in}}\mathrm{\Omega },

where \ddot{\mathit{u}} is the acceleration, \mathit{\text{σ}} is the total Cauchy stress, \mathit{b} is the body force, {{\mathit{v}}_{}}_{\mathrm{w}} is the water velocity, the subscript with comma means a partial differentiation with respect to the argument behind the comma. It is noted that the total stress \mathit{\text{σ}} is related to the effective stress \widehat{\mathit{\text{σ}}} as [⁶¹]

\text{σ}\text{=}\text{σ}+\tilde{p}\mathrm{1},

(8)\mathit{\sigma }=\mathit{\sigma }+\tilde{p}\mathit{1},\tilde{p}={S_{\mathrm{r}}}_{}{p_{\mathrm{w}}}_{}+(\mathrm{1}-{S_{\mathrm{r}}}_{}){p_{\mathrm{a}}}_{},

with \tilde{p}, {p_{\mathrm{w}}}_{} and {p_{\mathrm{a}}}_{} being the mixture, water and air pore pressures, respectively. Here, the frequently used passive air pore pressure condition with {p_{\mathrm{a}}}_{}=0 is adopted in the subsequent development, which is a reasonable assumption for the near surface landslides.

Equation (6) is supplemented by the boundary and initial conditions for the soil and rainfall infiltration:

(9)\mathit{u}=\overline{\mathit{u}}\mathrm{\text{on}}{\mathrm{\Gamma }}^{\mathrm{g}},\text{σ}\cdot \mathit{n}=\overline{\mathit{t}}\mathrm{\text{on}}{\mathrm{\Gamma }}^{\mathrm{t}},

(10)h=\overline{h}\mathrm{\text{on}}{{\mathrm{\Gamma }}^{\mathrm{h}}}^{},{v_{\mathrm{w}}}_{}\cdot \mathit{n}=\overline{q}\mathrm{\text{on}}{{\mathrm{\Gamma }}^{\mathrm{w}}}^{},

(11)\begin{array}{ll}\mathit{u}(\mathit{x},\mathrm{0})={\mathit{u}}_{\mathrm{0}}(\mathit{x}),\dot{\mathit{u}}(\mathit{x},\mathrm{0})={\dot{\mathit{u}}}_{\mathrm{0}}(\mathit{x}),h(\mathit{x},\mathrm{0})=h_{\mathrm{0}}(\mathit{x})& \mathrm{\text{in}}\mathrm{\Omega }\end{array},

where h denotes the water head, {\mathrm{\Gamma }}^{\mathrm{g}} and {\mathrm{\Gamma }}^{\mathrm{t}} represent the essential and natural boundaries, respectively, \mathit{n} denotes the outward boundary normal, {}^{}{\mathrm{\Gamma }}^{\mathrm{h}} is the prescribed water head boundary, and {}^{}{\mathrm{\Gamma }}^{\mathrm{w}} is the flux boundary corresponding to the rainfall infiltration, respectively.

The variational formulation of Eqs. (6) and (7) reads as follows:

(12){\displaystyle {\int }_{\Omega }\delta \mathit{u}\cdot \rho \ddot{\mathit{u}}\text{d}\mathrm{\Omega }}={\displaystyle {\int }_{{\Gamma }^t}\delta \mathit{u}\cdot \overline{\mathit{t}}\text{d}\Gamma }+{\displaystyle {\int }_{\Omega }\delta \mathit{u}\cdot \mathit{b}\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\nabla \delta \mathit{u}:\mathit{\text{σ}}\mathrm{\text{d}}\mathrm{\Omega }},

(13)$$\begin{gathered} {\displaystyle {\int }_{\mathrm{\Omega }}\delta p_{\mathrm{w}}n{S}_{r,s}{\dot{p}}_{\mathrm{w}}\mathrm{d}\mathrm{\Omega }={\int }_{\mathrm{\Omega }}\delta p_{\mathrm{w}}\mathrm{t}\mathrm{r}(\dot{ϵ})S_{\mathrm{r}}d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}^{\mathrm{w}}}\delta p_{\mathrm{w}}\overline{q}\mathrm{d}\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}(\nabla \delta p_{\mathrm{w}})\cdot v_{\mathrm{w}}\mathrm{d}\mathrm{\Omega }, \\ } \end{gathered}$$

where s is the matrix suction defined as s={p_{\mathrm{a}}}_{}-{p_{\mathrm{w}}}_{} [⁵⁸].

Damage and plasticity relationships

In this study the damage behavior of soil skeleton is characterized by the following exponential damage law [⁵⁴–⁵⁶]:

(14)\xi =\mathcal{D}(e_{\mathrm{t}})=\mathrm{1}-(\mathrm{1}-c_{\mathrm{1}}){\displaystyle \frac{e_{\mathrm{0}}}{e_{\mathrm{t}}}}-c_{\mathrm{1}}\mathrm{\text{exp}}[c_{\mathrm{2}}(e_{\mathrm{0}}-e_{\mathrm{t}})],

where \xi is the damage index, c_{\mathrm{1}} and c_{\mathrm{2}} are the damage material parameters, e_{\mathrm{0}} is the damage threshold, e_t is an equivalent and cumulative strain-like measure of deformation defied as follows:

(15)e_{\mathrm{t}}=\sqrt{\mathrm{2}\mathrm{\Lambda }(\widehat{ϵ}),}\wedge (\widehat{ϵ})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\widehat{ϵ}:\mathit{C}:\widehat{ϵ},\widehat{ϵ}={\int }_{\mathrm{0}}^t\dot{ϵ}(t)\mathrm{d}t,

in which \widehat{ϵ} and

\Lambda

are the equivalent deformation and energy measures. \mathit{C} is the elasticity tensor given by

(16)\mathit{C}=\lambda \mathit{1}\otimes \mathit{1}+\mathrm{2}\mu \mathit{I},

where \lambda and \mu are the Lame’s constants and \mathit{I} is the fourth order symmetric identity tensor.

With the aid of the damage index

\xi

in Eq. (14), the true effective stress \overline{\sigma } can be defined as

(17)\overline{\mathit{\sigma }}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\xi }}\widehat{\mathit{\sigma }},

where \widehat{\mathit{\sigma }} is the effective stress tensor given by Eq. (8).

The plasticity response of the soil skeleton is described by the classical Drucker-Prager model with the true effective stress \overline{\sigma }:

(18)F(\overline{\mathit{\sigma }})=\Vert \mathrm{d}\mathrm{e}\mathrm{v}(\overline{\mathit{\sigma }})\Vert +A\mathrm{t}\mathrm{r}(\overline{\mathit{\sigma }})-Y\le \mathrm{0},

where \overline{\sigma } is the true effective stress, \mathrm{\text{dev}}(•) and \mathrm{\text{tr}}(•) are the deviatoric and trace operators, A and Y are the material constants:

(19)A={\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{6}}\sin {\varphi }^{\prime }}{\mathrm{3}(\mathrm{3}+\sin {\varphi }^{\prime })}},Y=\mathrm{3}Ac\tan {\varphi }^{\prime },

where c and {\varphi }^{\prime } are the soil cohesion and the effective friction angle, respectively. In computation, the relationship of c=c^{\prime }+s\mathrm{\text{tan}}({\varphi }^{\prime }/\mathrm{2}) is employed [⁵⁸], where s and c^{\prime } stand for the suction and the effective cohesion.

The rainfall seepage through the soil skeleton follows the Darcy’s law:

(20){{\mathit{v}}_{\mathrm{w}}}_{}=-{{\mathit{D}}_{\mathrm{w}}}_{}\nabla \left({\displaystyle \frac{{p_{\mathrm{w}}}_{}}{{{\gamma }_{\mathrm{w}}}_{}}}+z\right),

in which z is the elevation water head, {{\gamma }_{}}_{\mathrm{w}} is the water gravity, and {{\mathit{D}}_{}}_{\mathrm{w}} is the seepage conductivity tensor [⁶²]:

(21){{\mathit{D}}_{\mathrm{w}}}_{}={\displaystyle \frac{{a_{\mathrm{w}}}_{}{k_{\mathrm{w}\mathrm{s}}}_{}}{{a_{\mathrm{w}}}_{}+{({c_{\mathrm{w}}}_{}s)}^{{b_{}}_{\mathrm{w}}}}}\mathit{1},

where {a_{}}_{\mathrm{w}}, {b_{}}_{\mathrm{w}}, {c_{\mathrm{w}}}_{} are material constants and {k_{}}_{\mathrm{w}\mathrm{s}} is the conductivity coefficient corresponding to the fully saturated state. Due to the rainfall infiltration, the degree of saturation {S_{}}_{\mathrm{r}} also depends the suction s [⁶³]:

(22){S_{\mathrm{r}}}_{}={S_{\mathrm{n}}}_{}+{\displaystyle \frac{{a_{\mathrm{s}}}_{}({S_{\mathrm{r}\mathrm{s}}}_{}-{S_{\mathrm{r}\mathrm{i}}}_{})}{{a_{\mathrm{s}}}_{}+{({c_{\mathrm{s}}}_{}s)}^{{b_{}}_{\mathrm{s}}}}},

with {a_{\mathrm{s}}}_{}, {b_{\mathrm{s}}}_{}, {c_{\mathrm{s}}}_{} being the soil constants. {S_{}}_{\mathrm{r}\mathrm{i}} and {S_{}}_{\mathrm{r}\mathrm{s}} represent the residual and maximum degrees of saturation, respectively.

Meshfree formulation with two-level gradient smoothing

Meshfree shape function

In Lagrangian meshfree approximation [⁸], the initial configuration {\mathrm{\Omega }}_{\mathrm{0}} and its boundary {\mathrm{\Gamma }}_{\mathrm{0}} are discretized by a set of meshfree nodes

{\{{\mathit{X}}_I\}}_{I=\mathrm{1}}^{NP}

, NP is the total number of meshfree nodes. The approximant of the displacement vector \mathit{u}, written by {\mathit{u}}^h, takes the following form:

(23){\mathit{u}}^h(\mathit{X},t)={\displaystyle \sum _{X\in {\mathcal{S}}_I}^{}{\Psi }_I(\mathit{X}){\mathit{d}}_I(t)},

where {\Psi }_I(\mathit{X}) and {\mathit{d}}_I are the meshfree shape function and nodal coefficient vector associated with the meshfree node {\mathit{X}}_I, {\mathcal{S}}_I is the influence domain of {\Psi }_I(\mathit{X}).

In this study, the MLS/RK meshfree approximation is considered and accordingly the shape function {\Psi }_I(\mathit{X}) reads [⁶,⁸]:

(24){\Psi }_I(\mathit{X})={\mathit{p}}^T({\mathit{X}}_I-\mathit{X})\mathit{c}(\mathit{X}){\varphi }_a({\mathit{X}}_I-\mathit{X}),

where {\varphi }_a({\mathit{X}}_I-\mathit{X}) is the non-negative kernel function measured by a support size “a”, which vanishes when \mathit{X} lies out of the compact support or influence domain of {\varphi }_a({\mathit{X}}_I-\mathit{X}) or {\Psi }_I(\mathit{X}), denoted by {\mathcal{S}}_I as just mentioned previously. The C² continuous cubic B-spline function is chosen as the kernel function and its tensor product formulation is used to construct the 3D kernel function. \mathit{p}(\mathit{X}) is the nth order basis vector:

(25)\mathit{p}(\mathit{X})={{\left\{\begin{array}{lllllll}\mathrm{1}& X_{\mathrm{1}}& X_{\mathrm{2}}& X_{\mathrm{3}}& X_{\mathrm{1}}^{\mathrm{2}}& \mathrm{...}& X_{\mathrm{3}}^n\end{array}\right\}}^{\mathrm{T}}}^{}.

The unknown coefficient vector \mathit{c}(\mathit{X}) in Eq. (24) is determined by enforcing the so-called consistency condition:

(26){\displaystyle \sum _{X\in {\mathcal{S}}_I}^{}{\Psi }_I(\mathit{X})\mathit{p}({\mathit{X}}_I-\mathit{X})}=\mathit{p}(\mathrm{0}).

Substituting Eq. (24) into Eq. (26) gives

(27)\mathit{ℳ}(\mathit{X})\mathit{c}(\mathit{X})=\mathit{p}(\mathrm{0}),

in which \mathit{ℳ}(\mathit{X}) is the moment matrix:

(28)\mathit{ℳ}(\mathit{X})={\displaystyle \sum _{X\in {\mathcal{S}}_I}^{}\mathit{p}({\mathit{X}}_I-\mathit{X}){\mathit{p}}^T({\mathit{X}}_I-\mathit{X}){\varphi }_a({\mathit{X}}_I-\mathit{X})},

Consequently we arrive at \mathit{c}(\mathit{X})={\mathit{ℳ}}^{-\mathrm{1}}(\mathit{X})\mathit{p}(\mathrm{0}) and the meshfree shape function {\Psi }_I(\mathit{X}) in Eq. (24) can be rephrased as

(29){\Psi }_I(\mathit{X})={\mathit{p}}^{}(\mathrm{0}){\mathit{ℳ}}^{-\mathrm{1}}(\mathit{X})\mathit{p}({\mathit{X}}_I-\mathit{X}){\varphi }_a({\mathit{X}}_I-\mathit{X}),

Two-level Lagrangian gradient smoothing formulation

The two-level strain smoothing meshfree formulation [⁵¹] has been introduced to resolve the discretization sensitivity issue in damage and failure modeling. This approach is adopted here through a two-level Lagrangian gradient smoothing to model the 3D large deformation rainfall induced landslides. The rate formulation is employed and the two-level strain smoothing meshfree formulation are essentially reflected by the two-level gradient smoothing operation on the Lagrangian meshfree shape functions evaluated at nodes.

The one-level gradient smoothing of Lagrangian meshfree shape function {\Psi }_I at a node {\mathit{X}}_L, denoted by {\tilde{\Psi }}_{I,J}({\mathit{X}}_L), is defined as [⁴²]

(30){\tilde{\Psi }}_{I,J}({\mathit{X}}_L)={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{0}}}^{}{\Phi }^{[\mathrm{1}]}({\mathit{X}}_L;\mathit{X}){\displaystyle \frac{\partial {\Psi }_I(\mathit{X})}{\partial X_J}}\mathrm{\text{d}}\mathrm{\Omega }},

where {\Phi }^{[\mathrm{1}]}({\mathit{X}}_L;\mathit{X}) is the first level weight function and it can be conveniently selected as [⁴¹]

(31){\mathrm{\Phi }}^{[\mathrm{1}]}({\mathit{X}}_L;\mathit{X})=\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{V_L}}& \mathit{X}\in {\mathrm{\Omega }}_{\mathrm{0}L}\\ \mathrm{0}& \mathit{X}\notin {\mathrm{\Omega }}_{\mathrm{0}L}\end{array},

where {\mathrm{\Omega }}_{\mathrm{0}L} denotes the nodal representative domain related to the node {\mathit{X}}_L as shown in Fig. 1. V_L is the volume of {\mathrm{\Omega }}_{\mathrm{0}L}. Substituting Eq. (31) into Eq. (30) and invoking the divergence theorem yields

(32){\tilde{\Psi }}_{I,J}({\mathit{X}}_L)\text{=}{\displaystyle \frac{\mathrm{1}}{V_L}}{\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{0}L}}^{}{\displaystyle \frac{\partial {\Psi }_I(\mathit{X})}{\partial X_J}}\mathrm{\text{d}}\mathrm{\Omega }}={\displaystyle \frac{\mathrm{1}}{V_L}}{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{0}L}}^{}{\Psi }_I(\mathit{X})N_J(\mathit{X})\mathrm{\text{d}}\mathrm{\Gamma }},

where {\mathrm{\Gamma }}_{\mathrm{0}L} denotes the surface of {\mathrm{\Omega }}_{\mathrm{0}L} and \mathit{N} stands for the outward normal of {\mathrm{\Gamma }}_{\mathrm{0}L}. It is noted that the one-level smoothed gradient of meshfree shape function, namely, {\tilde{\Psi }}_{I,J}({\mathit{X}}_L) in Eq. (32), constitutes the basis for the widely used method of stabilized conforming nodal integration [⁴¹,⁴²,⁵¹].

In numerical implementation, as shown in Fig. 1, the boundary integral in Eq. (32) can be conveniently computed as

(33){\tilde{\Psi }}_{I,J}({\mathit{X}}_L)={\displaystyle \frac{\mathrm{1}}{V_L}}{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{0}L}}^{}{\Psi }_I(\mathit{X})N_J\mathrm{\text{d}}\mathrm{\Gamma }}={\displaystyle \frac{\mathrm{1}}{V_L}}{\displaystyle \sum _{C=\mathrm{1}}^{NS}{A}_C{\Psi }_I({\mathit{X}}_C)N_J},

where {\mathit{X}}_C and A_C are the center and area of a generic boundary surface {\mathrm{\Gamma }}_{\mathrm{0}C}. NS represents the number of sub-surfaces of {\mathrm{\Gamma }}_{\mathrm{0}C} used for numerical integration. A_C and N_I\text{'}s can be conveniently expressed by the vertex coordinates of the surface {\mathrm{\Gamma }}_{\mathrm{0}C} as follows [⁶⁴]:

(34)\{\begin{array}{l}{A}_C=\sqrt{D_{\mathrm{12}}^{\mathrm{2}}+D_{\mathrm{23}}^{\mathrm{2}}+D_{\mathrm{31}}^{\mathrm{2}}}\\ N_{\mathrm{1}}={\displaystyle \frac{D_{\mathrm{23}}}{A_C}},N_{\mathrm{2}}={\displaystyle \frac{D_{\mathrm{31}}}{A_C}},N_{\mathrm{3}}={\displaystyle \frac{D_{\mathrm{12}}}{A_C}}\\ D_{IJ}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{M=\mathrm{0}}^{NV}[X_{IM}^C{X}_{J(M+\mathrm{1})}^C-X_{I(M+\mathrm{1})}^C{X}_{JM}^C]}\end{array},

where {\mathit{X}}_M^C stands for the vertex point for {\mathrm{\Gamma }}_C, and NV denotes the total number of vertices of {\mathrm{\Gamma }}_C. X_{IM}^C represents the Ith component of the position vector {\mathit{X}}_M^C.

In order to properly resolve the discretization sensitivity issue, a two-level smoothed gradient is introduced based upon a further gradient smoothing operation on the one-level smoothed gradient defined by Eq. (33) [⁵¹]:

(35)\begin{array}{l}{\overbrace{\mathit{\psi }}}_{I,J}({\mathit{X}}_K)={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{0}}}^{}{\mathrm{\Phi }}^{[\mathrm{2}]}({\mathit{X}}_K;\mathit{X}){\tilde{\mathrm{\Psi }}}_{I,J}(\mathit{X})\mathrm{\text{d}}\mathrm{\Omega }}\\ ={\displaystyle \sum _{X_K\in {\mathcal{S}}_L}^{}{V}_L{\mathrm{\Phi }}^{[\mathrm{2}]}({\mathit{X}}_K;{\mathit{X}}_L){\tilde{\Psi }}_{I,J}({\mathit{X}}_L)},\end{array}

where for convenience the second level weight function {\Phi }^{[\mathrm{2}]}({\mathit{X}}_K;\mathit{X}) is selected as a weighted form of the meshfree shape function:

(36){\mathrm{\Phi }}^{[\mathrm{2}]}({\mathit{X}}_K;\mathit{X})={\displaystyle \frac{\mathrm{1}}{V_L}}{\mathrm{\Psi }}_L({\mathit{X}}_K).

Thus, the two-level smoothed gradient of meshfree shape function finally becomes

(37){\overbrace{\mathrm{\Psi }}}_{I,J}({\mathit{X}}_K)={\displaystyle \sum _{X_K\in {\mathcal{S}}_L}^{}{\mathrm{\Psi }}_L({\mathit{X}}_K){\tilde{\mathrm{\Psi }}}_{I,J}({\mathit{X}}_L)}.

It is noted that both the one-level and two-level smoothed nodal gradients of meshfree shape function, i.e., {\overbrace{\mathrm{\Psi }}}_{I,J}({\mathit{X}}_K) and {\overbrace{\mathrm{\Psi }}}_{I,J}({\mathit{X}}_L), do not involve the time-consuming computation of the direct derivatives of meshfree shape function, which is preferable from the computational efficiency point of view. Meanwhile, due to the successive smoothing operation by the weight functions, the support size or influence domain of the two-level smoothed gradient of meshfree shape function {\overbrace{\mathrm{\Psi }}}_{I,J}({\mathit{X}}_K) is usually larger than those of the direct derivative and the one-level smoothed gradient, which is clearly illustrated in Fig. 2. Figure 3 presents an efficiency comparison for the computation of various gradients of meshfree shape function with respect to different discretizations of a 3D cube, where the actual normalized support size is 3.0. The results reveal that for a given final influence domain of the shape function gradient, the computation of the two-level smoothed gradient could be more efficient than the standard gradient of meshfree shape function, and the efficiency of the one- and two-level smoothed gradient calculations is comparable. This is because of the no direct derivative evaluation and the relative smaller support size in each gradient smoothing level.

Discrete meshfree equations with two-level smoothed gradient

Introducing the Lagrangian meshfree approximation into the primary dependent variables of the soil-water mixture, say, the displacement \mathit{u} and the water pore pressure {p_{}}_{\mathrm{w}}, leads to the following expressions:

(38){\dot{\mathit{u}}}^h(\mathit{X},t)={\displaystyle \sum _{X\in {\mathcal{S}}_I}^{}{\mathrm{\Psi }}_I(\mathit{X}){\dot{\mathit{d}}}_I(t)},{\ddot{\mathit{u}}}^h(\mathit{X},t)={\displaystyle \sum _{\mathit{X}\in {\mathcal{S}}_I}^{}{\mathrm{\Psi }}_I(\mathit{X}){\ddot{\mathit{d}}}_I(t)},

(39)p_w^h(\mathit{X},t)={\displaystyle \sum _{X\in {\mathcal{S}}_I}^{}{\mathrm{\Psi }}_I(\mathit{X})q_I(t)},{\dot{p}}_{\mathrm{w}}^h(\mathit{X},t)={\displaystyle \sum _{\mathit{X}\in {\mathcal{S}}_I}^{}{\mathrm{\Psi }}_I(\mathit{X}){\dot{q}}_I(t)}.

Here, equal order basis functions are employed for both displacement and pore pressure field approximations, the stability of this type of meshfree formulation within the stabilized conforming nodal integration framework has been studied through the consolidation analysis of pressure field [⁶⁵,⁶⁶]. At the same time, the two-level smoothed gradient of meshfree shape function is employed to construct the two-level smoothed deformation gradient:

(40){\overbrace{\mathit{F}}}^h({\mathit{X}}_K,t)=\mathit{1}+{\displaystyle \sum _{{\mathit{X}}_K\in {\mathcal{S}}_{IL}}^{}{\mathit{d}}_I(t)\otimes {\overbrace{\mathrm{\Psi }}}_{I,X}({\mathit{X}}_K)},

where {\mathcal{S}}_{IL}={\mathcal{S}}_I\cup {\mathcal{S}}_L with {\mathcal{S}}_I\cap {\mathcal{S}}_L\ne \varnothing, and it is noted the relationship of Eq. (37) is needed to compute Eq. (40).

According to Eqs. (38)–(40), the rate of deformation tensor and the two-level water pressure gradient evaluated at node {\mathit{X}}_K take the following forms:

(41){\dot{\epsilon }}^h({\mathit{X}}_K,t)={\displaystyle \sum _{X_K\in {\mathcal{S}}_{IL}}^{}{\overbrace{\mathit{B}}}_I^{}({\mathit{X}}_K){\dot{\mathit{d}}}_I(t)},

(42)\nabla p^h({\mathit{X}}_K,t)={\displaystyle \sum _{{\mathit{X}}_K\in {\mathcal{S}}_{IL}}^{}{\overbrace{\mathit{B}}}_I^{}({\mathit{X}}_K)q_I(t)},

where \dot{ϵ}={\{{\dot{ϵ}}_{\mathrm{11}}{\dot{ϵ}}_{\mathrm{22}}{\dot{ϵ}}_{\mathrm{33}}{\mathrm{2}\dot{ϵ}}_{\mathrm{12}}{\mathrm{2}\dot{ϵ}}_{\mathrm{13}}{\mathrm{2}\dot{ϵ}}_{\mathrm{23}}\}}^T, \nabla {p_{\mathrm{w}}}_{}^{}={\{\begin{array}{lll}{p}_{\mathrm{w},x_{\mathrm{1}}}^{}& p_{\mathrm{w},x_{\mathrm{2}}}^{}& p_{\mathrm{w},x_{\mathrm{3}}}^{}{\}}^{}\end{array}}^{\mathrm{T}}. The gradient matrices {\overbrace{\mathit{B}}}_I({\mathit{X}}_K) and {\overbrace{\mathit{B}}}_I^w({\mathit{X}}_K) are given by

(43){\overbrace{\mathit{B}}}_I({\mathit{X}}_K)=\left[\begin{array}{lll}{\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{1}}}({\mathit{X}}_K)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{2}}}({\mathit{X}}_K)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{3}}}({\mathit{X}}_K)\\ {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{2}}}({\mathit{X}}_K)& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{1}}}({\mathit{X}}_K)& \mathrm{0}\\ {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{3}}}({\mathit{X}}_K)& \mathrm{0}& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{1}}}({\mathit{X}}_K)\\ \mathrm{0}& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{3}}}({\mathit{X}}_K)& {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{2}}}({\mathit{X}}_K)\end{array}\right],

(44){\overbrace{\mathit{B}}}_I^{}({\mathit{X}}_K)={\overbrace{\mathrm{\Psi }}}_{I,\mathit{x}}({\mathit{X}}_K)\text{=}\left\{\begin{array}{l}{\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{1}}}({\mathit{X}}_K)\\ {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{2}}}({\mathit{X}}_K)\\ {\overbrace{\mathrm{\Psi }}}_{I,x_{\mathrm{3}}}({\mathit{X}}_K)\end{array}\right\},

where {\overbrace{\mathrm{\Psi }}}_{I,\mathit{x}}({\mathit{X}}_K) stands for the two-level smoothed gradient referring to the spatial coordinate \mathit{x}, and it is related to the two-level smoothed gradient referring to the material coordinate \mathit{X}, say, {\overbrace{\mathrm{\Psi }}}_{I,\mathit{X}}({\mathit{X}}_K), through the following relationship:

(45){\overbrace{\mathrm{\Psi }}}_{I,\mathit{x}}({\mathit{X}}_K)={\overbrace{\mathrm{\Psi }}}_{I,\mathit{X}}({\mathit{X}}_K)\cdot {\overbrace{\mathit{F}}}^h{}^{-\mathrm{1}}({\mathit{X}}_K,t)

Substituting Eqs. (38)–(42) into the weak form of Eq. (12) with a nodal integration leads to the following semi-discrete equation:

(46)\left[\begin{array}{ll}\mathit{M}& \mathrm{0}\\ \mathrm{0}& \mathit{M}\end{array}\right]\left[\begin{array}{l}\ddot{\mathit{d}}\\ \dot{\mathit{q}}\end{array}\right]+\left[\begin{array}{ll}\mathrm{0}& \mathrm{0}\\ -\mathit{P}& -\mathit{S}\end{array}\right]\left[\begin{array}{l}\dot{\mathit{d}}\\ \mathit{q}\end{array}\right]=\left[\begin{array}{l}{{\mathit{f}}^{\mathrm{e}\mathrm{x}\mathrm{t}}}^{}-{{\mathit{f}}^{\mathrm{i}\mathrm{n}\mathrm{t}}}^{}\\ {{\mathit{f}}^{\mathrm{w}-\mathrm{e}\mathrm{x}\mathrm{t}}}_{}^{}\end{array}\right],

in which the entries of various matrices and vectors are given by

(47)\{\begin{array}{l}{\mathit{M}}_{IJ}={\displaystyle \sum _{K=\mathrm{1}}^{NP}{\mathrm{\Psi }}_I({\mathit{X}}_K){\Psi }_J({\mathit{X}}_K)\rho ({\mathit{X}}_K)\mathit{1}\widehat{J}({\mathit{X}}_K)V_K}\\ {\mathit{f}}_I^{\mathrm{e}\mathrm{x}\mathrm{t}}={\displaystyle \sum _{S=\mathrm{1}}^{NB}{\mathrm{\Psi }}_I({\mathit{X}}_S){\overline{\mathit{t}}}_{\mathrm{0}}({\mathit{X}}_S)A_S}+{\displaystyle \sum _{K=\mathrm{1}}^{NP}{\mathrm{\Psi }}_I({\mathit{X}}_K)\mathit{b}({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K}\\ {\mathit{f}}_I^{\mathrm{int}}={\displaystyle \sum _{K=\mathrm{1}}^{NP}{\overbrace{\mathit{B}}}_I^T({\mathit{X}}_K)\mathit{\text{σ}}({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K}\end{array},

(48)\{\begin{array}{l}{M}_{IJ}^{\mathrm{w}}={\displaystyle \sum _{K=\mathrm{1}}^{NP}{\mathrm{\Psi }}_I^{}({\mathit{X}}_K)n({\mathit{X}}_K)S_{r,s}({\mathit{X}}_K){\mathrm{\Psi }}_J({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K^{}}\\ f_I^{w-\mathrm{w}\mathrm{x}\mathrm{t}}={\displaystyle \sum _{S=\mathrm{1}}^{NB}{\mathrm{\Psi }}_I({\mathit{X}}_S){\overline{q}}_{\mathrm{0}}({\mathit{X}}_S)A_S}+{\displaystyle \sum _{K=\mathrm{1}}^{NP}{\overbrace{\mathit{B}}}_I^{}({\mathit{X}}_K){\mathit{D}}_{}({\mathit{X}}_K)\nabla z({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K^{}}\\ P_{IJ}\text{=}{\displaystyle \sum _{K=\mathrm{1}}^{NP}{\mathrm{\Psi }}_I^{}({\mathit{X}}_K)S_{\mathrm{r}}({\mathit{X}}_K){\mathit{m}}_{}^{}{\overbrace{\mathit{B}}}_J^{}({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K^{}}\\ S_{IJ}^{}={\displaystyle \sum _{K=\mathrm{1}}^{NP}{\widehat{\mathit{B}}}_I^{\mathrm{w}\mathrm{T}}(X_K){\gamma }_{}^{-\mathrm{1}}(X_K)D_{}(X_K){\widehat{\mathit{B}}}_J^{\mathrm{w}}({\mathit{X}}_K)\widehat{J}({\mathit{X}}_K)V_K^{}}\end{array},

where \mathit{\text{σ}}={{\left\{\begin{array}{llllll}{\sigma }_{\mathrm{11}}^{}& {\sigma }_{\mathrm{22}}^{}& {\sigma }_{\mathrm{33}}^{}& {\sigma }_{\mathrm{12}}^{}& {\sigma }_{\mathrm{13}}^{}& {\sigma }_{\mathrm{23}}^{}\end{array}\right\}}^{\mathrm{T}}}^{}, \mathit{m}={{\left\{\begin{array}{llllll}\mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right\}}^{\mathrm{T}}}^{}. {\overline{\mathit{t}}}_{\mathrm{0}} and {\bar{q}}_{\mathrm{0}} are the traction and seepage flux referring to the initial configuration. NB represents the number of boundary integration points and A_S is the corresponding surface integration weight, and \widehat{J}({\mathit{X}}_K)=\mathrm{\text{det}}[{\widehat{\mathit{F}}}^h({\mathit{X}}_K)].

Further introducing the Newmark method and generalized trapezoidal rule for the temporal discretizations of the displacement field and the pore pressure field, respectively, into Eq. (46) yields the fully discrete equation:

(49)\left[\begin{array}{ll}\mathit{M}& \mathrm{0}\\ -\Delta t\gamma \mathit{P}& {\mathit{M}}^{}-\Delta t\lambda \mathit{S}\end{array}\right]\left\{\begin{array}{l}{\ddot{\mathit{d}}}_{n+\mathrm{1}}\\ {\dot{\mathit{q}}}_{n+\mathrm{1}}\end{array}\right\}=\left\{\begin{array}{l}{\mathit{f}}_{n+1}^{\mathrm{e}\mathrm{x}\mathrm{t}}-{\mathit{f}}_{n+1}^{\mathrm{int}}\\ {\mathit{f}}_{n+1}^{\mathrm{w}-\mathrm{e}\mathrm{x}\mathrm{t}}+\mathit{P}{\widetilde{\dot{\mathit{d}}}}_{n+\mathrm{1}}+\mathit{S}{\tilde{\mathit{q}}}_{n+\mathrm{1}}\end{array}\right\},

with

(50)\{\begin{array}{l}{\tilde{\mathit{d}}}_{n+\mathrm{1}}={\mathit{d}}_n^{}+\Delta t{\dot{\mathit{d}}}_n^{}+{\displaystyle \frac{\Delta t^{\mathrm{2}}}{\mathrm{2}}}(\mathrm{1}-\mathrm{2}\beta ){\ddot{\mathit{d}}}_n^{}\\ {\widetilde{\dot{\mathit{d}}}}_{n+\mathrm{1}}={\dot{\mathit{d}}}_n^{}+\Delta t(\mathrm{1}-\gamma ){\ddot{\mathit{d}}}_n^{}\\ {\tilde{\mathit{q}}}_{n+\mathrm{1}}={\mathit{q}}_n^{}+\Delta t(\mathrm{1}-\alpha ){\dot{\mathit{q}}}_n^{}\\ {\mathit{d}}_{n+\mathrm{1}}^{}={\tilde{\mathit{d}}}_{n+\mathrm{1}}+\beta \Delta t^{\mathrm{2}}{\ddot{\mathit{d}}}_{n+\mathrm{1}}^{}\\ {\dot{\mathit{d}}}_{n+\mathrm{1}}^{}={\widetilde{\dot{\mathit{d}}}}_{n+\mathrm{1}}^{}+\Delta t\gamma {\ddot{\mathit{d}}}_{n+\mathrm{1}}^{}\\ {\mathit{q}}_{n+\mathrm{1}}={\tilde{\mathit{q}}}_{n+\mathrm{1}}+\Delta t\alpha {\dot{\mathit{q}}}_{n+\mathrm{1}}\end{array},

where \beta and \gamma are the Newmark parameters and \alpha is the parameter of generalized trapezoidal rule. In the subsequent numerical examples, the central difference and forward-Euler explicit schemes with \beta =\mathrm{0}, \gamma =\mathrm{1}/\mathrm{2} and \alpha =\mathrm{1}/\mathrm{2} are employed.

Deformation and stress update

In a strain-driven computation, we would like to advance the field variables \{{\mathit{x}}_n,{\widehat{\mathit{\sigma }}}_n,{\xi }_n,e_n\} at t=t_n to their counterparts \{{\mathit{x}}_{n+\mathrm{1}},{\widehat{\mathit{\sigma }}}_{n+\mathrm{1}},{\xi }_{n+\mathrm{1}},e_{n+\mathrm{1}}\} at t=t_{n+\mathrm{1}} provided the displacement increment \Delta {\mathit{u}}^h. To achieve a second order accuracy, the intermediate configuration {\mathit{x}}_{n+\mathrm{1}/\mathrm{2}}={\mathit{x}}_n+\Delta {\mathit{u}}^h/\mathrm{2} is employed to compute the relative displacement gradient \mathit{G}:

(51)\mathit{G}={\displaystyle \frac{\partial \Delta {\mathit{u}}^h}{\partial {\mathit{x}}_{n+\mathrm{1}/\mathrm{2}}}}=\left[{\displaystyle \sum _{I=\mathrm{1}}^{NP}\Delta {\mathit{d}}_I(t)\otimes {\widehat{\mathrm{\Psi }}}_{I,\mathit{X}}(\mathit{X})}\right]\cdot {\left({\widehat{\mathit{F}}}_{n+\mathrm{1}/\mathrm{2}}^h\right)}^{-\mathrm{1}},

Consequently, the incremental strain and spin tensors \Delta \text{ε} and \Delta \text{ω} corresponding to Eqs. (3) and (4) within the time interval \Delta t=t_{n+\mathrm{1}}-t_n become:

(52)\Delta \text{ε}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{G}+{{\mathit{G}}^{\mathrm{T}}}^{}),

(53)\Delta \text{ω}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{G}+{{\mathit{G}}^{\mathrm{T}}}^{}).

For the damage evolution, according to Eq. (14) and (15) the damage index at the time step t_{n+\mathrm{1}} is evaluated as [⁵⁵]

(54){\xi }_{n+\mathrm{1}}\text{=}\{\begin{array}{l}{\xi }_n{e}_{n+\mathrm{1}}-r_n\le \mathrm{0}\\ \mathcal{D}(e_{n+\mathrm{1}})\mathrm{\text{otherwise}}\end{array},

with

(55)r_{n+\mathrm{1}}=\mathrm{\text{max}}\{r_n,e_{n+\mathrm{1}}\},r_{\mathrm{0}}=e_{\mathrm{0}},

(56)e_{n+\mathrm{1}}=\sqrt{\mathrm{2}\wedge \left({\widehat{ϵ}}_{n+\mathrm{1}}\right),}{\widehat{ϵ}}_{n+\mathrm{1}}={\widehat{ϵ}}_n+\mathrm{\Delta }ϵ\mathrm{.}

As mentioned earlier, the true effective stress \overline{\sigma } is rationally employed for Drucker-Prager plasticity model. The stress update follows the strain-based damage formulation [⁵⁵] and objective integration using a rotated configuration [⁵⁷]. With the aid of Eq. (53), the rotation tensors are obtained as [⁵⁷]

(57)\{\begin{array}{l}{\mathit{R}}_{n+\mathrm{1}}=\mathrm{\text{exp}}(\Delta \text{ω}){\mathit{R}}_n\\ {\mathit{R}}_{n+\mathrm{1}/\mathrm{2}}=\mathrm{\text{exp}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Delta \text{ω}\right){\mathit{R}}_n\\ {\mathit{R}}_{\mathrm{0}}=\mathit{1}\end{array},

(58)\{\begin{array}{l}\mathit{Q}={\mathit{R}}_{n+\mathrm{1}}{\mathit{R}}_n^{\mathrm{T}}\\ \tilde{Q}={\mathit{R}}_{n+\mathrm{1}}{\mathit{R}}_{n+1/2}^{\mathrm{T}}\end{array},

The stress predictor {\overline{\sigma }}_{n+\mathrm{1}}^{trial} then reads:

(59)\{\begin{array}{l}{\widehat{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}=\mathit{Q}{\overline{\mathit{\sigma }}}_n{\mathit{Q}}^T+\mathit{C}:\left[\tilde{\mathit{Q}}\mathrm{\Delta }ϵ{\tilde{\mathit{Q}}}^T\right]\\ {\overline{\mathit{\sigma }}}_n={\displaystyle \frac{{\widehat{\mathit{\sigma }}}_n}{\mathrm{1}-{\xi }_n}}\end{array}

Subsequently, if F\left({\overline{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\right)\le \mathrm{0}, the trial stress {\overline{\sigma }}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}} is set to be {\overline{\mathit{\sigma }}}_{n+\mathrm{1}}, otherwise, the plasticity corrector is carried out as follows:

(60){\overline{\mathit{\sigma }}}_{n+\mathrm{1}}={\overline{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}-\mathrm{\Delta }{\overline{\gamma }}_{n+\mathrm{1}}(\mathrm{2}\mu {\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{v}\left({\overline{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\right)}{\Vert {\overline{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\Vert }}+\mathrm{3}KA\mathrm{\Delta }{\lambda }_{n+\mathrm{1}}\mathit{1}),

with

(61)\mathrm{\Delta }{\overline{\gamma }}_{n+\mathrm{1}}={\displaystyle \frac{F\left({\overline{\mathit{\sigma }}}_{n+\mathrm{1}}^{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}}\right)}{\mathrm{2}\mu +\mathrm{9}{KA}^{\mathrm{2}}}},

where K is the bulk modulus. Finally, the nominal effective stress {\widehat{\mathit{\sigma }}}_{n+\mathrm{1}}=(\mathrm{1}-{\xi }_{n+\mathrm{1}}){\widehat{\mathit{\sigma }}}_{n+\mathrm{1}} is then returned to the next load step.

Numerical examples

Simulation of Kashiwabara slope failure experiment

The first problem considered here is a rainfall experiment performed by Kawamura et al. [⁶⁷] on the Kashiwabara volcanic soil slope. The experimental setup is shown in Fig. 4(a), where the soil container is 2000 mm in length, 700 mm in depth and 600 mm in width. A rainfall with an intensity of 100 mm/h is applied at the top of the container. Figure 4(b) describes the computational model for meshfree simulation. The material properties of the soil are given as: Young’s modulus E=24.5 MPa, Poisson’s ratio \nu =\mathrm{0.34}, soil skeleton density ρ_s=2340 kg/m³, effective cohesion c^{\prime }=\mathrm{1.08335} kPa, effective friction angle {\varphi }^{\prime }={\mathrm{40}}^{\circ }, initial soil porosity n_{\mathrm{0}}=\mathrm{78.1}\text{\%}. The damage coefficients are e_{\mathrm{0}}=\mathrm{1.0}, c_{\mathrm{1}}=\mathrm{0.58}, and c_{\mathrm{2}}=\mathrm{0.04}. The parameters used in Eqs. (21) and (22) are: {a_{\mathrm{w}}}_{}=\mathrm{800.0}, {b_{\mathrm{w}}}_{}=\mathrm{0.8}, {c_{\mathrm{w}}}_{}=\mathrm{0.009}, {k_{\mathrm{w}\mathrm{s}}}_{}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}\text{m/s}, {a_{\mathrm{s}}}_{}=\mathrm{1.0}, {b_{\mathrm{s}}}_{}=\mathrm{0.9}, {c_{\mathrm{s}}}_{}=\mathrm{0.00567}, S_{ri}=\mathrm{0.11}, S_{rs}=\mathrm{1.0}. The initial degree of saturation S_{r\mathrm{0}} is described based upon the experimental data about 20 s prior to the failure [⁶⁷], i.e., 46.0% near the slope surface and 31.7% near the left bottom corner of the cross section.

The meshfree discretization for this problem is shown in Fig. 5, where the total number of meshfree nodes is 7425. A linear basis function with a normalized support size of 1.5 is used to construct the meshfree shape functions for both displacement and pore pressure approximations. The step for time integration is \Delta t=\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{4}} s. The progressive saturation and failure evolutions in the slope are presented in Figs. 5 and 6, where the failure initiation and propagation in the slope are clearly simulated by the proposed two-level gradient smoothing meshfree method. Moreover, the slip line provides an important and critical way to measure the slope failure or landslide. Thus, the numerical and experimental slip lines are compared in Fig. 7. The results show a reasonably good agreement between the slip line by the proposed meshfree simulation and that observed in the experiment.

Simulation of Yangbaodi landslide

The second benchmark for meshfree simulation is the Yangbaodi landslide that is located in Shenzhen, China [⁶⁸]. Figure 8(a) gives an overall view of the Yangbaodi landslide and the geological profile with layered granite and granitic residual soil is described in Fig. 8(b). The landslide has a plane length of 170 m, width of about 40 m. A rainfall with an intensity of 20 mm/h is applied at the slope surface. The soil properties used herein are: Young’s modulus E=\mathrm{24.7} MPa, Poisson’s ratio \nu =\mathrm{0.35}, soil skeleton density ρ_s=2670 kg/m³, initial suction s_{\mathrm{0}}=\mathrm{0.2} kPa, initial soil porosity n_{\mathrm{0}}=\mathrm{50}\text{\%}, initial degree of saturation S_{\mathrm{r}0}=\mathrm{45.85}\%, effective friction angle {\varphi }^{\prime }={\mathrm{25}}^{\circ }. The soil is quite loosely and its effective cohesion c^{\prime } is set to be zero. The saturation and seepage parameters are given as: {a_{\mathrm{w}}}_{}=\mathrm{1000.0}, {b_{\mathrm{w}}}_{}=\mathrm{1.7}, {c_{\mathrm{w}}}_{}=\mathrm{0.01}, {k_{\mathrm{w}\mathrm{s}}}_{}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}m/s, {a_{\mathrm{s}}}_{}=\mathrm{1.0}, {b_{\mathrm{s}}}_{}=\mathrm{1.0}, {c_{\mathrm{s}}}_{}=\mathrm{0.008}, S_{\mathrm{r}\mathrm{i}}=\mathrm{0.12}, {S_{\mathrm{r}\mathrm{s}}}_{}=\mathrm{1.0}. Moreover, the parameters of e_{\mathrm{0}}=\mathrm{100}, c_{\mathrm{1}}=\mathrm{0.01}, and c_{\mathrm{2}}=\mathrm{5} are employed for the damage evolution.

Since the granite foundation is nearly rigid compared to the upper soil blanket, thus only the soil blanket is modeled in the meshfree computation. As shown in Figs. 8(c) and 8(d), two meshfree discretizations with 117738 and 195027 nodes are employed for this problem, respectively. The meshfree shape functions employ a linear basis function and the normalized support size is 1.2. The time increment is \Delta t=\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{4}} s. The progressive saturation and failure behaviors of Yangbaodi landslide are illustrated by the meshfree results with 195027 nodes in Figs. 9 and 10, where large deformation damage and failure behavior is properly simulated by the proposed 3D meshfree formulation. Figure 11 presents a comparison between the meshfree results using different discretizations and the in situ survey regarding the final slope surface after landslide. It is observed that the two meshfree discretizations produce convergent results, and the final slope surface is essentially captured by the proposed meshfree method, as demonstrates the capability of the present 3D two-level gradient smoothing meshfree method to deal with the rainfall induced complex landslide failure.

Conclusions

A 3D two-level gradient smoothing meshfree formulation was developed for landslide simulations driven by rainfalls. The two-level gradient smoothing was carried out in the context of Lagrangian meshfree approximation and the corresponding 3D implementation details was given as well. It was shown that for a given final domain of influence, the computation of two-level smoothed gradient of meshfree shape function is less costly than the standard gradient calculation due to the fact of no direct derivative evaluation and the relatively smaller support size with respect to each gradient smoothing level. Subsequently, the two-level smoothed gradient of meshfree shape function was employed to discretize the coupled equations of soil motion and rainfall seepage using a nodal integration, and formulate a regularized meshfree approach for 3D large deformation simulations of rainfall induced landslides. The soil failure is characterized by the exponential damage law and the pressure dependent plasticity model with the true effective stress measure. Meanwhile, the large deformation soil behavior of landslides was accommodated by the rate formulation. A rotationally neutralized objective integration algorithm was used to perform the objective stress update. Two typical numerical examples were presented to demonstrate the effectiveness of the proposed method. It was seen that the results by the proposed two-level gradient smoothing meshfree formulation agree quite well with the experimental and in situ observations, and thus adequately validate the capability of the present 3D meshfree method for rainfall induced landslide simulations.

References

[1] Zienkiewicz O C, Taylor R L, Fox D D. The Finite Element Method for Solid and Structural Mechanics. 7th ed. Oxford: Butterworth-Heinemann, 2013

[2] Lucy L B. A numerical approach to the testing of the fission hypothesis. Astronomical Journal, 1977, 82: 1013–1024

[3] Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 375–389

[4] Liu M B, Liu G R. Smoothed particle hydrodynamics (SPH): An overview and recent developments. Archives of Computational Methods in Engineering, 2010, 17(1): 25–76

[5] Belytschko T, Lu Y Y, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[6] Liu W K, Jun S, Zhang Y F. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids, 1995, 20(8–9): 1081–1106

[7] Belytschko T, Lu Y Y, Gu L. Crack propagation by element-free Galerkin methods. Engineering Fracture Mechanics, 1995, 51(2): 295–315

[8] Chen J S, Pan C, Wu C T, Liu W K. Reproducing kernel particle methods for large deformation analysis of non-linear structures. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1–4): 195–227

[9] Libersky L D, Randles P W, Carney T C, Dickinson D L. Recent improvements in SPH modeling of hypervelocity impact. International Journal of Impact Engineering, 1997, 20(6–10): 525–532

[10] Liu W K, Jun S. Multiple-scale reproducing kernel particle methods for large deformation problems. International Journal for Numerical Methods in Engineering, 1998, 41(7): 1339–1362

[11] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[12] Vidal Y, Bonet J, Huerta A. Stabilized updated Lagrangian corrected SPH for explicit dynamic problems. International Journal for Numerical Methods in Engineering, 2007, 69(13): 2687–2710

[13] Wang D, Li Z, Li L, Wu Y. Three dimensional efficient meshfree simulation of large deformation failure evolution in soil medium. Science China. Technological Sciences, 2011, 54(3): 573–580

[14] Ren B, Li S, Qian J, Zeng X. Meshfree simulations of spall fracture. Computer Methods in Applied Mechanics and Engineering, 2011, 200(5–8): 797–811

[15] Wu Y, Wang D, Wu C T. Three dimensional fragmentation simulation of concrete structures with a nodally regularized meshfree method. Theoretical and Applied Fracture Mechanics, 2014, 72: 89–99

[16] Drathi R, Das A J M, Rangarajan A. Meshfree simulation of concrete structures and impact loading. International Journal of Impact Engineering, 2016, 91: 194–199

[17] Wu C T, Wu Y, Crawford J E, Magallanes J M. Three-dimensional concrete impact and penetration simulations using the smoothed particle Galerkin method. International Journal of Impact Engineering, 2017, 106: 1–17

[18] Atluri S N, Shen S P. The Meshless Local Petrov-Galerkin (MLPG) Method. Henderson: Tech Science Press, 2002

[19] Li S F, Liu W K. Meshfree Particle Methods. New York: Springer, 2004

[20] Zhang X, Liu Y. Meshless Methods. Beijing: Tsinghua University Press, 2004 (in Chinese)

[21] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[22] Liu G R. Meshfree Methods: Moving Beyond the Finite Element Method. 2nd ed. Boca Raton: CRC Press, 2009

[23] Chen J S, Hillman M, Chi S W. Meshfree methods progress made after 20 years. Journal of Engineering Mechanics, 2017, 143(4): 04017001

[24] Bui H H, Fukagawa R, Sako K, Wells J C. Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics (SPH). Geotechnique, 2011, 61(7): 565–574

[25] Pastor M, Blanc T, Haddad B, Petrone S, Sanchez Morles M, Drempetic V, Issler D, Crosta G B, Cascini L, Sorbino G, Cuomo S. Application of a SPH depth-integrated model to landslide run-out analysis. Landslides, 2014, 11(5): 793–812

[26] Hu M, Liu M B, Xie M W, Liu G R. Three-dimensional run-out analysis and prediction of flow-like landslides using smoothed particle hydrodynamics. Environmental Earth Sciences, 2015, 73(4): 1629–1640

[27] Dai Z, Huang Y. A three-dimensional model for flow slides in municipal solid waste landfills using smoothed particle hydrodynamics. Environmental Earth Sciences, 2016, 75(2): 132

[28] Rabczuk T, Areias P M A. A new approach for modelling slip lines in geological materials with cohesive models. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(11): 1159–1172

[29] Zheng W, Zhuang X, Tannant D, Cai Y, Nunoo S. Unified continuum/discontinuum modeling framework for slope stability assessment. Engineering Geology, 2014, 179: 90–101

[30] Liu G, Zhuang X, Cui Z. Three-dimensional slope stability analysis using independent cover based numerical manifold and vector method. Engineering Geology, 2017, 225: 83–95

[31] Dolbow J, Belytschko T. Numerical integration of the Galerkin weak form in meshfree methods. Computational Mechanics, 1999, 23(3): 219–230

[32] Chen J S, Hillman M, Rüter M. An arbitrary order variationally consistent integration for Galerkin meshfree methods. International Journal for Numerical Methods in Engineering, 2013, 95(5): 387–418

[33] Duan Q, Gao X, Wang B, Li X, Zhang H, Belytschko T, Shao Y. Consistent element free Galerkin method. International Journal for Numerical Methods in Engineering, 2014, 99(2): 79–101

[34] Hematiyan M R, Khosravifard A, Liu G R. A background decomposition method for domain integration in weak-form meshfree methods. Computers & Structures, 2014, 142: 64–78

[35] Joldes G R, Wittek A, Miller K. Adaptive numerical integration in element-free Galerkin methods for elliptic boundary value problems. Engineering Analysis with Boundary Elements, 2015, 51: 52–63

[36] Wang D, Wu J. An efficient nesting sub-domain gradient smoothing integration algorithm with quadratic exactness for Galerkin meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2016, 298: 485–519

[37] Wei H, Chen J S, Hillman M. A stabilized nodally integrated meshfree formulation for fully coupled hydro-mechanical analysis of fluid-saturated porous media. Computers & Fluids, 2016, 141: 105–115

[38] Wu C T, Chi S W, Koishi M, Wu Y. Strain gradient stabilization with dual stress points for the meshfree nodal integration method in inelastic analyses. International Journal for Numerical Methods in Engineering, 2016, 107(1): 3–30

[39] Wu J, Deng J, Wang J, Wang D. A review of numerical integration approaches for Galerkin meshfree methods. Chinese Journal of Solid Mechanics, 2016, 37: 208–233 (in Chinese)

[40] Beissel S, Belytschko T. Nodal integration of the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1–4): 49–74

[41] Chen J S, Wu C T, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering, 2001, 50(2): 435–466

[42] Chen J S, Yoon S P, Wu C T. Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering, 2002, 53(12): 2587–2615

[43] Kwok O L A, Guan P C, Cheng W P, Sun C T. Semi-Lagrangian reproducing kernel particle method for slope stability analysis and post-failure simulation. KSCE Journal of Civil Engineering, 2015, 19(1): 107–115

[44] Guan P C, Chen J S, Wu Y, Teng H, Gaidos J, Hofstetter K, Alsaleh M. Semi-Lagrangian reproducing kernel formulation and application to modeling earth moving operations. Mechanics of Materials, 2009, 41(6): 670–683

[45] Lian Y, Zhang X, Liu Y. Coupling between finite element method and material point method for problems with extreme deformation. Theoretical and Applied Mechanics Letters, 2012, 2(2): 021003

[46] Zhang X, Krabbenhoft K, Sheng D, Li W. Numerical simulation of a flow-like landslide using the particle finite element method. Computational Mechanics, 2015, 55(1): 167–177

[47] Belytschko T, Bažant Z P, Yul-Woong H, Ta-Peng C. Strain-softening materials and finite-element solutions. Computers & Structures, 1986, 23(2): 163–180

[48] Chen J S, Wu C T, Belytschko T. Regularization of material instabilities by meshfree approximations with intrinsic length scales. International Journal for Numerical Methods in Engineering, 2000, 47(7): 1303–1322

[49] Chen J S, Zhang X, Belytschko T. An implicit gradient model by a reproducing kernel strain regularization in strain localization problems. Computer Methods in Applied Mechanics and Engineering, 2004, 193(27–29): 2827–2844

[50] Askes H, Pamin J, de Borst R. Dispersion analysis and element-free Galerkin solutions of second- and fourth-order gradient-enhanced damage models. International Journal for Numerical Methods in Engineering, 2000, 49(6): 811–832

[51] Wang D, Li Z. A two-level strain smoothing regularized meshfree approach with stabilized conforming nodal integration for elastic damage analysis. International Journal of Damage Mechanics, 2013, 22(3): 440–459

[52] Wang D, Li L, Li Z. A regularized Lagrangian meshfree method for rainfall infiltration triggered slope failure analysis. Engineering Analysis with Boundary Elements, 2014, 42: 51–59

[53] Rabczuk T, Belytschko T, Xiao S P. Stable particle methods based on Lagrangian kernels. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12‒14): 1035–1063

[54] Maxars J. Mechanical damage and fracture of concrete structures. In: Proceedings of the 5th International Conference of Fracture. Cannes, 1981, 4: 1499–1506

[55] Simo J C, Ju J W. Strain- and stress-based continuum damage models—II. Computational aspects. International Journal of Solids and Structures, 1987, 23(7): 841–869

[56] Ju J W. On energy-based coupled elastoplastic damage theories: Constitutive modeling and computational aspects. International Journal of Solids and Structures, 1989, 25(7): 803–833

[57] Simo J C, Hughes T J R. Computational Inelasticity. New York: Springer, 1998

[58] Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils. New York: John Wiley & Sons , 1993

[59] Song X, Borja R I. Mathematical framework for unsaturated flow in the finite deformation range. International Journal for Numerical Methods in Engineering, 2014, 97(9): 658–682

[60] Cho S E, Lee S R. Instability of unsaturated soil slopes due to infiltration. Computers and Geotechnics, 2001, 28(3): 185–208

[61] Borja R I, White J A. Continuum deformation and stability analyses of a steep hillside slope under rainfall infiltration. Acta Geotechnica, 2010, 5(1): 1–14

[62] Jacquard C. Experimental study in laboratory of a capillary barrier. Dissertation for the Doctoral Degree. Paris: Ecole Mines Paris, 1988 (in French)

[63] Bourgeois M. The concept of capillary barrier: study by numerical model. Dissertation for the Doctoral Degree. Paris: Ecole Mines Paris, 1986 (in French)

[64] Wei P, Xiao W. Area calculation of three dimensional polygon. Chinese Mathematics Bulletin, 1984, 2: 18–21 (in Chinese)

[65] Wang D, Xie P, Lu H. Meshfree consolidation analysis of saturated porous media with stabilized conforming nodal integration formulation. Interaction and Multiscale Mechanics, 2013, 6(2): 107–125

[66] Chi S W, Siriaksorn T, Lin S P. Von Neumann stability analysis of the u-p reproducing kernel formulation for saturated porous media. Computational Mechanics, 2017, 59(2): 335–357

[67] Kawamura S, Miura S, Ishikawa T, Yokohama S. Rainfall-induced failure of unsaturated volcanic slope subjected to freeze-thaw action and its mechanism. JSCE Journal of Geotechnical and Geoenvironmental Engineering, 2010, 66(3): 577–594

[68] Li W C, Li H J, Dai F C, Lee L M. Discrete element modeling of a rainfall-induced flowslide. Engineering Geology, 2012, 149–150: 22–34

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature