Investigation of a boundary simulation of continuity using the discrete solid element method

Abstract

The discrete solid element method is an efficient numerical method that simulates the large deformation, strong material nonlinearity, fracture, and dynamic problems of continuity. In the discrete solid element method model, the spring stiffness of the spherical elements on the boundary is different from that inside the discrete solid element method model based on the principle of conservation of energy. The spring stiffness of the spherical elements on the boundary of the discrete solid element method model is shown to have a significant effect on the macroscopic properties. According to the position of the spherical elements on the boundary of the discrete solid element method model, the spherical elements on the boundary are divided into three types, which are spherical elements on the surface position, on the edge position, and on the corner position. To accurately reflect the mechanical behavior of the material, the principle of energy conservation is used to strictly deduce the spring stiffness of the three types of spherical elements on the boundary, and the relationship between the spring stiffness and elastic constants is established. The numerical example shows that the calculation accuracy of the discrete solid element method in modeling the mechanical behavior of continuity is improved after the spring stiffness of the spherical elements on the boundary is revised. In addition, the applications of the discrete solid element method to dynamic buckling of the thin plate and buckling of the cracked thin plate are also given.

Keywords

Spring stiffness spherical elements on the boundary discrete solid element method strong nonlinear problem

Introduction

In the field of numerical calculation, the finite element method (FEM) is widely used to solve various mechanical and nonlinear problems due to its perfect theory and mature commercial software. The basis of the FEM is the variational principle and continuum mechanics. In terms of unit deformation solution, the FEM must satisfy the deformation coordination condition. Therefore, it’s difficult to simulate the strong nonlinear mechanical behavior of materials using the FEM, which requires special processing and complicated correction. Khoei et al.¹ used the arbitrary Lagrangian–Eulerian technique in the FEM to solve large deformation problems of solid mechanic. Lynn and Isobe² develop a FEM code that applied the adaptively shifted integration (ASI)–Guess technique to cope with strong non-linearity and discontinuity problems.

In order to overcome the shortcomings of the continuum mechanics method, the meshless method has been fully developed. The discrete element method (DEM) proposed by Cundall and Strack³ is one of the representatives of the meshless method. In the DEM calculation process, the damping is used to dissipate energy so that the calculation converges. Qi and Ji-Hong⁴ studied the calculation of damping force and the selection of damping coefficient. The Nonlinear viscous damping is used to calculate static problems. While for the dynamic problem, the selected damping should be able to reproduce the real response of the natural system under the dynamic load and the Rayleigh damping is used to calculate the dynamic problem.

In the DEM, the material is discretized into mass particles with the same or different sizes. The particles can move, touch, and collide in the analysis process of the DEM without a need to deliberately meet the deformation coordination and continuous condition. Therefore, the DEM is especially suitable for solving the large deformations and the strong nonlinear problems compared with the FEM. In recent years, DEM has been an effective numerical tool to solve the dynamic and fracture problems of the continuum.

Cheng et al.⁵ introduced a two-dimensional (2D) discrete element model of hexagonal arrangement to solve the impact problem. A DEM model is proposed for the large deformation and strong nonlinear analysis of member structures by Ye and Qi,⁶ and introduces the establishment and solution of the equations of motion and contact constitutive equations. Hentz et al.⁷ used the DEM to model concrete subjected to dynamic loading at high strain rates. The cohesive beam model in the DEM was chosen by André et al.⁸ and Meguro and Tagel-Din⁹ to simulate the static and dynamic mechanical behavior of continuous materials. Griffiths and Mustoe¹⁰ used a grillage of structural elements derived from the DEM concepts to simulate an elastic continuum. Sinaie¹¹ and Hentz et al.¹² applied the DEM to solve the fracture problem of concrete structures. These applications show that the DEM has been employed successfully in solving the elastoplastic, fracture, and dynamic and strong nonlinear problems of continuity.

FEM is not ideal for the simulation of large deformation problems, while the DEM is suitable for dealing with material damage and large deformation problems. For many problems in the structure, the region where the large deformation occurs is only a part of the material, and the material deformation is small in most regions. Therefore, the DEM simulation is applied to local and large deformation regions, and the other small deformation regions are simulated by FEM. The advantages of the two calculation methods are comprehensively combined, and the combined finite DEM is established. Munjiza et al.¹³ and Rousseau et al.¹⁴ developed a combined finite discrete element model for failure and collapse of reinforced concrete structure. Lei and Zhang¹⁵ used an algorithm combining three-dimensional (3D) discrete and FEMs to simulate the impact fracture behavior of a laminated glass plate. Azevedo and Lemos¹⁶ presented a DEM/FEM coupling algorithm which enables the use of rigid circular particles only on the fracture processes, the remaining part of the structure being modeled with finite elements.

For the nature of the random distribution of particles in soil and rocks, the random arrangement of spherical elements is suitable for simulating granular material. The DEM with the random arrangement of spherical elements was employed for simulating nonlinear mechanical behaviors, large deformations, strain softening, and dynamic fracture of soil grains and rock by Li et al.,¹⁷ Bono and McDowell,¹⁸ and McDowell.¹⁹ For the continuum problems addressed by the DEM, the regular arrangement of spherical elements is more suitable in terms of accuracy and efficiency compared with the random arrangement of spherical elements. The discrete solid element method (DSEM) with the cubic arrangement of spherical elements is proposed to calculate the extremely large deformation and strong material nonlinearity of continuity by Zhu and Feng.²⁰ The precision of the discrete gird is determined by the radius of the spherical elements. For the random arrangement, the radius of the spherical elements is determined by the size of the actual particles in the granular material. For the regular arrangement, the radius of the spherical elements is determined by the geometric size of the continuum structure.

In the DSEM, eight spherical elements constitute a DSEM cubic model. When the continuity is simulated by the cubic arrangement of spherical elements, the number of spherical elements connected to the spherical elements on the boundary of the DSEM model is different from that connected to the spherical elements inside the DSEM model, resulting in the spring stiffness of the spherical elements on the boundary being different from that inside the calculation model based on the principle of conservation of energy. If all the springs adopt the spring stiffness of the spherical elements inside the DSEM model, the displacement and the contact force of the spherical elements on the boundary will produce large errors. Especially in the case of controlling the number of spherical elements, when the DSEM is used to simulate the mechanical behavior of thin-walled structures, the displacement and the contact force of the spherical elements on the boundary may produce erroneous results.

