Nonlinear wave in granular systems based on elastoplastic dashpot model

The dynamic dashpot models are widely used in EDEM commercial software. However, most dashpot models suffer from a serious numerical issue in calculating the granular chain because the denominator of damping force includes the initial impact velocity. Moreover, the existing dynamic dashpot models extended from the original Hertz contact law overestimated the contact stiffness in the elastoplastic contact phase. These two reasons above result in most dynamic dashpot models confronting some issues in calculating the multiple collision of the granular chain. Therefore, this investigation aims to propose a new composite dynamic dashpot model for the dynamic simulation of granular matters. First, the entire contact process is divided into three different phases: elastic, elastoplastic, and full plastic phases. The Hertz contact stiffness is still used in the elastic contact phase when the contact comes into the elastoplastic or full plastic phase. Hertz contact stiffness in the dynamic dashpot model is replaced by linearizing the contact stiffness from the Ma‐Liu (ML) model in each time step. Second, the whole contact behavior is treated as a linear mass‐spring‐damper model, and the damping factor is obtained by solving the single‐degree‐freedom underdamped vibration equation. The new dynamic dashpot model is proposed by combining the contact stiffnesses in different contact phases and corresponding damping factors, which not only removes the initial impact velocity from the denominator of damping force but also updates the contact stiffness based on the constitutive relation of the contact body when the contact comes into the elastoplastic or full plastic phase. Finally, a granular chain is treated as numerical examples to check the reasonability and effectiveness of the new dynamic dashpot model by comparing it to the experimental data. The simulation shows that the solitary waves obtained using the new dashpot model are more accurate than the dashpot model used in EDEM software.


| INTRODUCTION
The dynamics of one-dimensional granular chains 1 exhibit a rich landscape of unforeseen behavior phenomena by simple physical laws, which supports the formation of highly nonlinear, strongly localized solitary waves. 2 This feature makes the granular chains a perfect candidate for various potential applications, [3][4][5] such as shock absorbers, impulse protectors, impact mitigation, acoustic switch, waveguide, and vibration filter, and so on. Although the granular chain has many known applications, [6][7][8][9] it is hard to understand its intrinsic dynamic properties caused by the strong nonlinearity of contact force between the particles. [10][11][12][13][14][15] The nonlinearity of the granular interaction is closely related to the local properties 16,17 of the contact body and energy dissipation 18,19 during impact.
In earlier works, the mechanics characteristic between the particles is described by the Hertz contact law. 20,21 To consider the energy dissipation during contact, 22 the equation of motion of granular chain is improved by introducing the dissipative term 23 related to the relative contact velocity between the particles or the coefficient of restitution 18 to represent the energy dissipation. 19 In this process, the energy dissipation during contact is not represented by the damping factor 24 from the continuous dashpot models.
Namely, the contact force between the particles is still calculated based on the original Hertz contact law. 25 This is one of the reasons why the existing dashpot models are seldom applied to the simulation of the granular chain. 2 More importantly, almost all of the existing dashpot models 26 derived based on the energy conservation [27][28][29][30][31] before and after impact suffered from the serious numerical singular problem 32 because the denominator of damping force [33][34][35] in the dashpot models includes the initial impact velocity. 36 In this scenario, when the initial impact velocity of the particle in the granular chain is equal to zero, the damping force in the dashpot model is infinite, 32 which violates the mechanics characteristic between the particles, and triggers the numerical issues. To sidestep the numerical singular issue, 31,37,38 some damping factors are assumed as constant according to the empirical evidence. 39,40 However, this approach is not a universal method regarding the different impact scenarios. The alternative method 41,42 is to develop the damping factor based on the spring-damper model, which removes the initial impact velocity from the denominator of the damping factor in the dashpot model. Therefore, the dashpot model 41 without the initial impact velocity in the denominator of damping force is used in the commercial EDEM software. Nevertheless, most existing dashpot models 26 are difficult to accurately describe the contact behavior involving elastoplastic or plastic deformation. 43 That is mainly because the Hertz contact stiffness 44,45 overestimates the contact stiffness in the elastoplastic or plastic contact phase. Accordingly, to overcome this problem, some scholars 46-51 developed a series of quasi-static elastoplastic contact models to describe the elastoplastic contact events using the piecewise function, including the elastic, elastoplastic, and full plastic phase. Since these elastoplastic contact models depend entirely on the constitutive relation of contact bodies, they can accurately reflect the relationship between the force and displacement in the elastoplastic phase.
Whereupon, these quasi-static elastoplastic models 16 had been used in calculating the solitary wave propagation of elastoplastic granular chain 17 in past decades. 16,17,[52][53][54] The conclusions obtained by the experimental data and numerical analysis illuminated that the elastoplastic deformation between the particles plays an essential role in propagating solitary waves in the granular chain. 17,54 In this process, the energy dissipation during contact is represented by the enclosed area composed of the loading and unloading path 55,56 However, since the contact behavior in the elastoplastic phase follows the different constitutive relation between the compression and recovery phase, the maximum contact deformation and residual deformation must be identified and saved at each impact for the sake of preparing for the following contact behavior. The situation makes the simulation process of the contact event complicated compared with the dynamic dashpot models. That is why the quasi-static elastoplastic contact models were not widely used in the EDEM commercial software.
In light of the research background above, two kinds of models, including quasi-static elastoplastic contact models 16,17,52,54,57,58 and dynamic dashpot models [39][40][41][42] without the numerical singular issue, are applied to simulate the solitary wave propagation in the granular chain. However, both of them confront with some issues. As for the dynamic dashpot models, 43 they can represent the energy dissipation by the damping factors and do not need to distinguish the compression and recovery phase in calculating the impact behavior. This approach extremely simplifies the simulation strategy of collision events and benefits for programming. Nevertheless, they cannot reflect the elastoplastic contact characteristic because of the Hertz contact stiffness. 41 On the contrary, although the quasi-static elastoplastic contact models can accurately describe the relationship between the loading and displacement in the elastoplastic contact phase, they made the calculation strategy of elastoplastic contact behavior in the granular chain more complicated. 16,52,54,57 It needs to be emphasized that the elastoplastic collision in the granular chain depends on the constitutive relation of contact bodies, 46   is the critical plastic deformation. ε is the dimensionless parameter that corresponds to a geometric relationship when the pressure on the contact surface approximately approaches uniformity. Its value is within the scope from 13 to 20 when the contact spheres are made of the same materials. 43 Moreover, the specified form of the coefficient in this contact model is expressed as = ( + ln( /2)) + , = (1 + ln( /2)) − 2 ln( /2) , where ψ is the dimensionless parameter that corresponds to a ratio between the Brinell hardness and the yielding strength of the material, which is in the range from 2.6 to 3.0. 43 During the unloading phase, the constitutive relation between the contact force and displacement is expressed as where δ max is the maximum contact deformation and δ r is the residual contact deformation.
Since both elastoplastic and plastic deformation happened before the restitution phase, the curvature radius is larger than before. Therefore, when the contact behavior aborts at the elastoplastic compression phase, the radius of curvature of the contact body can be written as where R ep e is the radius of curvature after impact in the elastoplastic phase and F max is the maximum contact force.
When the contact behavior ends at the plastic compression phase, the radius of curvature of the contact body can be written as where R p e is the radius of curvature after impact in the plastic phase.  models suffer from a seriously numerical problem when the initial impact velocity is very small or equal to zero. 32 The damping force in the existing dynamic dashpot models is infinite when the initial impact velocity equals zero, which is impossible and violates the realistic physical meaning. That is mainly because the denominator of damping force in most models includes the initial impact velocity. 32 The other dashpot models 41

Mδ Dδ Kδ
where M is the mass of the contact body, D is the damping coefficient, and K is the stiffness coefficient of spring.
The solution of the underdamped system is expressed as where A is the amplitude, ϕ is the phase angle, ω is the frequency, the damping, and ξ D Mω = /2 .
The amplitude and phase angle can be determined from the initial condition Equation (7) can be rewritten as The duration time of the entire contact process is expressed as  When an entire contact process is over, the deformation in Equation (9) is equal to zero, the deformation velocity is written as according to Equation (9) Based on the definition of Newton's coefficient of restitution (COR), one can obtain the following equation where c r is the COR and D is the damping coefficient.
Therefore, according to Equation (6), the other new dynamic dashpot model in the elastoplastic or plastic phase can be written as where K p is the linearized contact stiffness in the elastoplastic or plastic phase based on the ML model in Figure 1.
When the contact behavior is in the elastic phase, the contact stiffness K is assumed as K K δ = e 1 2 , 59 which is still linear with the contact deformation. Namely, the dynamic dashpot model can also be obtained according to Equation (6): It is noteworthy that the new dynamic dashpot model is almost the same as the dashpot model used in the EDEM software, which is expressed as 41 Equations (13) and (14) (13) and (14) can evaluate the dynamic collision behavior and accurately calculate the dynamic elastoplastic collision event.

| HORIZONTAL GRANULAR CHAIN
To exhibit the merit of the new dashpot model in calculating dynamic elastoplastic contact behavior, Figure 3 shows the schematic diagram of the experimental setup that is established by Daraio et al. 18 for a monodisperse granular chain of 70 stainless steel particles assembled horizontally, the particle is made of stainless steel 316. The particle size and material properties in the chain are identical to each other, as shown in Table 2.
The striker impacts along the ramp activated the solitary waves. pening, and then, the contact force between the particles can be calculated using Equation (14) in the elastic phase and Equation (13) Figure 4A is smaller than the blue error bar between the experimental data and the ML model in Figure 4B. This result means the maximum contact force obtained using the new dashpot model is more closed to the experimental data compared with the ML model in the elastoplastic phases in Figure 4C.
Significantly, in Figure 5B, the loading paths of the different particles are the same, but the unloading paths are different from each other. Namely, the ML model can obtain the residual deformation according to the unloading path. However, in Figure 5A, the hysteresis loop depends on the relative deformation velocity direction during impact. Hence, the compression and recovery paths are discrepant for each particle. Meanwhile, the energy dissipation obtained by the ML model is represented by the difference between the loading and unloading path in Figure 5B, which is homogenized to the hysteresis loop in the new dashpot model in Figure 5A. Although these two contact models, including the ML model and the new dashpot model, used different approaches to calculate the elastoplastic contact behavior, the maximum contact force and the maximum contact deformation can keep harmonizing with each other. Therefore, the relative impact velocity between the particles calculated using the ML model is closed to the results obtained using the new dynamic dashpot model in Figure 6.
As the solitary wave propagates forward ceaselessly, the kinetic energy is gradually dissipated, the contact behavior comes into the elastic  In Figure 7C, the solitary wave propagation obtained using the new dashpot model in Figure 4A is almost identical to the results using the dashpot model used in EDEM software in Figure 7A, which can keep consistent with the experimental data. However, in calculating the maximum contact force between the particles, the blue error bar between the new dynamic dashpot model and experimental data is smaller than the red one between the dashpot model used in the EDEM software and experimental data in Figure 7B. The absolute error percentage of the peak value of the solitary wave can be seen in  Figure 8. This conclusion illuminates that the new dynamic dashpot model possesses higher fidelity than the dashpot model used in the EDEM software. It can accurately describe the relationship between the contact force and deformation in three different contact phases, i.e., elastic, elastoplastic, and plastic phases.

| CONCLUSION
The motivation of this investigation lies in that most existing dynamic dashpot models suffer from a serious numerical issue when the initial impact velocity is very small or equal to zero. Moreover, since the dynamic simulation process of the elastoplastic collision behavior is complicated when using the static elastoplastic contact

CONFLICT OF INTEREST
The authors declare that there are no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.