3D mode discrete element method with the elastoplastic model

Abstract

The three-dimensional mode-deformable discrete element method (3MDEM) is an extended distinct element approach under the assumptions of small strain, finite displacement, and finite rotation of blocks. The deformation of blocks is expressed by the combination of the deformation modes in 3MDEM. In this paper, the elastoplastic constitutive relationship of blocks is implemented on the 3MDEM platform to simulate the integrated process from elasticity to plasticity and finally to fracture. To overcome the shortcomings of the conventional criterion for contact fracturing, a new criterion based on plastic strain is introduced. This approach is verified by two numerical examples. Finally, a cantilever beam is simulated as a comprehensive case study, which went through elastic, elastoplastic, and discontinuous fracture stages.

Keywords

mode discrete element method, elastoplastic, numerical method, discontinuum, contact

Introduction

The discrete element method (DEM) is a discontinuum analysis method that was first proposed by Cundall [1] to study two-dimensional slope stability problems in jointed rock masses. In DEM, bodies may be rigid or deformable. To calculate the deformation of these bodies, DEM was extended [2-5] to deformable elements. The bodies were divided into internal elements to enable the calculation of deformability. However, to obtain a reasonable result for engineering problems in large-scale structures, this method may require a large number of elements in addition to small time steps. Another disadvantage of this method is that a body of complex shapes must be divided into many zones, even if only a simple deformation pattern is required.

A complex deformation pattern may be achieved by the superposition of several mode shapes for the entire body. Williams and Mustoe [6,7] rewrote the matrix equation of motion for an element in terms of a set of orthogonal modes that may or may not be two-dimensional eigenmodes. Any number of these modes may be combined to obtain the complexity required by the deformation pattern. However, it is difficult to find a group of orthogonal modes to describe the deformation of bodies. In some conditions, the coupling effect between motion and deformation cannot be neglected. A somewhat similar scheme was devised by Shi [8] in his discontinuous deformation analysis (DDA). DDA uses series approximations to supply an increasingly complex set of strain patterns superimposed for each block. However, the use of direct strain modes may be inconsistent [7]. This also applies to the “simply deformable” element of Cundall [9].

Ideas based on mode decomposition have been used to simplify deformation within a discrete element or block [6,7,10]. The three-dimensional mode-deformable discrete element method (3MDEM) was built to efficiently model bodies of complicated shapes [11]. The deformation of blocks can be expressed by the combination of a series of deformation modes that can be decoupled under given orthogonal conditions. The simulation results [12] demonstrate that the present method are reasonably accurate under small deformation conditions compared to the finite element method. Similar results can be obtained by the present method and the commercial three-dimensional deformable distinct element code (3DEC) in the elastic stage under finite deformation conditions.

The advantage of traditional continuous media mechanics methods (such as the FEM) is elastic and elastoplastic calculations for a continuum, while discontinuous media mechanics methods are mainly about the simulation from continuum to discontinuum. It is obviously important to combine the plastic and discontinuum behaviors in the simulation. In other words, the plastic strain before failure should be included in the simulation.

In this paper, the constitutive relationship is extended from the elastic to the elastoplastic stage, and a simulation from continuum to discontinuum is performed to include an elastic-plastic fracture on the proposed 3MDEM platform (Fig. 1). The elastic-plastic step is modeled by the constitutive relationship of blocks, while the solution to the plastic-fracture step is face-to-face contact failure.

3MDEM

The motion and deformation of an individual block can be decomposed into three types of motion under small deformation: block translation, rotation, and deformation (Fig. 2). This section will only describe block rotation and deformation. Translation needs no further analysis due to its simplicity. For each individual block, the density, volume, and moment of inertia are assumed to be constant for small deformations.

Generally, in dealing with multi-body systems, three types of coordinate systems are required. The global coordinate system is fixed and represents a unique standard for all blocks in the system. This coordinate system is named frame a, with a base vector of ea and coordinates x, y, and z. In addition to this frame, a body frame named frame b is assigned to each block. The base vector of frame b is eb, and the coordinates of this frame are described by X, Y, and Z. This frame is located in the gravity center of each block, and the orientation of the coordinates is chosen to be parallel to the principal axis of inertia of each block. The third frame is superimposed with frame b of each block at the beginning of each time step, and its orientation remains steady in the time step. This frame, named frame b’, is the reference frame.

