DEFEM Method and Its Application in Pebble Flows

Abstract

Based on the concept of embedded discrete elements (EDEs), the discrete element-embedded finite element model (DEFEM) is extended in this work. The new method can be used to calculate the motion and stress variation of particles. This work discusses its application in granular flow simulation for particle motions with small deformations. The updated Lagrangian finite element method is used to obtain the coupling solution of the internal stress and the overall motion of particles in the DEFEM. The computation of deformation displacement is based on the concepts of displacement decomposition (translational and rotational motions and deformation displacement). The deformation displacement is the difference between particles and template particles [rigid body, using the discrete element method (DEM) to calculate translational and rotational displacements]. It is used to calculate the dynamic stress distribution of particles and the internal force of the node. Therefore, it has a wide scope of application (for example, it can be extended to non-spherical particles). The software validation proves the accuracy of this method. The application of the DEFEM in the accumulation process of particles is given. The motion characteristics and deformation of particles are discussed, and the stress distribution and force chain structure in particle accumulation are obtained.

1. Introduction

A particle system is a complex system composed of many discrete materials widely existing in nature and applied to industrial production [1]. In particular, it is of great engineering significance and academic value to discuss and analyze the relationship between motion characteristics and the internal stress of particles.

The discrete element method (DEM) is widely used to study the macroscopic and microscopic physical properties of granular materials. A particle is regarded as a discrete element if its motion satisfies Newton’s second law in the DEM. The DEM can quickly obtain the overall motion information of particles (position, velocity, etc.). Kačianauskas et al. [2] proposed a parallel three-dimensional DEM simulation of polydisperse materials described by normal size distribution. Fang et al. [3] developed a CUDA-GPU parallel algorithm based on super-quadric elements, which is applicable to and reliable for the large-scale engineering applications of non-spherical granular systems. Sun et al. [4] simulate the interactions of the bionic subsoilers and an ordinary subsoiler (O-S) with the soil based on the DEM. However, particles are regarded as discrete elements, and the whole particles are considered as a whole because of ignoring the internal changes of particles. One cannot obtain the internal information of particles (internal stress distribution, deformation, etc.).

When calculating the physical information inside the particles, the finite element method (FEM) is widely used [5, 6, 7]. It is a numerical solution method for elastic mechanics problems that developed rapidly with the advancement of computer power. It is applied in continuum mechanics to obtain the deformation, stress, natural frequency, and vibration mode of the structure [8]. The analyzed objects have been extended from elastic materials to plastic, viscoelastic, viscoplastic, and composite materials and from the continuum to a discontinuity [9]. The FEM is widely used in the internal stress simulation of objects. Fang et al. [3] studied the particle–wall collision process using by the FEM, which could solve local stress and strain rate. Kabir et al. [10] used the explicit FEM to simulate the flow phenomenon of particles. The results of shear behavior, particle kinetic energy, and particle stresses within the shear cell with time were given. Wagner et al. [11] proposed a new particle flow simulation method based on the extended FEM (x-FEM), which simulates moving particles without re-meshing. Krok et al. [12] conducted a systematic finite element analysis of the thermo-mechanical behavior of pharmaceutical powders during the molding process using the finite element solver ABAQUS.

The DEM and FEM have their advantages, so the combination of DEM and FEM is widely recognized and studied [13, 14, 15]. Guo and Zhao [16] proposed a multiscale framework to simulate the mechanical behavior of granular media based on DEM and FEM. A DEM assembly with the memory of its loading history is embedded in the Gauss integral points of the finite element mesh. The DEM assembly receives the global deformation at its Gauss point from the FEM as input boundary conditions in this new multiscale framework. Onate and Rojek [17] conducted a dynamic analysis of geological mechanics problems based on the combination of DEM and FEM. The combined models can employ spherical rigid and finite elements to discretize different parts of the system. Zárate and Oñate [18] proposed a new numerical method to predict the occurrence and evolution of fractures in continuous media, which combines the FEM with DEM. Munjiza and John [19] further studied the sensitivity to the mesh size of the combined single and smeared crack model in the context of the combined finite–discrete element method. Azevedo and Lemos [20] applied the hybrid method to analyze large structures. This method uses DEM to discrete the fracture zone and a discretization based on the FEM for the surrounding areas. Argilaga et al. [21] proposed a multiscale model based on an FEM × DEM approach. The method uses discrete elements in a standard finite element framework, and it has proven to be an effective way to treat real-scale engineering problems. However, it should be emphasized that the coupling of the discrete element method and finite element method still needs further research. For example, the solution of coupled motion, internal stress, and other contents need to be further discussed and analyzed.