The main purpose of this article is to derive the spring stiffness of the spherical elements on the boundary of the DSEM model using the principle of conservation of energy and to establish the relationship between the spring stiffness and elastic constants. Through revising the spring stiffness of the spherical elements on the boundary, the displacement and the contact force errors of the spherical elements on the boundary are reduced, and the accuracy of the simulation of the mechanical behavior of continuity obtained using the DSEM is improved. First, the physical model and the calculation principle of the DSEM are introduced, including the motion equation, elastoplastic force–displacement equations, and yield equation. Second, the spherical elements on the boundary are classified according to their positions, including the spherical elements on the surface position, on the edge position, and on the corner position. The principle of energy conservation is used to derive the spring stiffness of various spherical elements on the boundary. Finally, the spring stiffness of the spherical element on the boundary is verified by the thin plate example, with prominent boundary problems. Moreover, the ability of simulating strong nonlinear mechanical behaviors using the DSEM is demonstrated by the dynamic buckling of the thin plate and the buckling of the cracked thin plate.

DSEM

Physical model

In the DSEM, the material is discretized into mass spherical elements. As shown in Figure 1, eight spherical elements compose the basic cube model in the DSEM. The two spherical elements that are located on the edge and the diagonal line of the DSEM basic cube model are connected by springs. Therefore, there are two different groups of springs in the DSEM model, which are the springs on the edge and the diagonal line. Each group of springs includes one normal spring $(kn)$ and two shear springs ( $k s_{1}$ and $k s_{2}$ ). The discretization scheme is formed by the spherical elements, and the springs can effectively simulate the complicated behavior of continuity, which includes extremely large deformation, high material nonlinearity, and fracture.

Calculation principle

The main calculation flowchart of the DSEM is shown in Figure 2(a). The motion equation (Newton’s second law of motion) is applied to calculate the velocity and displacement of the spherical elements, and the force–displacement equation is applied to calculate the contact force (spring force) between spherical elements. The calculation flow is show in Figure 2(b). Given the relative displacement between spherical elements (obtained from the previous time step), the contact force between spherical elements is updated by the force–displacement equation. Then, according to the contact force and the motion equation, the velocity and displacement of the spherical elements are updated.

The motion equation, as presented in equation (1), is solved by using the explicit central finite difference scheme. Then, the displacement of the spherical element i at $t_{0} + Δ t$ can be obtained as

\sum F_{i} = m_{i} {\overset{\cdot\cdot}{u}}_{i}

(1)

u_{i} (t_{0} + Δ t) = u_{i} (t_{0}) + ({\overset{\cdot}{u}}_{i} (t_{0} - \frac{Δ t}{2}) + \frac{\sum F_{i}^{t_{0}}}{m_{i}} \cdot Δ t) \cdot Δ t

(2)

where $m_{i}$ is the spherical element mass, ${\overset{\cdot\cdot}{u}}_{i}$ is the spherical element acceleration, $\sum F_{i}$ is the sum of forces acting on the spherical element i, $u_{i} (t_{0} + Δ t)$ is the spherical element displacement at $t_{0} + Δ t$ , $u_{i} (t_{0})$ is the spherical element displacement at initial moment $t_{0}$ , ${\overset{\cdot}{u}}_{i} (t_{0} - (Δ t / 2))$ is the spherical element velocity at $t_{0} - (Δ t / 2)$ , and $Δ t$ is the time step.

In the DSEM, the relationship between the contact force increment and the relative displacement between spherical elements is described by the force–displacement equation. The elastoplastic force–displacement equations for a perfect plastic material in the DSEM are presented in equation (3)

\begin{matrix} [\begin{matrix} Δ F n_{1} \\ Δ {Fs}_{1, 1} \\ Δ {Fs}_{1, 2} \\ Δ F n_{2} \\ Δ {Fs}_{2, 1} \\ Δ {Fs}_{2, 2} \end{matrix}] = [\begin{matrix} k n_{1} \\ k s_{1, 1} \\ k s_{1, 2} \\ k n_{2} \\ k s_{2, 1} \\ k s_{2, 2} \end{matrix}] \\ ([\begin{matrix} Δ u_{n 1} \\ Δ u_{s 1, 1} \\ Δ u_{s 1, 2} \\ Δ u_{n 2} \\ Δ u_{s 2, 1} \\ Δ u_{s 2, 2} \end{matrix}] - d λ [\begin{matrix} \frac{\partial f}{\partial F n_{1}} \\ \frac{\partial f}{\partial F s_{1, 1}} \\ \frac{\partial f}{\partial F s_{1, 2}} \\ \frac{\partial f}{\partial F n_{2}} \\ \frac{\partial f}{\partial F s_{2, 1}} \\ \frac{\partial f}{\partial F s_{2, 2}} \end{matrix}]) \end{matrix}

(3)

where $k n_{1}$ , $k s_{1, 1}$ , and ${ks}_{1, 2}$ are the normal and shear spring stiffness on the edge of the cube model; $Δ F n_{1}$ , $Δ F s_{1, 1}$ , and $Δ F s_{1, 2}$ are the corresponding normal and shear contact force increments on the edge; $Δ u_{n 1}$ , $Δ u_{s 1, 1}$ , and $Δ u_{s 1, 2}$ are the corresponding normal and shear relative displacements on the edge; $k n_{2}$ , $k s_{2, 1}$ , and $k s_{2, 2}$ are the normal and shear spring stiffness on the diagonal line of the cube model; $Δ F n_{2}$ , $Δ {Fs}_{2, 1}$ , and $Δ {Fs}_{2, 2}$ are the corresponding normal and shear contact force increments on the diagonal line; $Δ u_{n 2}$ , $Δ u_{s 2, 1}$ , and $Δ u_{s 2, 2}$ are the corresponding normal and shear relative displacements on the diagonal line and $d λ$ is the non-negative scale factor.