Additional details of 3MDEM are presented in Ref. [12]. Only a brief introduction to the translation, rotation, and deformation equations are provided in the following section.

Translation equation

The translation of a block can be represented by the displacement of its gravity center:

(1){\displaystyle \int {\rho }^a{\ddot{\mathit{r}}}_{o}dV={\displaystyle \int \mathit{f}dV}}+{\displaystyle \int {\nabla }^b·\sigma dV}.

On the force boundary,

(2)\mathit{F}={\mathit{F}}_{\mathrm{0}}.

When the Green formula expressed in index form is applied,

(3)m^a{\ddot{r}}_{oi}={\displaystyle \int f_{i}dV}+{\displaystyle \int F_{oi}d}\mathrm{\Gamma },

where m is the mass of the block, BoldItalic is the body force that includes gravity, and BoldItalic₀ is the boundary force that includes contact force.

Rotation equation

The rotation of frame b relative to frame a could be used to denote the rotation of the block:

(4){\displaystyle \int {\mathit{R}}_{\mathrm{0}}^b\times \rho }\left[{\dot{\mathit{\omega }}}^b\times {\mathit{R}}_{\mathrm{0}}^b+{\mathit{\omega }}^b\times ({\mathit{\omega }}^b\times {\mathit{R}}_{\mathrm{0}}^b)\right]dV={\displaystyle \int {\mathit{R}}_{\mathrm{0}}^b\times {\mathit{f}}^{b}dV}+{\displaystyle \int {\mathit{R}}_{\mathrm{0}}^b\times ({\nabla }^b·{\sigma }^b)dV}.

Substituting the boundary conditions into Eq. (4) and employing the Green formula, one can obtain:

(5)\begin{array}{l}{I}_{\mathrm{1}}{\dot{\omega }}_{\mathrm{1}}-(I_{\mathrm{2}}-I_{\mathrm{3}}){\omega }_{\mathrm{2}}{\omega }_{\mathrm{3}}=M_{\mathrm{1}}\hfill \\ I_{\mathrm{2}}{\dot{\omega }}_{\mathrm{2}}-(I_{\mathrm{3}}-I_{\mathrm{1}}){\omega }_{\mathrm{3}}{\omega }_{\mathrm{1}}=M_{\mathrm{2}}\hfill \\ I_{\mathrm{3}}{\dot{\omega }}_{\mathrm{3}}-(I_{\mathrm{1}}-I_{\mathrm{2}}){\omega }_{\mathrm{1}}{\omega }_{\mathrm{2}}=M_{\mathrm{3}}\hfill \end{array},

where I₁, I₂, and I₃ are the moments of inertia of the block in frame b; and M₁, M₂, and M₃ are moments around the three coordinate axes of frame b.

Deformation equation

The deformation of a block can be represented by the sum of deformation modes:

(6){\displaystyle {\int }_V\rho {\mathrm{\Phi }}_{ik}{\mathrm{\Phi }}_{ij}{\ddot{\alpha }}_{j}dV={\displaystyle {\int }_V{\mathrm{\Phi }}_{jk}{f}_j^{b}dV}+{\displaystyle {\int }_S{\mathrm{\Phi }}_{jk}{F}_{oj}^{b}d{\mathrm{\Gamma }}_{\mathrm{2}}}}-{\displaystyle {\int }_V{\mathrm{\Phi }}_{ki,j}^T{\sigma }_{ij}^{B}dV}-{\displaystyle {\int }_V\rho {{\mathrm{\Phi }}_{jk}}^a{\ddot{r}}_{\mathrm{0}j}^{b}dV}-{\displaystyle {\int }_V\rho {\mathrm{\Phi }}_{jk}\left[{\dot{\tilde{\omega }}}_{ji}^b{R}_{\mathrm{0}i}^b\right]dV}.

Let:

(7)\begin{array}{l}M={\displaystyle {\int }_V\rho {\mathrm{\Phi }}_{ik}{\mathrm{\Phi }}_{ij}dV}\hfill \\ S_A={\displaystyle {\int }_V{\mathrm{\Phi }}_{jk}{f}_j^{b}dV}+{\displaystyle {\int }_S{\mathrm{\Phi }}_{jk}{F}_{oj}^{b}d{\mathrm{\Gamma }}_{\mathrm{2}}}\hfill \\ S_I={\displaystyle {\int }_V{\mathrm{\Phi }}_{ki,j}^T{\sigma }_{ij}^{B}dV}\hfill \\ S_T={\displaystyle {\int }_V\rho {{\mathrm{\Phi }}_{jk}}^a{\ddot{r}}_{\mathrm{0}j}^{b}dV}+{\displaystyle {\int }_V\rho {\mathrm{\Phi }}_{jk}\left[{\dot{\tilde{\omega }}}_{ji}^b{R}_{\mathrm{0}i}^b\right]dV},\hfill \end{array}

Where M is the general mass matrix, S_A is the general external force matrix, S_I is the general internal force matrix, and S_T is the couple matrix of rigid movement.

Substitute Eq. (7) into Eq. (6) to obtain:

(8)M\ddot{\alpha }=S_A-S_I-S_T.

If the deformation modes satisfy the following conditions:

(9)\begin{array}{l}{\displaystyle \int {\varphi }^N\mathrm{d}V}=\mathrm{0};\hfill \\ \begin{array}{ll}{\displaystyle \int {\varphi }^N·{\varphi }^L\mathrm{d}V}=\mathrm{0},& \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}L\ne N,\end{array}\hfill \end{array}

then all deformation modes in Eq. (9) will be decoupled from each other, and the couple vector of rigid movement will be zero. Equation (8) can be simplified to:

(10)m^L{\ddot{\alpha }}^L=S_A^L-S_I^L,

where m^L, S_A^L, and S_I^L are the L^th row diagonal element of matrix M, S_A, and S_I, respectively.

Calculation process

3MDEM is based on a dynamic (time domain) algorithm that solves the equations of motion of the block system by an explicit finite difference method[12]. At each time step, the force, moment, stress, and force boundary conditions are used as inputs of the block systems; the equations of translation, rotation, and deformation for the blocks are then solved. Integrating motion and deformation provides the new block positions and the contact displacement increments (or velocities). The contact force-displacement law and displacement relationship law are then used to obtain the new contact forces to be applied to the blocks in the next time step. The computation flow diagram of 3MDEM is illustrated in Fig. 3.

Elastoplastic constitutive model

Stress-strain relationship

Six mode coefficients make up the Cauchy strain tensor ϵ for the simple deformation mode distinct element method. The stress tensor for an elastic block is:

(11)\mathrm{\Delta }\sigma =[D_e]\mathrm{\Delta }ϵ,

where D_e is the elastic matrix.

The stress at the end of the step is shown as Eq. (12):

(12)\sigma (t+\mathrm{\Delta }t)=\sigma (t)+\mathrm{\Delta }\sigma .

When the stress in Eq. (12) is compiled in frame b, and the third equation in Eq. (7) is integrated, the general internal force S_I can then be obtained. The strain increment for the following step is calculated by Eq. (8).

If the plasticity of the blocks is considered, the stress state of the block can be represented by that of centroid when determining stress yield, and the plastic flow rule of the block is the same as the rule for the centroid in this paper. Using the von Mises criterion as an example, the yield function can be expressed by Eq. (13):

(13)F=J_{\mathrm{2}}-{\sigma }_s^{\mathrm{2}},

where J_{\mathrm{2}} is the 2nd deviatoric stress invariant and {\sigma }_s is the yield stress.

For associated flow, the plastic potential function is the same as the yield function in Eq. (13). The plastic matrix can be derived using Eq. (14), where H is the hardness parameter.

(14)[D_p]={\displaystyle \frac{{\left\{{\displaystyle \frac{\partial F}{\partial [\sigma ]}}\right\}}^{\mathrm{T}}[D_e][D_e]\left\{{\displaystyle \frac{\partial F}{\partial [\sigma ]}}\right\}}{H+{\left\{{\displaystyle \frac{\partial F}{\partial [\sigma ]}}\right\}}^{\mathrm{T}}[D_e]\left\{{\displaystyle \frac{\partial F}{\partial [\sigma ]}}\right\}}}.

Deviation of Eq. (13) is

(15)\left\{{\displaystyle \frac{\partial F}{\partial [\sigma ]}}\right\}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{J_{\mathrm{2}}}}}{\{s_x,s_y,s_z,\mathrm{2}{\tau }_{xy},\mathrm{2}{\tau }_{xz},\mathrm{2}{\tau }_{yz}\}}^{\mathrm{T}},