The discrete element-embedded finite element model (DEFEM) is proposed, and it can be used to calculate particle motion and heat transfer [22, 23, 24]. This method is applied to the calculation of temperature gradient and deformation in particles. In this paper, the DEFEM is extended based on the concept of the embedded discrete element (EDE). The computation of deformation displacement is based on the concepts of displacement decomposition (translational and rotational motions and deformation displacement). The DEFEM mainly adopts the Lagrange finite element method to obtain the coupling solution about the stress and motion of particles. The finite element software verifies the relevant algorithms. The motion characteristics and deformation of particles are discussed, and the stress distribution and force chain structure in particle accumulation are obtained.

2. The discrete element-embedded finite element model

In this paper, the DEFEM is extended and used to calculate the translational and rotational motion and the stress of particles. In this method, the EDE is used to obtain the contact force distribution of the particle. The coupling solution of the internal stress and the overall motion of the particle are obtained by combining the advantages of the DEM and FEM.

2.1 Introduction of the EDE

There is a set of EDEs around the particles, which cover the outermost boundary of the particles (grid cells). In this way, the contact surfaces between the particles are transformed into contact between the EDEs. Among them, the concept of the EDE can refer to our team’s previous articles [22]. The EDE covers boundary elements of particles, and the soft-sphere model in the DEM calculates the force distribution of particles. The obtained force distribution of particles is used as the boundary condition in the FEM [23, 24]. In this paper, the two-dimensional spherical particles are divided into triangular grid elements, and the other methods for three-dimensional cases or other grid element types are similar.

2.2 Updated Lagrangian finite element method

The numerical method is used to simulate the impact and collision process, and the discretization of the object is essential. In the FEM, the object is divided into grid elements, and the motion equation satisfied is usually established by the node of each element.

The Euler method and the Lagrange method are mainly used for solving the impact and collision problem in the FEM. The Euler method fixes the computational grid in the spatial coordinates and remains unchanged in the process of deformation and motion. This method can avoid the distortion of the grid. Still, it is difficult to track the deformation boundary of the object accurately, and it requires specific processing methods to identify the shape and position of the deformation boundary. The Lagrange method is used in the DEFEM. Its grid and object remain coincident in the whole moving process; that is, the grid and object do not move relative to each other. The Lagrangian method can track and process complex deformation boundaries. Since the small deformation problem is discussed in this paper, the Lagrange method ignores the mesh distortion. In this paper, the updated Lagrangian scheme is used in the DEFEM. The Eulerian coordinates are used for derivative and integral simulation, and the stress and strain are expressed in the form of the Eulerian metric.

2.2.1 Weak form of control equation

The material derivative of the object’s momentum is equal to the sum of the external forces acting on the object. The material derivative and momentum of the object momentum are defined as follows:

$$f_i\left(t\right)=\underset{V}{\int }\rho b_i\left(x,t\right)\mathit{dV}+\underset{A}{\int }{t}_i\left(x,t\right)\mathit{dA},$$E1

$$p_i\left(t\right)=\underset{V}{\int }\rho v_i\left(x,t\right)\mathit{dV},$$E2

where $f_i$ is the material derivative of the object momentum, $b_i$ is the force (physical force) acting on the unit mass of the object, $t_i$ is the force (surface force) acting on the surface of the object, $p_i$ is the momentum of the object, $\rho$ is the density of the object, and $v_i$ is the velocity of the object.

According to the momentum equation, the Reynold’s transport theorem, and the Gaussian theorem, the following equation can be obtained:

$$\underset{V}{\int }\left(\rho \frac{{\mathit{Dv}}_i}{\mathit{Dt}}-\rho b_i-\frac{\partial {\sigma }_{\mathit{ji}}}{\partial x_j}\right)\mathit{dV}=0.$$E3

The momentum equation described by the updated Lagrangian is obtained as follows:

$$\rho \frac{{\mathit{Dv}}_i}{\mathit{Dt}}=\frac{\partial {\sigma }_{\mathit{ji}}}{\partial x_j}+\rho b_i\mathit{\text{or}}\rho \dot{\mathbf{v}}=\nabla \cdot \mathbf{\sigma }+\rho \mathbf{b}.$$E4

The boundary conditions are given as follows:

$$\begin{array}{cc}{\left.\left(n\cdot \sigma \right)\right|}_{A_t}=\overline{t},& {\left.v\right|}_{A_v}=\overline{v}\end{array},$$E5

where $A_t$ is the specified surface force boundary in the current configuration and $A_v$ is the specified velocity boundary in the current configuration.

