**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(t) = \int\_V \rho b\_i(\mathbf{x}, t)dV + \int\_A t\_i(\mathbf{x}, t)dA,\tag{1}$$

$$p\_i(t) = \int\_V \rho v\_i(\mathbf{x}, t)dV,\tag{2}$$

where *fi* is the material derivative of the object momentum, *bi* is the force (physical force) acting on the unit mass of the object, *ti* is the force (surface force) acting on the surface of the object, *pi* is the momentum of the object, *ρ* is the density of the object, and *vi* 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:

$$\int\_{V} \left( \rho \frac{Dv\_i}{Dt} - \rho b\_i - \frac{\partial \sigma\_{ji}}{\partial \mathbf{x}\_j} \right) dV = \mathbf{0}.\tag{3}$$

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

$$
\rho \frac{Dv\_i}{Dt} = \frac{\partial \sigma\_{ji}}{\partial \mathbf{x}\_j} + \rho b\_i \text{ or } \rho \dot{\mathbf{v}} = \nabla \cdot \boldsymbol{\sigma} + \rho \mathbf{b}. \tag{4}
$$

The boundary conditions are given as follows:

$$(\boldsymbol{n}\cdot\boldsymbol{\sigma})|\_{A\_t} = \overline{t}, \quad \boldsymbol{v}|\_{A\_v} = \overline{\boldsymbol{\sigma}},\tag{5}$$

where *At* is the specified surface force boundary in the current configuration and *Av* 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:

$$\int\_{V} \delta v\_{i} \left( \frac{\partial \sigma\_{ji}}{\partial \mathbf{x}\_{j}} + -\rho b\_{i} - \rho \ddot{u}\_{i} \right) dV = \mathbf{0} \tag{6}$$

where *δvi* is the virtual velocity, which satisfies the following equation:

$$\delta v\_i(\mathbf{X}) \in \mathbb{R}\_0, \quad \mathbb{R}\_0 = \left\{ \delta v\_i^0(\mathbf{X}) |\delta v\_i(\mathbf{X}) \in \mathcal{C}^0(\mathbf{X}), \delta v\_i|\_{A\_v} = 0 \right\} \tag{7}$$

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\_{ji}}{\partial \mathbf{x}\_{j}} dV = \int\_{A\_{t}} \delta v\_{i} \tilde{t}\_{i} dA - \int\_{V} \frac{\partial (\delta v\_{i})}{\partial \mathbf{x}\_{j}} \sigma\_{ji} dV \tag{8}
$$

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(\delta v\_{i})}{\partial \mathbf{x}\_{j}} \sigma\_{jt} dV - \int\_{A\_{t}} \delta v\_{i} \tilde{t}\_{i} dA - \int\_{V} \delta v\_{i} \rho b\_{i} dV + \int\_{V} \delta v\_{i} \rho \ddot{u}\_{i} dV = \mathbf{0} \tag{9}
$$

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

$$
\delta \dot{p} = \delta \dot{p}\_{\text{int}} - \delta \dot{p}\_{\text{ext}} + \delta \dot{p}\_{\text{kin}} = \mathbf{0} \tag{10}
$$

where *δp*\_ *int*, *δp*\_*ext*, and *δp*\_ *kin* are internal force virtual power, external force virtual power, and inertia force virtual power. Their definitions are as follows:

$$\delta \dot{p}\_{\rm int} = \int\_{V} \frac{\partial (\delta v\_{i})}{\partial \mathbf{x}\_{j}} \sigma\_{\dot{\rm l}} dV, \qquad \delta \dot{p}\_{\rm ext} = \int\_{A\_{\rm t}} \delta v\_{i} \overline{\mathbf{f}}\_{i} dA + \int\_{V} \delta v\_{i} \rho b\_{i} dV, \qquad \delta \dot{p}\_{\rm kin} = \int\_{V} \delta v\_{i} \rho \ddot{u}\_{i} dV \tag{11}$$

#### *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, ⋯*XnN* ; the node coordinate is *x*1ð Þ*t* , *x*2ð Þ*t* , ⋯*xnN* ð Þ*t* in the current configuration. In the FEM, the spatial coordinates of particle *X* at time *t* is *xi*ð Þ *X*, *t* , which can be obtained by the following equation:

$$\mathbf{x}\_{i}(\mathbf{X},t) = \mathbf{N}\_{I}(\mathbf{X})\mathbf{x}\_{iI}(t) \tag{12}$$

where *NI*ð Þ *X* is the shape function of node I, and the repeated subscripts represent the summation in the range of values. Similarly, the coordinate *Xi* of each node in the initial configuration can also be approximately expressed by the coordinate *XiI* of the unit node as follows:

$$X\_i(t) = N\_I(X)X\_{iI} \tag{13}$$

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(X,t) = x\_i(X,t) - X\_i = N\_l(X)u\_{il}(t) \tag{14}$$

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

$$
\dot{u}\_i(X,t) = N\_I(X)\dot{u}\_{il}(t), \quad \ddot{u}\_i(X,t) = N\_I(X)\ddot{u}\_{il}(t) \tag{15}
$$

The approximate expression of virtual velocity is as follows:

$$
\delta v\_i(X) = N\_I(X)\delta v\_{iI} \tag{16}
$$

where *δvi*ð Þ *X* is virtual velocity and *δviI* is the virtual velocity of node I. Substituting Eqs. (16) and (14) into Eq. (9) yields the following equation:

$$
\delta v\_{il} \left( \int\_{V} \frac{\partial N\_{l}}{\partial \mathbf{x}\_{j}} \sigma\_{ji} dV - \int\_{V} N\_{l} \rho b\_{i} dV - \int\_{A\_{l}} N\_{l} \overline{t}\_{i} dA + \int\_{V} N\_{l} \rho N\_{l} \ddot{u}\_{l} dV \right) = \mathbf{0} \tag{17}
$$

From the boundary condition of virtual velocity, Eq. (7), and due to the arbitrariness of virtual velocity *δviI*, the following equation can be obtained:

$$\int\_{V} \frac{\partial N\_{l}}{\partial \mathbf{x}\_{j}} \sigma\_{jl} dV - \int\_{V} N\_{l} \rho b\_{l} dV - \int\_{A\_{t}} N\_{l} \overline{t}\_{i} dA + \int\_{V} N\_{l} \rho N\_{l} \ddot{u}\_{lj} dV = 0 \quad \forall I \notin A\_{v} \tag{18}$$

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

$$\mathcal{M}\_{\text{II}}\ddot{u}\_{\text{i}\text{I}} + f\_{\text{i}\text{I}}^{\dot{\text{out}}} = f\_{\text{i}\text{I}}^{\text{ext}} \quad \forall \text{I} \notin \mathcal{A}\_{v} \tag{19}$$

where *MIJ* is the mass matrix of the system, *f* i*nt iI* is the internal force of the node, and *f ext iI* is the external force of the node. Their definitions are given as follows:

$$\begin{cases} M\_{\text{II}} = \int\_{V} N\_{l} \rho N\_{l} dV, \quad f\_{\text{il}}^{\text{int}} = \int\_{V} \frac{\partial N\_{l}}{\partial \mathbf{x}\_{j}} \sigma\_{j\text{t}} dV \\\\ f\_{\text{il}}^{\text{ext}} = \int\_{V} N\_{l} \rho b\_{i} dV + \int\_{A\_{t}} N\_{l} \mathbf{f}\_{i} dA \end{cases} \tag{20}$$

#### **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:

$$
\mu = \mathfrak{u}\_{move} + \mathfrak{u}\_{rotate} + \mathfrak{u}\_{defom} \tag{21}
$$

where *u* is the absolute displacement of the node, *umove* is the translational motion displacement, *urotate* is the rotational displacement, and *udeform* 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 *umove* and rotational displacement *urotate* 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+\Delta t} = \chi^{t-\Delta t} + \dot{\chi}^{t+\Delta t/2} \Delta t, \qquad \theta^{t+\Delta t} = \theta^{t-\Delta t} + \dot{\theta}^{t+\Delta t/2} \Delta t \tag{22}$$

where *<sup>χ</sup>t*þΔ*<sup>t</sup>* and *<sup>χ</sup>t*�Δ*<sup>t</sup>* are the translational displacements of particle *<sup>i</sup>* at time *<sup>t</sup>* <sup>þ</sup> <sup>Δ</sup>*<sup>t</sup>* and *<sup>t</sup>* � <sup>Δ</sup>*t*, *<sup>θ</sup>t*þΔ*<sup>t</sup>* and *<sup>θ</sup>t*�Δ*<sup>t</sup>* are the rotation angles of particle *<sup>i</sup>* at time *<sup>t</sup>* <sup>þ</sup> <sup>Δ</sup>*<sup>t</sup>* and *<sup>t</sup>* � <sup>Δ</sup>*t*, *χ*\_ *<sup>t</sup>*þΔ*t=*<sup>2</sup> is the translational velocity of particle *<sup>i</sup>* at time *<sup>t</sup>* <sup>þ</sup> <sup>Δ</sup>*t=*2, and \_ *θ <sup>t</sup>*þΔ*t=*<sup>2</sup> is the angular velocity of particle *i* at time *t* þ Δ*t=*2. The formulas of translational velocity and angular velocity are as follows:

$$\dot{\chi}^{t+\Delta t/2} = \dot{\chi}^{t-\Delta t/2} + \frac{\Delta t}{m\_i} \sum\_{j=1}^{n} F\_i^t, \qquad \dot{\theta}^{t+\Delta t/2} = \dot{\theta}^{t-\Delta t/2} + \frac{\Delta t}{J\_i} \sum\_{j=1}^{n} T\_i^t \tag{23}$$

where *mi* is the mass of particle *i*, *Ji* is the rotational inertia of particle *i*, *n* is the number of particles in contact with particle *i*, *F<sup>t</sup> <sup>i</sup>* is the contact force between particle *i* and other contact particles at time *t*, *T<sup>t</sup> <sup>i</sup>* is the torque between particle *i* and other contact particles at time *t*, P*<sup>n</sup> <sup>j</sup>*¼<sup>1</sup>*F<sup>t</sup> <sup>i</sup>* and <sup>P</sup>*<sup>n</sup> <sup>j</sup>*¼<sup>1</sup>*T<sup>t</sup> <sup>i</sup>* are the resultant force and resultant torque of particle *i* at time *t*, *χ*\_ *<sup>t</sup>*þΔ*t=*<sup>2</sup> and *χ*\_ *<sup>t</sup>*�Δ*t=*<sup>2</sup> are the translational velocities of particle *<sup>i</sup>* at time *<sup>t</sup>* <sup>þ</sup> <sup>Δ</sup>*t=*2 and *<sup>t</sup>* � <sup>Δ</sup>*t=*2, and \_ *θ <sup>t</sup>*þΔ*t=*<sup>2</sup> and \_ *θ <sup>t</sup>*�Δ*t=*<sup>2</sup> are the angular velocities of particle *i* at time *t* þ Δ*t=*2 and *t* � Δ*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 *umove* and rotational displacement *urotate* of the particle after the decomposition of the overall displacement. Then, the deformation displacement *udeform* 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**.

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,

#### **Figure 1.**

*Numerical process diagram of deformation displacement.*

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**).