where s is the stress deviator and \tau is the shear stress.

Substituting Eq. (15) into Eq. (14), the plastic matrix [D_p] for the von Mises criterion is derived as Eq. (16):

(16)[D_p]={\displaystyle \frac{\mathrm{3}{G}^{\mathrm{2}}}{J_{\mathrm{2}}(\mathrm{3}G+H)}}\left[\begin{array}{cccccc}{s}_x^{\mathrm{2}}& s_x{s}_y& s_x{s}_z& s_x{\tau }_{xy}& s_x{\tau }_{xz}& s_x{\tau }_{yz}\\ & s_y^{\mathrm{2}}& s_y{s}_z& s_y{\tau }_{xy}& s_y{\tau }_{xz}& s_y{\tau }_{yz}\\ & & s_z^{\mathrm{2}}& s_z{\tau }_{xy}& s_z{\tau }_{xz}& s_z{\tau }_{yz}\\ & & & {\tau }_{xy}^{\mathrm{2}}& {\tau }_{xy}{\tau }_{xz}& {\tau }_{xy}{\tau }_{yz}\\ & & & & {\tau }_{xz}^{\mathrm{2}}& {\tau }_{xz}{\tau }_{yz}\\ & & & & & {\tau }_{yz}^{\mathrm{2}}\end{array}\right].

The stress increment for the nth step in the plasticity stage is calculated by Eq. (17):

(17)\mathrm{\Delta }{\sigma }^n=[D_{ep}^n]\mathrm{\Delta }{ϵ}^n=([D_e]-[D_p^n])\mathrm{\Delta }{ϵ}^n.

The stress at the end of the step is:

(18){\sigma }^{n+\mathrm{1}}={\sigma }^n+\mathrm{\Delta }{\sigma }^n.

Substituting Eq. (17) into Eq. (18), the general internal force increment \mathrm{\Delta }{S}_I can be obtained:

(19)\mathrm{\Delta }{S}_I={\displaystyle \int \mathrm{\Delta }{\sigma }_{ij}^b{\mathrm{\Phi }}_{ki,j}^{\mathrm{T}}\mathrm{d}V}=\left\{\begin{array}{l}{\displaystyle \int \mathrm{\Delta }{\sigma }_{XX}\mathrm{d}V}\hfill \\ {\displaystyle \int \mathrm{\Delta }{\sigma }_{YY}\mathrm{d}V}\hfill \\ {\displaystyle \int \mathrm{\Delta }{\sigma }_{ZZ}\mathrm{d}V}\hfill \\ {\displaystyle \int (\mathrm{\Delta }{\sigma }_{XY}+\mathrm{\Delta }{\sigma }_{YX})\mathrm{d}V}\hfill \\ {\displaystyle \int (\mathrm{\Delta }{\sigma }_{XZ}+\mathrm{\Delta }{\sigma }_{ZX})\mathrm{d}V}\hfill \\ {\displaystyle \int (\mathrm{\Delta }{\sigma }_{YZ}+\mathrm{\Delta }{\sigma }_{ZY})\mathrm{d}V}\hfill \end{array}\right\}.

The general internal force S_I can be easily derived from \mathrm{\Delta }{S}_I. The strain increment, including the elastic and plastic strains, can then be calculated using Eq. (10).

Loading and unloading conditions

3MDEM is mainly used for discontinuous calculations, and the blocks are associated by contact. When the contact relationship satisfies the contact failure criterion, the blocks are no longer in contact, and the continuum becomes a discontinuum. Contact failure leads to the unloading of local blocks. Therefore, determining loading or unloading is important in an elastoplastic constitutive relationship in Fig. 4

The loading and unloading condition is defined by Eq. (20):

(20)\begin{array}{ll}{\displaystyle \frac{\partial F}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}>\mathrm{0},\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g};\hfill \\ {\displaystyle \frac{\partial F}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}=\mathrm{0},\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{n}\mathrm{e}\mathrm{u}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g},\mathrm{a}\mathrm{n}\mathrm{d}\hfill \\ {\displaystyle \frac{\partial F}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}=\mathrm{0}\hfill & \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}.\hfill \end{array}

For the von Mises model, substituting Eq. (13) into Eq. (20) results in:

(21){\displaystyle \frac{\partial F}{\partial {\sigma }_{ij}}}\mathrm{d}{\sigma }_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{J}}}(s_x,s_y,s_z,\mathrm{2}{\tau }_{xy},\mathrm{2}{\tau }_{xz},\mathrm{2}{\tau }_{yz}){(\mathrm{d}{\sigma }_x,\mathrm{d}{\sigma }_y,\mathrm{d}{\sigma }_z,\mathrm{d}{\tau }_{xy},\mathrm{d}{\tau }_{xz},\mathrm{d}{\tau }_{yz})}^{\mathrm{T}}=(s_x\mathrm{d}{\sigma }_x+s_y\mathrm{d}{\sigma }_y+s_z\mathrm{d}{\sigma }_z+\mathrm{2}{\tau }_{xy}\mathrm{d}{\tau }_{xy}+\mathrm{2}{\tau }_{xz}\mathrm{d}{\tau }_{xz}+\mathrm{2}{\tau }_{yz}\mathrm{d}{\tau }_{yz})/\mathrm{2}\sqrt{J}

and this expression shows the unloading condition for the von Mises model.

Contact failure criterion

The contact failure criterion is useful for elastic blocks in 3MDEM. If the plastic constitutive relationship of the blocks is considered, their stiffness becomes much smaller. The deformation of the blocks can be large with slight changes in the stress state. Stress is constant even under the ideal elastoplastic constitutive model. Furthermore, stress is insensitive in the plastic stage, and numerical error can affect the determination of contact failure. Therefore, tension and Mohr coulomb criteria based on stress judgment have a high tendency for error.

The application of the elastoplastic constitutive relationship for blocks simulates the transition of the block system from continuum to discontinuum, which include the elastic, plastic, and fracture stages. Face-to-face contact can be assumed as a non-failure in the elastic stage because contact failure always occurs in the plastic stage. Therefore, a criterion based on plastic strain is presented in this paper. When tension or shear plastic strain of both blocks in face-to-face contact is larger than the ultimate plastic strain, contact failure occurs- as shown in Eq. (22):

(22)\begin{array}{lll}\min ({ϵ}_{xx}^i,{ϵ}_{xx}^j)>{ϵ}_{tt}& \mathrm{o}\mathrm{r}& \min (\left|{ϵ}_{xy}^i\right|,\left|{ϵ}_{xy}^j\right|)>{ϵ}_{ts}\end{array},

where ϵ_xx is the normal plastic strain of the contact block, ϵ_xy is the shear plastic strain of the contact block, i and j are the block numbers in contact, and ϵ_tt and ϵ_ts are the ultimate tension and shear plastic strain, respectively.

In the example shown in Fig. 6, if the strain of Block 1 is smaller than the ultimate strain while the strain of Block 2 and Block 3 reach the ultimate strain, the contact between Block 2 and Block 3 that satisfies the failure criterion will be broken.

Examples for the elastoplastic constitutive model

Axial tension and pure shear

For axial tension analysis, a face force of 2.0 × 10⁶ N/m² is applied on the left and right surfaces of a block with the dimensions of 1 m × 1 m × 1 m (Fig. 7). The block has a Young’s modulus of 24 GPa, a Poisson’s ratio of 0.2, and a density of 2400 kg/m³. The von Mises criterion is chosen with a yield stress of 1.732 MPa and a hardening parameter of 1 GPa. A comparative calculation by the FEM software ABAQUS is conducted using the same parameters.

The stress and strain relationship is shown in Fig. 8. When the tension stress reaches 1.732 MPa or when the von Mises stress is 1.732 MPa, the block will yield. As a result of adopting the hardening model, stress increases slowly to 2.0 MPa. When the block yields, strain increases rapidly because of plastic flow. The result of the 3MDEM is similar with that of the FEM. Therefore, the 3MDEM is considered a fairly accurate method to describe the material nonlinearity problem under tension load.

For pure shear analysis, a face force of 1.2 × 10⁶ N/m² is applied on the left, right, top, and bottom surfaces of a block with dimensions of 1 m × 1 m × 1 m (Fig. 9). The block has a Young’s modulus of 24 GPa, a Poisson’s ratio of 0.2, and a density of 2400 kg/m³ with a yield stress of 1.732 MPa and a hardening parameter of 1 GPa. A comparative calculation by the FEM software ABAQUS is conducted using the same parameters.