f is the yield equation in the DSEM, as presented in equation (4)

\begin{matrix} f = \frac{{(F n_{1})}^{2} λ_{n, 1}}{k n_{1}} + \frac{{({Fs}_{1, 1})}^{2} λ_{s}}{{ks}_{1, 1}} + \frac{{({Fs}_{1, 2})}^{2} λ_{s}}{{ks}_{1, 2}} \\ + \frac{2 {(F n_{2})}^{2} λ_{n, 2}}{k n_{2}} + \frac{2 {({Fs}_{2, 1})}^{2} λ_{s}}{{ks}_{2, 1}} + \frac{2 {({Fs}_{2, 2})}^{2} λ_{s}}{{ks}_{2, 2}} - \frac{4 V σ_{s}^{2}}{36 G} \end{matrix}

(4)

where $λ_{n, 1}$ and $λ_{n, 2}$ are the coefficients of the distortion energy density corresponding to the normal spring on the edge and the diagonal line of the cube model, respectively; $λ_{s}$ is the coefficient of the distortion energy density corresponding to the shear spring both on the edge and the diagonal line; V is the volume of the zone of influence of a spherical element; $σ_{s}$ is the yield stress; and G is the shear modulus.

Spherical elements on the boundaryof the DSEM model

The spring stiffness of the spherical elements on the boundary has an important influence on the mechanical behavior of continuity modeled using the DSEM. As shown in Figure 3, according to the position of the spherical elements on the boundary of the DSEM model, the spherical elements on the boundary are divided into three types, which are the spherical elements on the surface position, on the edge position, and on the corner position. The number of the surrounding spherical elements connected to the spherical elements inside the DSEM model is different from that connected to the spherical elements on the boundary. The spherical elements inside the DSEM model are in contact with 18 surrounding spherical elements, as shown in Figure 4. The spherical elements on the surface of the boundary are in contact with 13 surrounding spherical elements, as shown in Figure 5. The spherical elements on the edge of the boundary are in contact with 9 surrounding spherical elements, as shown in Figure 6. The spherical elements on the corner of the boundary are in contact with 6 surrounding spherical elements, as shown in Figure 7. Therefore, the spring stiffness deduced from 18 spherical elements cannot be applied to spherical elements on the model boundary. The spring stiffness of the spherical elements on the boundary is different from that inside in the DSEM model based on the principle of conservation of energy. When the material is discretized into regularly arranged spherical elements in the DSEM, the spring stiffness of spherical elements on the boundary cannot adopt the spring stiffness of spherical elements inside the DSEM model. The spring stiffness of spherical elements on the boundary has an important influence on the macroscopic properties, the spring stiffness of spherical elements on the boundary need to be deduced to accurately simulate the mechanical behavior of the continuum boundary.

Spring stiffness deduction of spherical elements on the surface

In the DSEM, the spherical elements and springs form a network system that captures the complicated mechanical behavior of continuity. To accurately reflect the mechanical response of the material, the spring stiffness in the force–displacement equation is critical. The relationship between the spring stiffness and the elastic Modulus E and Poisson’s ratio $ν$ is derived by the energy conservation principle. For the system consisting of spherical elements and springs, the strain energy stored in the material can be expressed as the elastic potential energy of the normal and shear springs. As shown in Figure 8, the spherical elements on the surface of boundary are in contact with 13 surrounding spherical elements. Therefore, the total strain energy stored per unit volume of spherical elements on the surface of boundary, that is, the strain energy density, can be written as follows

Π_{i}^{″} = \frac{1}{V} \sum_{j = 1}^{13} \frac{1}{2} [\frac{1}{2} kn {(u_{n, ij})}^{2} + \frac{1}{2} k s_{1} {(u_{s 1, ij})}^{2} + \frac{1}{2} k s_{2} {(u_{s 2, ij})}^{2}]

(5)

where $u_{n, ij}$ , $u_{s 1, ij}$ , and $u_{s 2, ij}$ are the normal and shear relative displacements between the spherical element i on the surface and the surrounding spherical element j.

The relative displacements between spherical elements can be expressed by the strain and the central distance, as presented in equation (6)

u_{n, ij} = l_{0} \frac{\partial u_{n, ij}}{\partial x_{l}}; u_{s 1, ij} = l_{0} \frac{\partial u_{s 1, ij}}{\partial x_{l}}; u_{s 2, ij} = l_{0} \frac{\partial u_{s 2, ij}}{\partial x_{l}}

(6)

where $l_{0}$ is the central distance between spherical elements and $\partial u_{n, ij} / \partial x_{l}$ , $\partial u_{s 1, ij} / \partial x_{l}$ , and $\partial u_{s 2, ij} / \partial x_{l}$ are line strains along the normal and shear detection in the local coordinate system.

The transformation matrix is given by

{\begin{matrix} n \\ s_{1} \\ s_{2} \end{matrix}} = [\begin{matrix} l_{1} l_{2} m_{2} m_{1} l_{2} \\ - l_{1} m_{2} l_{2} - m_{1} m_{2} \\ - m_{1} 0 l_{1} \end{matrix}] {\begin{matrix} x \\ y \\ z \end{matrix}}

(7)

where $l_{1} = \cos α$ , $l_{2} = \cos β$ , $m_{1} = \sin α$ , $m_{2} = \sin β$ , $α$ is the rotation angle between the projection of local coordinate axis n in the global coordinate system $x - y$ plane and the positive direction of global coordinate axis x, and $β$ is the rotation angle between the local coordinate axis n and the global coordinate system $x - y$ plane.

Thus, the strain in the local coordinate system can be written in terms of global strain quantities by transformation equations, as presented in equation (8)

\begin{matrix} \frac{\partial u_{n, ij}}{\partial x_{l}} = l_{1}^{2} l_{2}^{2} ε_{x} + m_{2}^{2} ε_{y} + l_{2}^{2} m_{1}^{2} ε_{z} + l_{1} l_{2} m_{2} γ_{xy} + l_{1} l_{2}^{2} m_{1} γ_{xz} \\ + l_{2} m_{1} m_{2} γ_{yz} \\ \frac{\partial u_{s 1, ij}}{\partial x_{l}} = - l_{1}^{2} l_{2} m_{2} ε_{x} + l_{2} m_{2} ε_{y} - l_{2} m_{1}^{2} m_{2} ε_{z} \\ - (l_{1} m_{2}^{2} - l_{1} l_{2}^{2}) \frac{γ_{xy}}{2} - l_{1} l_{2} m_{1} m_{2} \frac{γ_{xz}}{2} \\ + (l_{2}^{2} m_{1} - m_{1} m_{2}^{2}) \frac{γ_{yz}}{2} \\ \frac{\partial u_{s 2, ij}}{\partial x_{l}} = - l_{1} l_{2} m_{1} ε_{x} + l_{1} l_{2} m_{1} ε_{z} - m_{1} m_{2} \frac{γ_{xy}}{2} \\ - (l_{2} m_{1}^{2} - l_{1}^{2} l_{2}) \frac{γ_{xz}}{2} + l_{2} m_{2} \frac{γ_{yz}}{2} \end{matrix}

(8)

where $ε_{x}$ , $ε_{y}$ , $ε_{z}$ , $γ_{xy}$ , $γ_{xz}$ , and $γ_{yz}$ are the global strain components.

Substituting equations (6) and (9) into equation (5) yields

Π_{i}^{″} = \frac{1}{4 V} \sum_{n = 1}^{13} [\begin{matrix} {knl}_{0}^{2} {(l_{1}^{2} l_{2}^{2} ε_{x} + m_{2}^{2} ε_{y} + l_{2}^{2} m_{1}^{2} ε_{z} + l_{1} l_{2} m_{2} γ_{xy} + l_{1} l_{2}^{2} m_{1} γ_{xz} + l_{2} m_{1} m_{2} γ_{yz})}^{2} \\ + {ksl}_{0}^{2} {(- l_{1}^{2} l_{2} m_{2} ε_{x} + l_{2} m_{2} ε_{y} - l_{2} m_{1}^{2} m_{2} ε_{z} - (l_{1} m_{2}^{2} - l_{1} l_{2}^{2}) \frac{γ_{xy}}{2} - l_{1} l_{2} m_{1} m_{2} \frac{γ_{xz}}{2} + (l_{2}^{2} m_{1} - m_{1} m_{2}^{2}) \frac{γ_{yz}}{2})}^{2} \\ + {ksl}_{0}^{2} {(- l_{1} l_{2} m_{1} ε_{x} + l_{1} l_{2} m_{1} ε_{z} - m_{1} m_{2} \frac{γ_{xy}}{2} ε_{y} - (l_{2} m_{1}^{2} - l_{1}^{2} l_{2}) \frac{γ_{xz}}{2} + l_{2} m_{2} \frac{γ_{yz}}{2})}^{2} \end{matrix}]

(9)

Assume that the radius of a spherical element is r and that the initial central distances of two spherical elements on the edge and the diagonal line of the model are $l_{0} = 2 r$ and $l_{0} = 2 \sqrt{2} r$ . As shown in Figure 8, there are five surrounding spherical elements with the initial central distance of $l_{0} = 2 r$ from the spherical element on the surface, and they are connected by springs on the edge. There are eight surrounding spherical elements with the initial central distance of $l_{0} = 2 \sqrt{2} r$ from the spherical element on the surface, and they are connected by springs on the diagonal line. Therefore, the contact force on one spherical element on the surface includes five groups of contact forces on the edge of the cube model and eight groups of contact forces on the diagonal line of the cube model, as shown in Figure 9. Each group of contact forces includes one normal and two shear contact forces. Substituting $k n_{1}$ , $k s_{1, 1}$ , ${ks}_{1, 2}$ , $k n_{2}$ , $k s_{2, 1}$ , $k s_{2, 2}$ , and $l_{0}$ into equation (9), the strain energy density in the DSEM can be expressed as

\begin{matrix} Π_{i}^{″} = ε_{x}^{2} (\frac{3 k n_{2}}{8 r} + \frac{k n_{1}}{4 r} + \frac{k s_{2, 1}}{8 r} + \frac{k s_{2, 2}}{4 r}) \\ + ε_{y}^{2} (\frac{k n_{2}}{4 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{4 r}) \\ + ε_{z}^{2} (\frac{3 k n_{2}}{8 r} + \frac{k n_{1}}{4 r} + \frac{k s_{2, 1}}{8 r} + \frac{k s_{2, 2}}{4 r}) \\ + ε_{x} ε_{y} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 1}}{4 r}) + ε_{x} ε_{z} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 2}}{4 r}) \\ + ε_{y} ε_{z} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 1}}{4 r}) + γ_{xy}^{2} (\frac{k n_{2}}{8 r} + \frac{k s_{2, 1}}{8 r} + \frac{3 k s_{1, 1}}{32 r} + \frac{k s_{2, 2}}{16 r}) \\ + γ_{xz}^{2} (\frac{k n_{2}}{8 r} + \frac{3 k s_{2, 2}}{16 r} + \frac{3 {ks}_{1, 2}}{32 r}) \\ + γ_{yz}^{2} (\frac{k n_{2}}{8 r} + \frac{k s_{2, 1}}{8 r} + \frac{k s_{1, 1}}{16 r} + \frac{k s_{2, 2}}{16 r} + \frac{{ks}_{1, 2}}{32 r}) \end{matrix}

(10)

According to the elastic mechanics used by Xu,²¹ the strain energy density can be expressed as

\begin{matrix} Π_{i} = \frac{E}{2 (1 + v)} [\frac{1 - v}{1 - 2 v} (ε_{x}^{2} + ε_{y}^{2} + ε_{z}^{2}) \\ + \frac{2 v}{1 - 2 v} (ε_{x} ε_{y} + ε_{x} ε_{z} + ε_{y} ε_{z}) + \frac{1}{2} (γ_{xy}^{2} + γ_{xz}^{2} + γ_{yz}^{2})] \end{matrix}

(11)

The strain energy densities expressed by equations (10) and (11) are equivalent, and the following assumptions of Jefferson et al.²² are used: $k s_{1, 1} = k s_{1, 2}$ and $k s_{2, 1} = k s_{2, 2}$ . Therefore, the following relationship can be obtained

\begin{matrix} \frac{3 k n_{2}}{8 r} + \frac{k n_{1}}{4 r} + \frac{3 k s_{2, 1}}{8 r} = \frac{E (1 - ν)}{2 (1 + ν) (1 - 2 ν)} \\ \frac{k n_{2}}{4 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{4 r} = \frac{E (1 - ν)}{2 (1 + ν) (1 - 2 ν)} \\ \frac{k n_{2}}{4 r} - \frac{k s_{2, 1}}{4 r} = \frac{E ν}{(1 + ν) (1 - 2 ν)} \\ \frac{k n_{2}}{8 r} + \frac{3 k s_{2, 1}}{16 r} + \frac{3 k s_{1, 1}}{32 r} = \frac{E}{4 (1 + ν)} \end{matrix}

(12)

Solving equation (12), the spring stiffness of a spherical element on boundary surface is given by

\begin{matrix} k n_{1} = \frac{4 Er (1 - ν)}{(2 ν - 1) (1 + ν)}; k s_{1, 1} = k s_{1, 2} = \frac{4 Er (3 - 2 ν)}{3 (2 ν - 1) (1 + ν)} \\ k n_{2} = \frac{- 2 Er}{(2 ν - 1) (1 + ν)}; k s_{2, 1} = k s_{2, 2} = \frac{2 Er}{(1 + ν)} \end{matrix}

(13)

Spring stiffness deduction of spherical elements on the edge

As shown in Figure 10, the spherical element on the edge of boundary is connected with nine surrounding spherical elements. There are four surrounding spherical elements with the initial central distance of $l_{0} = 2 r$ from the spherical element on the edge, and they are connected by springs on the edge. There are five surrounding spherical elements with the initial central distance of $l_{0} = 2 \sqrt{2} r$ from the spherical element on the edge, and they are connected by springs on the diagonal line. Therefore, the contact force on one spherical element on the edge includes four groups of contact forces on the edge of the cube model and five groups of contact forces on the diagonal of the cube model, as shown in Figure 11. The spring stiffness of a spherical element on the edge is deduced using the same method as was used to deduce the spring stiffness of a spherical element on the surface. The total strain energy stored per unit volume of the spherical element on the edge can be written as follows

\begin{matrix} Π_{i}^{″} = ε_{x}^{2} (\frac{k n_{2}}{4 r} + \frac{k n_{1}}{4 r} + \frac{k s_{2, 1}}{8 r} + \frac{k s_{2, 2}}{8 r}) + ε_{y}^{2} (\frac{3 k n_{2}}{16 r} + \frac{k n_{1}}{8 r} + \frac{3 k s_{2, 1}}{16 r}) + ε_{z}^{2} (\frac{3 k n_{2}}{16 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{16 r} + \frac{k s_{2, 2}}{8 r}) \\ + ε_{x} ε_{y} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 1}}{4 r}) + ε_{x} ε_{z} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 2}}{4 r}) + ε_{y} ε_{z} (\frac{k n_{2}}{4 r} - \frac{k s_{2, 1}}{4 r}) + γ_{xy}^{2} (\frac{k n_{2}}{8 r} + \frac{k s_{2, 1}}{16 r} + \frac{3 k s_{1, 1}}{32 r} + \frac{k s_{2, 2}}{32 r}) \\ + γ_{xz}^{2} (\frac{k n_{2}}{8 r} + \frac{3 k s_{2, 2}}{32 r} + \frac{3 k s_{1, 2}}{32 r}) + γ_{yz}^{2} (\frac{k n_{2}}{8 r} + \frac{k s_{2, 1}}{16 r} + \frac{k s_{1, 1}}{16 r} + \frac{k s_{2, 2}}{32 r} + \frac{k s_{1, 2}}{32 r}) \end{matrix}

(14)

The spring stiffness of a spherical element on boundary edge is given by

\begin{matrix} k n_{1} = \frac{2 Er (1 - ν)}{(2 ν - 1) (1 + ν)}; k s_{1, 1} = k s_{1, 2} = \frac{2 Er (3 + 2 ν)}{3 (2 ν - 1) (1 + ν)} \\ k n_{2} = \frac{- 2 Er}{(2 ν - 1) (1 + ν)}; k s_{2, 1} = k s_{2, 2} = \frac{2 Er}{(1 + ν)} \end{matrix}

(15)

Spring stiffness deduction of spherical elements on the corner

As shown in Figure 12, the spherical element on the corner of boundary is connected with six surrounding spherical elements. There are three surrounding spherical elements with the initial central distance of $l_{0} = 2 r$ from the spherical element on the corner, and they are connected by springs on the edge. There are three surrounding spherical elements with the initial central distance of $l_{0} = 2 \sqrt{2} r$ from the spherical element on the corner, and they are connected by springs on the diagonal line. Therefore, the contact force on one spherical element on the corner includes three groups of contact forces on the edge of cube model and three groups of contact forces on the diagonal of cube model, as shown in Figure 13. The total strain energy stored per unit volume of the spherical element on the corner can be written as follows

\begin{matrix} Π_{i}^{″} = ε_{x}^{2} (\frac{k n_{2}}{8 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{16 r} + \frac{k s_{2, 2}}{16 r}) + ε_{y}^{2} (\frac{k n_{2}}{8 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{8 r}) + ε_{z}^{2} (\frac{k n_{2}}{8 r} + \frac{k n_{1}}{8 r} + \frac{k s_{2, 1}}{16 r} + \frac{k s_{2, 2}}{16 r}) \\ + ε_{x} ε_{y} (\frac{k n_{2}}{8 r} - \frac{k s_{2, 1}}{8 r}) + ε_{x} ε_{z} (\frac{k n_{2}}{8 r} - \frac{k s_{2, 2}}{8 r}) + ε_{y} ε_{z} (\frac{k n_{2}}{8 r} - \frac{k s_{2, 1}}{8 r}) + γ_{xy}^{2} (\frac{k n_{2}}{16 r} + \frac{k s_{2, 1}}{32 r} + \frac{k s_{1, 1}}{16 r} + \frac{k s_{2, 2}}{32 r}) \\ + γ_{xz}^{2} (\frac{k n_{2}}{16 r} + \frac{k s_{2, 2}}{16 r} + \frac{k s_{1, 2}}{16 r}) + γ_{yz}^{2} (\frac{k n_{2}}{16 r} + \frac{k s_{2, 1}}{32 r} + \frac{k s_{1, 1}}{32 r} + \frac{k s_{2, 2}}{32 r} + \frac{k s_{1, 2}}{32 r}) \end{matrix}

(16)

The spring stiffness of the spherical element on the boundary corner is given by

\begin{matrix} k n_{1} = \frac{- 4 Er}{(2 ν - 1)}; k s_{1, 1} = k s_{1, 2} = \frac{- 4 Er}{(2 ν - 1) (1 + ν)} \\ k n_{2} = 0; k s_{2, 1} = k s_{2, 2} = \frac{8 Er ν}{(2 ν - 1) (1 + ν)} \end{matrix}

(17)

Numerical examples

Elastoplastic bending of thin plate

In this example, the elastoplastic bending of a thin plate subjected to uniformly distributed loads is simulated. The purpose is to study the effect of the spring stiffness of the spherical elements on the boundary on the mechanical behavior of continuity modeled using the DSEM. The FEM solution is used as the reference for evaluating the DSEM results. The geometry information, boundary conditions, and material properties are shown in Figure 14. The geometric dimensions of the thin plate are $0.3 m \times 0.3 m \times 0.006 m$ . One end of the thin plate is fixed, and the upper surface is subjected to a uniform load of $P = 4.8 \times 10^{4} Pa$ . The elastic modulus is $E = 2.1 \times 10^{11} Pa$ , Poisson’s ratio is $v = 0.24$ , and the density of the model is $ρ = 7850 kg / m^{3}$ . The yield stress is $σ_{s} = 235 MPa$ and the constitutive relation of materials is an ideal elastoplastic. The DSEM and FEM models are shown in Figures 15 and 16. In the DSEM model, the radius of the spherical elements is 1.5 mm, and the total number of spherical elements is 30,603. The thin plate is discretized into three-layer spherical elements along the z-direction. Both the upper and lower spherical elements on the surface of the boundary, and the percentage of spherical elements on the boundary in all the spherical elements, is approximately 67%. The resolution of the FEM model is $3 mm \times 3 mm \times 3 mm$ .

The simulation results of the FEM and the DSEM are shown in Figures 17 –24. The DSEM model can reproduce the same displacement and plastic distribution as the FEM model. This outcome means that the DSEM can be regarded as a valid representation of an isotropic elastoplastic material. To quantitatively study the calculation accuracy of the mechanical behavior of continuity simulated using the DSEM after revising the spring stiffness of the spherical elements on the boundary, two lines, Line I $(y = 12 mm, z = 6 mm)$ and Line II $(x = 60 mm, z = 6 mm)$ , are selected to record the x-direction, y-direction, and z-direction displacements predicted using the DSEM model. The full comparison of the displacement field predicted by FEM and DSEM along Line I $(y = 12 mm, z = 6 mm)$ on the boundary of thin plate are shown in Figures 25 –27. It can be seen that through revising the spring stiffness of the spherical elements on the boundary, the displacement errors of the boundary of the DSEM model are effectively reduced, and the accuracy of the mechanical behavior of continuity simulated using the DSEM is improved. After the spring stiffness of the spherical elements on the boundary is revised, the maximum displacement error in the x-direction decreases from 12.92% to 5.23%. The maximum displacement error in the y-direction decreases from 14.41% to 5.76%, and the maximum displacement error in the z-direction decreases from 11.38% to 4.31%. The full comparison of the displacement field along Line II $(x = 60 mm, z = 6 mm)$ on the boundary of thin plate predicted by the FEM and the DSEM are shown in Figures 28 –30. After the spring stiffness of the spherical elements on the boundary is revised, the maximum displacement error in the x-direction decreases from 12.46% to 4.56%. The maximum displacement error in the y-direction decreases from 13.11% to 4.96%, and the maximum displacement error in the z-direction decreases from 10.08% to 3.84%.

The concept of a contact force between the spherical elements in the DSEM is different from the stress concept in traditional mechanical theory. Therefore, the internal force is selected for comparison. The thin plate is only subject to the uniform load in the z-direction; thus, the internal force in the x- and y-directions is approximately $0 N$ . The results of the internal force in the z-direction predicted by the FEM and the DSEM are shown in Table 1. The internal force predicted by the DSEM is approximately identical to that of the FEM model. It can be seen that through revising the spring stiffness of the spherical elements on the boundary, the internal force errors of the DSEM model are reduced effectively. The maximum error of the internal force of the DSEM model is reduced from 3.81% to 1.64%; thus, the calculation accuracy of the mechanical behavior of continuity modeled by the DSEM is improved.

Table 1.

Results of the internal force in the z-direction.

Cross section (mm)	FEM results (N)	DEM results (N)		Error %
Cross section (mm)	FEM results (N)	Before revising thespring stiffness	After revising thespring stiffness	Before revising the spring stiffness	After revising thespring stiffness
x = 0	4320	4480.27	4386.10	3.71	1.53
x = 50	3586	3716.89	3639.07	3.65	1.48
x = 100	2894	3002.24	2937.70	3.74	1.51
x = 150	2160	2238.41	2190.89	3.63	1.43
x = 200	1426	1480.33	1449.39	3.81	1.64
x = 250	734	761.45	745.16	3.74	1.52

FEM: finite element method; DEM: discrete element method.

Dynamic buckling of thin plate

In this section, the DSEM is used to simulate the dynamic buckling of a thin plate. The purpose of this approach is to show the potential of the DSEM when solving the extremely large deformation and high material nonlinearity problems of continuity. The experimental results by So and Chen²³ are used as the reference for evaluating the DSEM results. The geometry information, boundary conditions, and material properties are specified in Figure 31. The length of the thin plate is $0.3 m$ , and the size of the cross section is $0.14 m \times 0.001 m$ . One end of the thin plate is fixed, and the other end has an attached mass of $100 kg$ .The initial velocity of the attached mass is $V_{0} = 12 m / s$ , so that this end will produce a compressive displacement along the length of the thin plate. The time step used for time integration is $4.78 \times 10^{- 5} ms$ , and the total calculation time is $14 ms$ . The elastic modulus is $E = 2.1 \times 10^{11} Pa$ , Poisson’s ratio is $v = 0.24$ , and the density of the model is $ρ = 7850 kg / m^{3}$ . The yield stress is $σ_{s} = 235 MPa$ , and the constitutive relation of materials is an ideal elastoplastic. The DSEM model is shown in Figure 32. The radius of the spherical elements is 0.5 mm. The DSEM model of the thin plate is discretized into two-layer spherical elements along the z-direction. Both the two-layer spherical elements are on the boundary of the DSEM model. Under such a condition, the thin plate can progressively form several wrinkles in the modes of a sine wave and can exhibit plastic post-buckling.

Figures 33 –35 show the front and side views of the typical deformation of the dynamic buckling of the thin plate obtained using the DSEM when the calculation times are $1 ms$ , $6 ms$ , and $14 ms$ . Under a given boundary condition, geometry information and loading system, the thin plate first undergoes wavy buckling along the length direction when $t = 1 ms$ . When $t = 6 ms$ , the first layer of obvious wrinkles is generated at the bottom of the thin plate. Then, when $t = 14 ms$ , the wrinkles in the modes of a sine wave superimpose at the bottom of the thin plate, and the second and third layers of wrinkles are gradually generated. The final deformation of the buckling mode of the DSEM model agrees well with the experimental results. Although the thin plate experiences large deformation and strong material nonlinearity, the DSEM can still effectively simulate the dynamic buckling process of the thin plate with multiple layers of wrinkles. The ability of the DSEM to solve strong nonlinear problems is verified.

Figure 36 shows the curve of the load against axial displacement for the loading end of the thin plate. When the axial load is applied, a series of peak loads of decreasing intensity are recorded during the progressive buckling of the thin plate. Through revising the spring stiffness of the spherical element on the boundary, the initial stiffness, the peak loads, and the corresponding displacement generated using the DSEM model more closely match the experimental results. After revising the spring stiffness, the first critical buckling load generated using the DSEM model decreases from $55 kN$ to $50 kN$ , and the experimental result is $48 kN$ . It is sensible to consider that the following peak loads that are lower than the first critical buckling load are caused by the severe distortion of the thin plate, which includes the formation of the wrinkles in the modes of the sine wave. Figure 37 shows the curve of the time against the axial displacement of the loading end of the thin plate. Although the reaction load fluctuates with time, the relation between the axial displacement and time is linear up to an axial displacement of $8 cm$ after revising the spring stiffness. This finding means that the speed of the attached mass is approximately constant within an axial displacement of $8 cm$ , and this behavior is identical to the experimental results.

Buckling of cracked thin plate

In this section, the DSEM is used to simulate the buckling of a cracked thin plate. The purpose of this approach is to show the ability of the DSEM when solving the fracture problems of continuity. The experimental results by Seifi and Kabiri²⁴ are used as the reference for evaluating the DSEM results. The geometry information, boundary conditions, and material properties are specified in Figure 38. The thin plate has dimensions of $0.24 m \times 0.24 m \times 0.001 m$ , with a center crack with length a; $θ$ is the orientation angle of the crack with respect to the x axis. The edge placed at the $y = - 0.12 m$ line is clamped, and the tensile load F is applied to the edge $y = 0.12 m$ . The elastic modulus is $E = 70 GPa$ , Poisson’s ratio is $v = 0.3$ , and the density of the model is $ρ = 2700 kg / m^{3}$ . The yield stress is $σ_{s} = 137 MPa$ .

Two cases are studied, as shown in Figure 39. The crack length is $a = 0.12 m$ , and the orientation angle is $θ = 0^{0}$ in the first case. In the second case, the crack length is $a = 0.168 m$ , and the orientation angle is $θ = 30^{0}$ . The DSEM model is shown in Figure 40. The radius of the spherical elements is 0.5 mm. The DSEM model of the thin plate is discretized into two-layer spherical elements along the z-direction. Both the two-layer spherical elements are on the boundary of the DSEM model. The crack is simulated by breaking the springs between the spherical elements.

The buckling forms of the thin plates are depicted in Figure 41. The buckling shapes in the tension correspond to the instability phenomenon with diffuse out-of-plane displacements localized around the cracked area. The DSEM results agree well with the experimental results of Seifi and Kabiri.²⁴ The comparison of the buckling load of the cracked thin plate is depicted in Table 2. Through revising the spring stiffness of the spherical elements on the boundary, the buckling load is reduced from $7334.90 N$ to $6775.76 N$ in the first sample, while the experimental result is $6361.03 N$ . In the second sample, the buckling load is reduced from $6775.95 N$ to $6233.57 N$ , while the experimental result is $5813.27 N$ . The errors are reduced from 15.31% to 6.52% and from 16.56% to 7.23% in the first and second samples, respectively.

Table 2.

Comparison of the buckling load.

Sample	Experimentalresult (N)	DSEM result (N)		Error %
Sample	Experimentalresult (N)	Before revising the spring stiffness	After revising the spring stiffness	Before revising the spring stiffness	After revising the spring stiffness
1	6361.03	7334.90	6775.76	15.31	6.52
2	5813.27	6775.95	6233.57	16.56	7.23

DSEM: discrete solid element method.

Conclusion

The DEM has been developed as a general and effective numerical calculation method for simulating the mechanical behavior of solid and granular materials, including the transformation from a continuum to a discontinuum. Because the mechanical behavior of the structure is calculated using the discrete grid system composed of spherical elements and springs in the DSEM, the spring stiffness has an important influence on the mechanical behavior of the structure. The spring stiffness of the spherical elements inside the DSEM model is different from that on the boundary when the continuity is modeled using the DSEM. In this article, the spring stiffness of the spherical elements on the boundary is deduced. The spring stiffness is verified using three examples. Through revising the spring stiffness of the spherical elements on the boundary, the displacement and contact force errors of the DSEM model are reduced effectively, and the calculation accuracy of continuity simulated using the DSEM is improved. In addition, through the simulation of dynamic buckling of a thin plate and the buckling of a cracked thin plate, the ability of the DSEM to address strong nonlinear, dynamic, and fracture problems is demonstrated. The following conclusions can be obtained:

According to the different positions of the spherical elements on the boundary of the DSEM model, the spherical elements on the boundary are divided into three types, which are spherical elements on the surface, on the edge, and on the corner. Based on the principle of energy conservation, the spring stiffnesses of the three types of spherical elements on the boundary are deduced.

The elastoplastic bending of a thin plate with a prominent boundary problem is simulated using the DSEM. In the DSEM model, the number of spherical elements on the boundary accounts for 67% of all spherical elements. Through revising the spring stiffness of the spherical elements on the boundary, the maximum displacement error in the x-direction decreases from 12.92% to 5.23%, the maximum displacement error in the y-direction decreases from 14.41% to 5.76%, and the maximum displacement error in the z-direction decreases from 11.38% to 4.31%. The maximum error of the internal force of the DSEM model is reduced from 3.81% to 1.64%.

The example of dynamic buckling of a thin plate subjected to compression shows that the DSEM can effectively solve the strong nonlinear and dynamic problems accompanied by multiple layers of wrinkles. After revising the spring stiffness of the spherical elements on the boundary, the error of the first critical buckling load decreases from 14.58% to 4.2% compared with the experimental results.

The example of the buckling of a cracked thin plate subjected to tension shows that the DSEM can effectively address fracture problems. Through revising the spring stiffness of the spherical elements on the boundary, the buckling load errors of a cracked thin plate decrease from 15.31% to 6.52% and from 16.56% to 7.23% in the first and second samples, respectively.

Footnotes

Handling Editor: Vesna Spasojević Brkić

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was financially supported by the Fundamental Research Funds for the Central Universities, by the Colleges and Universities in Jiangsu Province Plans to Graduate Research and Innovation (KYLX15_0089), by a project funded by the Priority Academic Program Development of the Jiangsu Higher Education Institutions, and the Natural Science Foundation of China under grant number 51538002.

ORCID iDs

Baochen Zhu

Ruoqiang Feng

References

Khoei

Anahid

Shahim

An extended arbitrary Lagrangian–Eulerian finite element method for large deformation of solid mechanics. Finite Elem Anal Des 2008; 44: 401–416.

Lynn

Isobe

Finite element code for impact collapse problems of framed structures. Int J Numer Meth Eng 2007; 69: 2538–2563.

Cundall

Strack

ODL

. A discrete numerical model for granular assemblies. Géotechnique 1979; 29: 47–65.

Ji-Hong

YE.

Application of elastic discrete element method in member structures. Eng Mech 2017; 7: 17–26.

Cheng

Liu

New discrete element models for elastoplastic problems. Acta Mech Sinica 2009; 25: 629–637.

Progressive collapse simulation based on DEM for single-layer reticulated domes. J Constr Steel Res 2017; 128: 721–731.

Hentz

Donzé

Daudeville

Discrete element modelling of concrete submitted to dynamic loading at high strain rates. Comp Struct 2004; 82: 2509–2524.

André

Iordanoff

Charles

et al . Discrete element method to simulate continuous material by using the cohesive beam model. Comp Meth Appl Mech Eng 2012; 213–216: 113–125.

Meguro

Tagel-Din

Applied element simulation of RC structures under cyclic loading. J Struct Eng 2001; 127: 1295–1305.

10.

Griffiths

Mustoe

GGV

. Modelling of elastic continua using a grillage of structural elements based on discrete element concepts. Int J Numer Meth Eng 2001; 50: 1759–1775.

11.

Sinaie

Application of the discrete element method for the simulation of size effects in concrete samples. Int J Solids Struct 2016; 108: 244–253.

12.

Hentz

Daudeville

Donzé

FV.

Identification and validation of a discrete element model for concrete. J Eng Mech 2004; 130: 709–719.

13.

Munjiza

Bangash

John

NWM

The combined finite-discrete element method. Eng Fract Mech 2004; 71: 469–483.

14.

Rousseau

Frangin

Marin

et al . Multidomain finite and discrete elements method for impact analysis of a concrete structure. Eng Struct 2009; 31: 2735–2743.

15.

Lei

Zang

An approach to combining 3D discrete and finite element methods based on penalty function method. Comput Mech 2010; 46: 609–619.

16.

Azevedo

Lemos

JV.

Hybrid discrete element/finite element method for fracture analysis. Comput Method Appl M 2006; 195: 4579–4593.

17.

Zhao

3D polycrystalline discrete element method (3PDEM) for simulation of crack initiation and propagation in granular rock. Comput Geotech 2017; 90: 96–112.

18.

Bono JPD and McDowell GR. Particle breakage criteria in discrete element modelling. Géotechnique 2016; 66: 1014–1027.

19.

McDowell

GR.

Discrete element modelling of soil particle fracture. Géotechnique 2002; 52: 131–135.

20.

Zhu

Feng

Wang

3D discrete solid element method for elastoplastic problems of continuity. J Eng Mech 2018; 144, https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29EM.1943-7889.0001459

21.

ZL.

Elastic mechanics. Beijing, China: Higher Education Press, 2006.

22.

Jefferson

Haritos

Mcmeeking

RM.

The elastic response of a cohesive aggregate —a discrete element model with coupled particle interaction. J Mech Phys Solids 2002; 50: 2539–2575.

23.

Chen

JT.

Experimental study of dynamic crushing of thin plates stiffened by stamping with V-grooves. Int J Impact Eng 2007; 34: 1396–1412.

24.

Seifi

Kabiri

AR.

Lateral load effects on buckling of cracked plates under tensile loading. Thin Wall Struct 2013; 72: 37–47.