The solution region satisfies Eq. (4). However, it is difficult to directly solve the above equation in an actual situation, so it is necessary to find an approximate solution with certain accuracy. The weighted residual method is commonly used to obtain approximate solutions to differential equations. It allows the equations and boundary conditions to have quantities at each node but requires that the weighted integrals of these quantities on the region and the boundary are zero.

The following equation can be obtained by using the weighted residual method:

$$\underset{V}{\int }\delta v_i\left(\frac{\partial {\sigma }_{\mathit{ji}}}{\partial x_j}+-\rho b_i-\rho {\ddot{u}}_i\right)\mathit{dV}=0$$E6

where $\delta v_i$ is the virtual velocity, which satisfies the following equation:

$$\begin{array}{cc}\delta v_i\left(\mathbf{X}\right)\in {\mathbb{R}}_0,& {\mathbb{R}}_0=\left\{\left.\delta {v_i}^0\left(\mathbf{X}\right)\right|\delta v_i\left(\mathbf{X}\right)\in C^0\left(\mathbf{X}\right),{\left.\delta v_i\right|}_{A_v}=0\right\}\end{array}$$E7

By using partial integration and Eq. (7), the first term of Eq. (6) can be transformed into the following form:

$${\int }_V\delta v_i\frac{\partial {\sigma }_{\mathit{ji}}}{\partial x_j}\mathit{dV}={\int }_{A_t}\delta v_i{\overline{t}}_i\mathit{dA}-{\int }_V\frac{\partial \left(\delta v_i\right)}{\partial x_j}{\sigma }_{\mathit{ji}}\mathit{dV}$$E8

Substituting Eq. (8) into Eq. (6), the weak form of the momentum equation is obtained. It is also called the principle of virtual power, as follows:

$${\int }_V\frac{\partial \left(\delta v_i\right)}{\partial x_j}{\sigma }_{\mathit{ji}}\mathit{dV}-{\int }_{A_t}\delta v_i{\overline{t}}_i\mathit{dA}-{\int }_V\delta v_i\rho b_i\mathit{dV}+{\int }_V\delta v_i\rho {\ddot{u}}_i\mathit{dV}=0$$E9

Eq. (9) can also be translated into the following form:

$$\delta \dot{p}=\delta {\dot{p}}_{\mathit{int}}-\delta {\dot{p}}_{\mathit{ext}}+\delta {\dot{p}}_{\mathit{kin}}=0$$E10

where $\delta {\dot{p}}_{\mathit{int}}$, $\delta {\dot{p}}_{\mathit{ext}}$_, and $\delta {\dot{p}}_{\mathit{kin}}$ are internal force virtual power, external force virtual power, and inertia force virtual power. Their definitions are as follows:

$$\delta {\dot{p}}_{\mathit{int}}={\int }_V\frac{\partial \left(\delta v_i\right)}{\partial x_j}{\sigma }_{\mathit{ji}}\mathit{dV},\delta {\dot{p}}_{\mathit{ext}}={\int }_{A_t}\delta v_i{\overline{t}}_i\mathit{dA}+{\int }_V\delta v_i\rho b_i\mathit{dV},\delta {\dot{p}}_{\mathit{kin}}={\int }_V\delta v_i\rho {\ddot{u}}_i\mathit{dV}$$E11

2.2.2 Finite element discretization

The entire solution region is divided into N nodes and several unit regions. In the initial configuration, the coordinate of the node is $X_1,X_2,\cdots X_{n_N}$; the node coordinate is $x_1\left(t\right),x_2\left(t\right),\cdots x_{n_N}\left(t\right)$ in the current configuration. In the FEM, the spatial coordinates of particle X at time t is $x_i\left(X,t\right)$, which can be obtained by the following equation:

$$x_i\left(X,t\right)=N_I\left(X\right)x_{\mathit{iI}}\left(t\right)$$E12

where $N_I\left(X\right)$ is the shape function of node I, and the repeated subscripts represent the summation in the range of values. Similarly, the coordinate $X_i$ of each node in the initial configuration can also be approximately expressed by the coordinate $X_{\mathit{iI}}$ of the unit node as follows:

$$X_i\left(t\right)=N_I\left(X\right)X_{\mathit{iI}}$$E13

The displacement of the node can approximately express the displacement of any point X in the element, and the specific form is given as follows:

$$u_i\left(X,t\right)=x_i\left(X,t\right)-X_i=N_I\left(X\right)u_{\mathit{iI}}\left(t\right)$$E14

Similarly, the approximate expressions of velocity and acceleration of any point X in the element are given as follows:

$${\dot{u}}_i\left(X,t\right)=N_I\left(X\right){\dot{u}}_{\mathit{iI}}\left(t\right),{\ddot{u}}_i\left(X,t\right)=N_I\left(X\right){\ddot{u}}_{\mathit{iI}}\left(t\right)$$E15

The approximate expression of virtual velocity is as follows:

$$\delta v_i\left(X\right)=N_I\left(X\right)\delta v_{\mathit{iI}}$$E16

where $\delta v_i\left(X\right)$ is virtual velocity and $\delta v_{\mathit{iI}}$ is the virtual velocity of node I.

Substituting Eqs. (16) and (14) into Eq. (9) yields the following equation:

$$\delta v_{\mathit{iI}}\left({\int }_V\frac{\partial N_I}{\partial x_j}\sigma {}_{\mathit{ji}}dV-{\int }_V{N}_I\rho b_i\mathit{dV}-{\int }_{A_t}{N}_I{\overline{t}}_i\mathit{dA}+{\int }_V{N}_I\rho N_J{\ddot{u}}_{\mathit{iJ}}\mathit{dV}\right)=0$$E17

From the boundary condition of virtual velocity, Eq. (7), and due to the arbitrariness of virtual velocity $\delta v_{\mathit{iI}}$, the following equation can be obtained:

$$\begin{array}{cc}{\int }_V\frac{\partial N_I}{\partial x_j}\sigma {}_{\mathit{ji}}dV-{\int }_V{N}_I\rho b_i\mathit{dV}-{\int }_{A_t}{N}_I{\overline{t}}_i\mathit{dA}+{\int }_V{N}_I\rho N_J{\ddot{u}}_{\mathit{iJ}}\mathit{dV}=0& \forall I\notin A_v\end{array}$$E18

Eq. (18) can also be translated into the following form:

$$\begin{array}{cc}{M}_{\mathit{IJ}}{\ddot{u}}_{\mathit{iJ}}+f_{\mathit{iI}}^{\mathrm{i}\mathit{nt}}=f_{\mathit{iI}}^{\mathit{ext}}& \forall I\notin A_v\end{array}$$E19

where $M_{\mathit{IJ}}$ is the mass matrix of the system, $f_{\mathit{iI}}^{\mathrm{i}\mathit{nt}}$ is the internal force of the node, and $f_{\mathit{iI}}^{\mathit{ext}}$ is the external force of the node. Their definitions are given as follows:

$$\left\{\begin{array}{l}{M}_{\mathit{IJ}}={\int }_V{N}_I\rho N_J\mathit{dV},f_{\mathit{iI}}^{\mathrm{i}\mathit{nt}}={\int }_V\frac{\partial N_I}{\partial x_j}\sigma {}_{\mathit{ji}}dV\\ f_{\mathit{iI}}^{\mathit{ext}}={\int }_V{N}_I\rho b_i\mathit{dV}+{\int }_{A_t}{N}_I{\overline{t}}_i\mathit{dA}\end{array}\right.$$E20

2.3 Introduction of the deformation displacement

The DEFEM can obtain the translational and rotational motion and the deformations of particles. The material constitutive model describes the response of the material of the object under an external force, which is mainly used to calculate the stress change of the object. The stress distribution of particles can be obtained using the deformation displacement and constitutive relation. During motion, particles have collisions and deformations. The displacement and rotation, small deformation, and the elastic material particle stress are discussed in this paper. The analysis of other constitutive models is similar to that of elastic materials.

In the DEFEM, it is necessary to decompose the overall displacement of particles when obtaining the deformation displacement, as shown in the following formula:

$$u=u_{\mathit{\text{move}}}+u_{\mathit{\text{rotate}}}+u_{\mathit{\text{deform}}}$$E21

where $u$ is the absolute displacement of the node, $u_{\mathit{\text{move}}}$ is the translational motion displacement, $u_{\mathit{\text{rotate}}}$ is the rotational displacement, and $u_{\mathit{\text{deform}}}$ is the deformation displacement. The displacement of any point can be decomposed into the above three parts, in which the deformation displacement is used to calculate the strain and stress of the element.

The translational motion displacement $u_{\mathit{\text{move}}}$ and rotational displacement $u_{\mathit{\text{rotate}}}$ can be solved by using the DEM. In the DEFEM, it is considered that there is a virtual template particle, which is a rigid object, and it only moves and rotates without deformation. The translational motion and rotation of the template particles satisfy Newton’s equations. The translational displacement and rotation angle of the template particles at each time interval can be obtained by using the DEM and the central difference method, as shown in the following equation:

$${\chi }^{t+\mathrm{\Delta }t}={\chi }^{t-\mathrm{\Delta }t}+{\dot{\chi }}^{t+\mathrm{\Delta }t/2}\mathrm{\Delta }t,{\theta }^{t+\mathrm{\Delta }t}={\theta }^{t-\mathrm{\Delta }t}+{\dot{\theta }}^{t+\mathrm{\Delta }t/2}\mathrm{\Delta }t$$E22

where ${\chi }^{t+\mathrm{\Delta }t}$ and ${\chi }^{t-\mathrm{\Delta }t}$ are the translational displacements of particle i at time $t+\mathrm{\Delta }t$ and $t-\mathrm{\Delta }t$, ${\theta }^{t+\mathrm{\Delta }t}$ and ${\theta }^{t-\mathrm{\Delta }t}$ are the rotation angles of particle i at time $t+\mathrm{\Delta }t$ and $t-\mathrm{\Delta }t$, ${\dot{\chi }}^{t+\mathrm{\Delta }t/2}$ is the translational velocity of particle i at time $t+\mathrm{\Delta }t/2$, and ${\dot{\theta }}^{t+\mathrm{\Delta }t/2}$ is the angular velocity of particle i at time $t+\mathrm{\Delta }t/2$. The formulas of translational velocity and angular velocity are as follows:

$${\dot{\chi }}^{t+\mathrm{\Delta }t/2}={\dot{\chi }}^{t-\mathrm{\Delta }t/2}+\frac{\mathrm{\Delta }t}{m_i}\sum _{j=1}^n{F}_i^t,{\dot{\theta }}^{t+\mathrm{\Delta }t/2}={\dot{\theta }}^{t-\mathrm{\Delta }t/2}+\frac{\mathrm{\Delta }t}{J_i}\sum _{j=1}^n{T}_i^t$$E23

where $m_i$ is the mass of particle i, $J_i$ is the rotational inertia of particle i, $n$ is the number of particles in contact with particle i, $F_i^t$ is the contact force between particle i and other contact particles at time $t$, $T_i^t$ is the torque between particle i and other contact particles at time $t$, ${\sum }_{j=1}^n{F}_i^t$ and ${\sum }_{j=1}^n{T}_i^t$ are the resultant force and resultant torque of particle i at time $t$, ${\dot{\chi }}^{t+\mathrm{\Delta }t/2}$ and ${\dot{\chi }}^{t-\mathrm{\Delta }t/2}$ are the translational velocities of particle i at time $t+\mathrm{\Delta }t/2$ and $t-\mathrm{\Delta }t/2$, and ${\dot{\theta }}^{t+\mathrm{\Delta }t/2}$ and ${\dot{\theta }}^{t-\mathrm{\Delta }t/2}$ are the angular velocities of particle i at time $t+\mathrm{\Delta }t/2$ and $t-\mathrm{\Delta }t/2$.

The translational velocity and angular velocity of the template particles can be calculated by Eq. (23). The translational displacement and rotation displacement of the template particles can be obtained by Eq. (22), which is also considered to be the translational displacement $u_{\mathit{\text{move}}}$ and rotational displacement $u_{\mathit{\text{rotate}}}$ of the particle after the decomposition of the overall displacement. Then, the deformation displacement $u_{\mathit{\text{deform}}}$ of the particle can be obtained by using Eq. (21).

2.4 Numerical process of the discrete element-embedded finite element model

In this paper, the DEFEM is extended to calculate the translational, rotational, and stress variations of particles. In this method, the EDE is used to allocate the resultant contact force and torque to the boundary node of particles.

The concept of template particles is used in the whole simulation process, which only moves and rotates without deformation. The deformation displacement of the particles is obtained by calculating the difference between the coordinates of the particles and the template particles simultaneously. The related process is mainly based on displacement decomposition, as shown in Eq. (21). The deformation displacement here mainly has two functions in the DEFEM: [8] The stress distribution of particles is calculated according to the constitutive relation (2) The internal force of the node is calculated as a part of the finite element equation for the next time, as shown in Eq. (19). It participates in the whole simulation cycle, as shown in Figure 1.

Figure 1.

Numerical process diagram of deformation displacement.

The numerical process of the DEFEM is as follows: [8] Initialization and conditionalization: The initial time and initial conditions are assigned, and the state of the template particle is the same as that of the particle (node coordinates, node velocity, acceleration, etc.). [21] The resultant force and moment (contact force with other particles or the wall) of particles are calculated according to the soft-sphere model in the DEM. [20] The DEM is used to calculate the next time information of the template particle (node coordinates, node velocity, acceleration, etc.). [3] The resultant force of particles is transformed into the force of the EDE, and finally, the force of the particle boundary node is obtained. [5] The information about the next time step of the particle is calculated by the updated Lagrangian FEM. The coordinate difference between the particle and the template particle is denoted as the deformation displacement, and the stress distribution of the particle at this moment and the internal force of the next time step are calculated. [16] Return to step [21] until all particles are calculated. [6] Update the simulation time. [1] Output: If the program calculates the simulation deadline or maximum steps, stop the program. Otherwise, return to step [21] to calculate the particle information at the next moment (Figure 2).

Figure 2.

Schematic diagram of the simulation process of the DEFEM.

2.5 Software validation of particle–wall collision

In this paper, the ABAQUS software is used to simulate the particle–wall collision process, and it is compared with the numerical results by using the DEFEM. The ABAQUS is a finite element software applied to engineering simulation, which can simulate the stress and strain of large structures. This paper simulates the collision process between a spherical particle and the wall using the ABAQUS. The particle with a diameter of 0.05 m vertically hit the wall at a velocity of 0.5 m/s. The elastic modulus of the particles is 0.02 GPa, and the Poisson’s ratio is 0.3. The elastic modulus of the wall is 2 GPa, and the Poisson’s ratio is 0.3. After the contact force of particles is counted, the simulation results of ABAQUS and DEFEM are compared, as shown in Figure 3(a). It can be observed that the changing trend of the two results is the same, and the numerical value is similar. The contact force of the two numerical methods is re-zero near 0.003 seconds. In addition, the displacement of particles during collision and rebound is tracked, as shown in Figure 3(b). The results obtained by the two methods are in good agreement. It can be seen from these two results that the simulation results of the DEFEM are consistent with those of the ABAQUS.

Figure 3.

(a) Comparison of the results of DEFEM and ABAQUS software for contact force of the particle; (b) comparison of DEFEM and ABAQUS software for the displacement of the particle.

Compared with the DEM, one of the advantages of the DEFEM is to calculate the dynamic deformation displacement and stress distribution of particles. It shows the comparison of deformation displacement calculated using the ABAQUS and the DEFEM, as shown in Figure 4. It can be observed that the distribution of the two methods is similar. The black vector in Figure 4(b) represents the deformation displacement at the grid node of the particle. Due to the collision between particles and the wall, the deformation near the wall is greater than that away from the wall, which also conforms to the relevant physical law.

Figure 4.

(a) Static deformation displacement distribution calculated by ABAQUS software; (b) dynamic deformation displacement distribution calculated by the DEFEM (black vector represents the deformation displacement of nodes).

As shown in Figure 5, the static stress distribution of the particle calculated using the ABAQUS software has the same variation trend as the dynamic stress distribution calculated by the DEFEM, but there are some differences in values. Because the FEM is used to solve the problem in the ABAQUS software, it is considered that the contact force between particles only exists at the contact point. As shown in Figure 5(b), face contact is selected in the DEFEM. The multiple EDEs cover the boundary of particles. Therefore, the contact between particles can be transformed into the contact between EDEs, involving multiple nodes. The judgment and simulation process of the contact force are different between the two methods, so the stress distribution calculated by the two methods is different, but these differences are acceptable within the allowable error range.

Figure 5.

(a) Static stress distribution calculated by the ABAQUS software; (b) dynamic stress distribution calculated by the DEFEM (black vector represents the contact force of the node).

Through the comparison between the ABAQUS software and the DEFEM, it can be observed that the DEFEM can well track the contact force, stress distribution, and movement in the process of particle collision, which also proves the accuracy of the DEFEM.

3. Simulation of particle accumulation by DEFEM

In this paper, the particle accumulation process is simulated by using the DEFEM. The motion characteristics of particles in the whole process of accumulation in pebble beds are discussed and analyzed. The stress distribution and force chain structure during particle accumulation are also obtained.

3.1 Geometric parameter settings

The particle accumulation process in a two-dimensional pebble bed is simulated. As shown in Figure 6, the bottom inclination angle of the pebble bed is changed in a certain range. The bottom diameter of the pebble bed is 0.2 m, and the bottom inclination angle is 15 degrees or 45 degrees, as shown in Figure 6. In the initial state, the particles with a regular arrangement move downward vertically at a certain velocity. The particles collide with the wall of the pebble bed or other particles and finally form a stable accumulation structure at the bottom of the pebble bed. Two regions in the pebble bed are defined: the central and the near-wall regions. The particle motion characteristics in two different regions are discussed. In the case of this paper, the near-wall region is within 0.2 m from the left and right walls, and the central region is within 0.2 m from the central axis of the pebble bed (the selection of values with careful consideration of particle diameter and pebble bed size). The velocity, stress, and deformation during particle accumulation are discussed and analyzed under the different pebble bed geometry structures. The specific parameters in the numerical simulation are shown in Table 1.

Figure 6.

Diagram of pebble bed structure.

Parameters	Values	Units
Elastic modulus	2 × 10^7	Pa
Poisson’s ratio	0.3
Density	1000	kg/m^3
Diameter	0.05	m
Initial velocity	0.5	m/s
Number of particles	72
Time step	5 × 10^−6	s
Friction coefficient	0.18

Table 1.

Parameters of the particle accumulation process.

3.2 Motion characteristics in the near-wall and central regions

The particles initially fall vertically at a certain velocity, forming a stable accumulation structure in the pebble bed. The particle motion characteristics are different in the accumulation process of the different shapes of pebble beds. Due to the effect of gravity, particles may collide with the wall or other particles. If the particles collide, the velocity of the particles changes due to the effect of contact force. The average velocity and average velocity difference of particles in different regions are defined as follows:

$$v_{\mathit{ave}}=\frac{1}{m}\sum _{i=0}^m{v}_i$$E24

$$v_{\mathit{\text{diff}}}=\frac{1}{m}\sum _{i=0}^m\left(v_i-\frac{1}{n}\sum _{j=0}^n{v}_j\right)$$E25

where $v_{\mathit{ave}}$ is the average velocity of particles in the region, $v_i$ is the velocity of particle $i$ in the region, $n$ is the number of particles in the region. $v_{\mathit{\text{diff}}}$ is the average velocity difference of particles in the region, $v_j$ is the velocity of particle $j$ contacting particle $i$, and $m$ is the total number of other particles in contact with particle $i$.

The difference in the average particle velocity between the central and near-wall regions at the bottom inclination of 45 degrees is shown in Figure 7. It can be seen that the horizontal velocity of particles is zero in the initial state. Moreover, the vertical velocity of particles increases steadily at the beginning of the accumulation process due to the effect of gravity. Then, it can be observed that some particles collide with the side wall, so the average velocity and average velocity difference of particles near the wall region fluctuate first.

Figure 7.

(a) Average velocity of particles in the central region and near-wall region with the bottom inclination of 45 degrees; (b) average velocity difference of particles in the central region and near-wall region with the bottom inclination of 45 degrees.

In terms of vertical velocity, the peak velocity of particles in the central region is much larger than that in the near-wall region. This is because compared with the particles in the near-wall region, the particles in the central region can keep free falling motion for a longer time before colliding with the wall or other particles. The change in the vertical velocity of particles in the central region after the collision is more obvious; that is, the vertical velocity difference of central velocity particles is significantly greater than that of the near-wall region. In terms of the horizontal velocity, it is more about the impact of particles. Compared with the near-wall region, the probability of particle collision in the central region is greater. Therefore, the fluctuation of the average horizontal velocity and velocity difference in the central region is significantly greater than that in the near-wall region. Finally, the particles reach the state of stable accumulation; that is, the particles remain stationary (the average velocity and velocity difference of the particles are zero).

3.3 Deformation of particles in the near-wall and central regions

A major advantage of the DEFEM compared with the DEM is that it can calculate the deformation of particles. This part discusses the whole deformation of particles in the accumulation process of the pebble bed. The accumulation process of two-dimensional particles is simulated in this paper, so the deformation of particles is represented by the change of area. The average coefficient of area ratio in the defined region is given as follows:

$$S_{\mathit{ave}}=\frac{1}{m}\sum _{i=0}^m\frac{S_i}{S_0}$$E26

where $S_{\mathit{ave}}$ is the average coefficient of area ratio, $m$ is the number of particles in the region, $S_i$ is the area of particles at the current time, and $S_0$ is the area of the initial moment of the particle (no deformation state). If the particle does not deform, the average coefficient of area ratio is equal to 1 ($S_{\mathit{ave}}=1$); if the particle is compressed, the average coefficient of area ratio is less than 1 ($S_{\mathit{ave}}<1$).

The average coefficient of area ratio in the central and near-wall regions with the bottom inclinations of 45 degrees and 15 degrees is shown in Figure 8. It can be seen that some particles are compressed (the average coefficient of area ratio is less than 1; $S_{\mathit{ave}}<1$) during the accumulation process. When the bottom inclination is 45 degrees, the absolute value of the peak about the average coefficient of area ratio is less than that of 15 degrees, and the changing trend is the same as that of the average stress of the particles. When the particle has severe deformation, the average coefficient of the area ratio and absolute internal value of stress also increases.

Figure 8.

(a) The average coefficient of area ratio in the central region and near-wall region with the bottom inclination of 45 degrees; (b) the average coefficient of the area ratio in the central region and near-wall region with the bottom inclination of 15 degrees.

3.4 Stress of particles in the near-wall and central regions

Particles may have strain and stress in the accumulation process because they may collide with other particles or walls. In this paper, the von Mises stress is used, and its formula is as follows:

$$\sigma =\sqrt{\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]}$$E27

Where $\sigma$ is the von Mises stress and ${\sigma }_1$, ${\sigma }_2$, and ${\sigma }_3$ are the first, second, and third principal stresses, respectively.

The average stress of particles in the region is defined as:

$${\sigma }_{\mathit{ave}}=\frac{1}{m}\sum _{i=0}^m{\sigma }_i$$E28

where ${\sigma }_{\mathit{ave}}$ is the average stress in the region, ${\sigma }_i$ is the stress of particle $i$ in the region, and $m$ is the number of particles in the region.

The average stress of particles in the central and near-wall regions with the bottom inclinations of 45 degrees and 15 degrees is shown in Figure 9.

Figure 9.

(a) Average stress of particles with the bottom inclination of 45 degrees; (b) average stress of particles with the bottom inclination of 15 degrees.

[8] Different effects of bottom inclination.

It can be seen from Figure 9 that the peak stress is smaller when the bottom inclination is 45 degrees. The initial positions of particles at different bottom inclination angles are the same. When the bottom inclination is 45 degrees, the particles have contact with the wall on both sides earlier in the process of falling, which can also be verified as shown in Figure 9(b). When the bottom inclination is 15 degrees, the stress remains zero at the beginning, and then the stress begins to become non-zero around 0.185 seconds gradually. If the bottom inclination angle is smaller, the acceleration time of particles is longer due to the action of gravity, and the velocity of collision with the wall is larger. Hence, the contact force and stress are larger. In the stable state formed during the accumulation process, it can be seen that there is a great difference in the average stress of particles between the near-wall region and central region when the bottom inclination is 45 degrees. When the bottom inclination is 15 degrees, there is little difference between the two. When the bottom inclination is larger, the pressure on the side wall of the pebble bed is larger, and the average stress difference between the particles in the near-wall region and central region is larger.

[21] Differences between the near-wall region and central region.

Due to the structure of the pebble bed, particles come into contact with the side wall first, so the stress value of particles in the near-wall region changes first compared with that in the central region. However, the collision of particles in the central region is more intense, so the peak stress in the central region is greater than that in the near-wall region. Moreover, the velocity of particles changes after being affected by the wall, and the velocity-changed particles are more likely to collide with the particles in the central region. Therefore, there is a certain phase difference between the average stress change of particles in the near-wall region and the central region; that is, the average stress change of particles in the central region always lags behind that in the near-wall region. The average stress of particles in the central region is greater than that in the near-wall region.

3.5 Dynamic stress distribution and force chain of particles

In the final stable stage, the change of particle position is not drastic. Particles have a stable accumulation structure and force chain due to gravity. Under different bottom inclinations, the internal stress distribution and force chain are given, as shown in Figure 10. It can be seen that the accumulated particles form a stable force chain structure, and the stress distribution and force chain distribution have a certain corresponding relationship. The stress of particles far from the force chain is smaller, and the stress near the force chain is larger in the stable force chain structure. The red vector in Figure 10 represents the contact force of the node.

Figure 10.

(a) Dynamic stress distribution and force chain of particles with the bottom inclination of 45 degrees; (b) dynamic stress distribution and force chain of particles with the bottom inclination of 15 degrees.

4. Conclusions

Based on the concept of the EDE, the DEFEM is extended for simulating translational and rotational motions and small deformations of particles.

When the particle contacts other particles, it is assumed that there is a group of EDEs on the outermost boundary element of the particle. The updated Lagrangian finite element method is used to obtain the coupling solution of internal stress and the overall motion of particles in the DEFEM. It is necessary to decompose the overall displacement of particles into three parts: the translational and rotation displacement and the deformation displacement. The concept of template particles is proposed in the whole simulation process, which only moves and rotates without deformation. The translational displacement and rotation displacement of the template particles can be obtained by using the DEM. The difference between the updated coordinates of the particles and the template particles is considered the deformation displacement in the DEFEM. And it is used to calculate the stress distribution of particles and the internal force of the node. Therefore, this method has a wide range of applications, such as the simulation of non-spherical particles.

The accuracy of this method is proved by the software validation. The application of the DEFEM in the particle accumulation process is given. The motion characteristics and deformation of particles are discussed, and the stress distribution and force chain structure in particle accumulation is obtained in this paper. This paper extends the application of the DEFEM in translational and rotational motions of particles with deformations. This new method can also be used to solve multi-physical field coupling problems such as thermal and mechanical coupling in the future.

Acknowledgments

The authors are grateful for the support of this research by the National Science and Technology Major Project (2011ZX06901-003). One of the authors (T. Hibiki) appreciates the support through the Global STEM Professorship from Hong Kong government.