The stress-and-strain relationship is shown in Fig. 10. When the shearing stress becomes 1.0 MPa or when the von Mises stress is 1.732 MPa, the block will yield. As a result of adopting the hardening model, stress increases slowly to 1.2 MPa. When the block yields, strain increases rapidly because of plastic flow. The result of the 3MDEM is similar with that of the FEM. Therefore, 3MDEM is considered an accurate method to describe the material nonlinearity problem under shearing load.

Tensile and shearing facture

Two blocks are chosen for tension fracture analysis. The shaded block on the left is a rigid body fixed in all directions, and the block on the right (1 m × 1 m × 1 m) has a Young’s modulus of 24 GPa, a Poisson’s ratio of 0.2, and a density of 2400 kg/m³. The von Mises criterion is chosen with a yield stress of 1.732 MPa and a hardening parameter of 1 GPa. The ultimate tension strain is 1000 × 10^-6. The normal contact stiffness is 50 GPa, and the shear contact stiffness is 0. The damping constant is 0.2. A uniform load of 3000 kN/m² is applied on the right surface (Fig. 11).

The displacement-, strain-, and stress-time interval curves in the centroid are shown in Fig. 12. At 0.09 s, the displacement suddenly increases in a parabolic curve because the block moved with constant acceleration after contact failure. At the same time, plastic strain is unchanged as the stress unloads.

The strain-time interval curve shows that the block yields at 0.06 s, and the strain is unchanged after reaching 1000 × 10^-6 because the block recovered elastically after face-to-face contact failure and unloading. The stress-time interval curve shows that the maximum stress occurs at approximately 1.5 MPa, which is half of the uniform load. The stress at 1.5 MPa represents the stress state of the centroid. When the block is cut along the mid-vertical face, the left and right blocks each bear half of the inertia force. Therefore, the stress of the centroid in the acceleration mode is half of that in the equilibrium state.

Two blocks with the same sizes and material parameters as the blocks for the tension fracture analysis are chosen for the shear fracture analysis. The normal contact stiffness is 50 GPa, and the shear contact stiffness is 20.8 GPa. A vertical uniform load of 1600 kN/m² is applied to the right end (Fig. 13).

The displacement-, strain-, and stress-time interval curves in the centroid are shown in Fig. 14. At 0.12 s, the displacement suddenly increases because shearing failure of face-to-face contact occurs. At the same time, the plastic strain is unchanged as the stress unloads. Compared with tension failure, shearing failure involves a more complex process. After the shearing fracture, the right block rotates and impacts the fixed block.

The strain-time interval curve shows that the block yields at 0.06 s. The strain increases to over 2000 × 10^-6 at 0.13 s and remains constant from then because the block separated from the fixed block and recovered elastically after contact failure.

Failure process of cantilever beam

To analyze the failure process of a cantilever beam, a cantilever beam (300 m × 60 m × 15 m) with a fixed left end (Fig. 15) with a Young’s modulus of 24 GPa, a Poisson’s ratio of 0.2, and a density of 2400 kg/m³ is chosen. The von Mises criterion is chosen with a yield stress of 1.732 MPa and a hardening parameter of 1 GPa. The ultimate tension and shearing strain are 1000 × 10^-6 and 2000 × 10^-6, respectively.

Four layers of blocks are meshed in the y direction (Fig. 16), and each block is 15 m × 15 m × 15 m. The normal contact stiffness is 5 GPa, and the shear contact stiffness is 2.08 GPa. The damping constant is 0.5. A uniform load is applied to the top surface, which increases linearly from 0 N/m² to 0.35 × 10⁶ N/m² in 30 s (Fig. 17).

A comparison of the yield zones from the FEM and the 3MDEM is shown in Fig. 18. When the load is 0.5 P, calculations by both the 3MDEM and the FEM are in a continuum, and the yield zones are similar. However, the calculations under continuous loading show different results. The yield zone from the FEM does not expand, while the plastic strain at the fixed end increases rapidly until the strain no longer converges. The result of the 3MDEM calculations shows that when the load is 0.6 P, the strain of Block 4 at the fixed end reaches its ultimate strain, and face-to-face contact failure occurs. This is the stage of transformation from continuum to discontinuum.

The stress determined from the FEM and the 3MDEM is shown in Fig. 19; the stress obtained using the FEM is shown as a solid line, and that from the 3MDEM is shown as a dashed line. When the load is 0.5 P, the dashed line agrees well with the solid line, except at the end of stress concentration. The stress concentration at the fixed end is not fully expressed because the six first-order modes are not sufficient to describe the deformation. When the load reaches 0.6 P, the calculation of the 3MDEM becomes continuous. Block 4 at the left end breaks from the fixed block, and local unloading occurs (Fig. 20). The stress decreases to zero in the top-left area and redistributed in the nearby zone.

The failure process of the cantilever beam is shown in Fig. 21, where the displacement is scaled 10 times. Before the load reaches 0.6 P, the 3MDEM and FEM results are very similar. The blocks at the fixed end yield from the start, and a large bending deformation occurs. When the load reaches 0.6 P, contact failure happens in the top-left zone in the 3MDEM calculation, and the fracture of the interface gradually expands downward. The cantilever beam is separated from the fixed block due to the interface damage resulting from tensile-shear failure. The result based on the FEM shows that the plastic strain at the fixed end increases continuously, and the deformation becomes wider until the strain no longer converges. The discontinuum after fracture cannot be simulated using the FEM.

The example on the cantilever beam shows that the result based on the 3MDEM is similar to the FEM in the continuous stage and therefore reasonable. The simulation from the continuum to discontinuum stage and from the elastic to elastoplastic to fracture stage is realized.

Conclusions

The 3MDEM is an effective numerical method to simulate mechanical behaviors while coupling inelasticity, finite motion and fracturing. In this paper, an elastoplastic constitutive relation is implemented on the 3MDEM platform. A complete simulation can be performed, and it can include the integrated process from continuum to discontinuum or from elasticity to plasticity and then to fracturing. The elastic-to-plastic stage is realized by the elastoplastic constitutive relation of blocks, while the plastic-to-fracture stage is realized by the face-to-face contact failure. To overcome the shortcomings of the conventional contact criterion based on force and stress, a new criterion based on plastic strain is presented. The effectiveness and accuracy of this method are verified with several numerical examples.

References

[1] Cundall P A. Computer model for simulating progressive large scale movements in blocky systems. Proc. Symp. Int. Soci. Rock Mech., Nancy, France, Vol. 1, paper No II-8, 1971

[2] Cundall P A. Formulation of a three-dimensional distinct element model—Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1988, 25(3): 107-116

[3] Hart R, Cundall P A, Lemos J. Formulation of a three-dimensional distinct element model—Part II. Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1988, 25(3): 117-125

[4] Cundall P A. UDEC-A generalized discrete element program of modeling jointed rock. European Research Office, U.S. Army, 1980

[5] Itasca, 3DEC- 3 Dimensional discrete element code, Version 3.0, User’s manual, Itasca Consulting Group, Inc. USA. 1996

[6] Williams J R, Hocking G, Mustoe G G W. The theoretical basis of the discrete element method, Proceedings of the NUMETA conference, Swansea, 1985, 7-11

[7] Williams J R, Mustoe G G W. Modal methods for the analysis of discrete systems. Computers and Geotechnics, 1987, 4(1): 1-19

[8] Shi G H. Discontinuous Deformation Analysis- A New Numerical Model for the Statics and Dynamics of Block Systems, Lawrence Berkeley Laboratory, Report to DOE OWTD, ContractAC03-76SF0098, September 1988; also DissertationTip, University of California, Berkeley, <month>August</month> 1989

[9] Cundall P A, Marti J, Beresford P J, Last N C, Asgian M I. Computer Modeling of Jointed Rock Masses, U.S. Army Engineer Waterways Experiment Station, Vicksburg, Mississippi, Tech. Report N-78-4, <month>August</month> 1978

[10] Genhua S. Discontinuous deformation analysis: A new numerical model for the statics and dynamics of deformable block structures, Engineering Computations, 1992. 9

[11] Chong Z, Feng J, Yanli H. 3-D simple deformable distinct element method. Chinese Journal of Geotechnical Engineering, 2007, 29(4): 6-10 (in Chinese)

[12] Jin F, Zhang C, Hu W, Wang J. 3-D mode discrete element method: elastic model. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(1): 59-66

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg