The establishment of the normal contact force model for a one‐dimensional sphere chain subjected to impact load

The normal contact force determines the behavior of a particle system. To investigate the normal contact force in a one‐dimensional sphere chain subjected to impact load, by comparing the simulation results of the existing typical normal contact force models embedded in the discrete element program, an improved normal contact force model was proposed in this paper. The improved model consists of two parts: the Cundall model for loading and the Daniel model for unloading. Moreover, a systematic test was designed to verify the accuracy and applicability of the improved model. The results showed that the calculated contact force curves agree well with the experimental results. Furthermore, the improved model is implemented in the solution algorithm without need for complex numerical methods and parameters fitting, leading to more efficient simulations.

terms. For instance, Tsuji et al. 10 suggested a damping term proportional to the fourth root of the sphere displacement. The damping term was modified to be proportional to the 1.5 power of the displacement by Hu et al. 11 and 1-1.5 power of the displacement by Hunt et al. 12 Kruggel-Emden et al. 13 proposed a nonlinear damping term that is exponential to the displacement of spheres, and the order was suggested to be 0.5-1. Zdancevicius et al. 14  Kuwabara and Kono 29 studied the internal friction or viscosity of solids, and proposed a dissipative term for the normal contact force. Brilliantov et al. 30 and Zheng et al. 31 assumed that the displacement field of a viscoelastic sphere is equal to that of an elastic sphere or consistent with that of the static case, and formulated different expressions of the dissipative term. Brilliantov et al. 32 developed a perturbation approach, which makes the derivation of dissipative force more rigorous without quasi-static approximation.
While many contributions have been made for developing and improving normal contact-force model, a physically realistic model for the sphere chain subjected to impact load has not been yet established. In this paper, five typical normal contact force models are implemented in a unified DEM algorithm used in the simulations of the collisions of sphere particles. An improved normal contact-force model was proposed by comparing simulation results of different existing models. A systematic test method is also proposed to investigate the normal force of the one-dimensional sphere chain under the impact load and verify feasibility and applicability of the improved model for onedimension sphere chain. In particular, three typical existing normal contact models and the improved model for a onedimensional sphere chain subjected to impact load were verified by comparing test and numerical results.

| DEM ALGORITHM DESIGN
The DEM program was developed using a modular structure on the platform of VC++. The modular design allows dividing a complex program into several independent simple programs.
The program used in this investigation is designed to include three main modules: pre-processor, main processor, and visualization post-processor. The pre-processing module reads the geometric characteristics, mechanical parameters, coordinates of the points of interest, and granular material's parameters to generate a particle system based on the input data. The main processor module represents the solver. First, the contact force is calculated at every time step according to the applied force.
Then, every particle state is iteratively updated at acceleration, velocity, and displacement levels. Finally, particle-contact information, including contact force at every time step, is saved to have the data used in the visualization. The visualization module uses the OpenGL graphics library, which has a visual output function. During the simulations, according to the contact information, the post-processor can be used for online motion visualization and plotting. The detailed algorithm is shown in Figure 1. Five typical normal contact force models (as listed in Table 1) are implemented in the DEM computer program for the purpose of comparative analysis.

| NORMAL CONTACT FORCE MODELS' COMPARISON
The collision of two sphere particles under impact was simulated by the finite element method (FEM) and DEM using different normal contact force models. ANSYS software was used for FEM modeling, and 164 solid elements were selected to establish sphere particles (as shown in Figure 2). Since a spherical particle is prone to mesh distortion during mesh section, a 1/8 model is established and meshes are divided first, and then the whole particle is generated by coordinate axis mirroring. The contact type is automatic surface-tosurface contact and the particle system is constrained laterally, that is, only normal deformation is allowed. Additionally, the bottom of the particle system is modeled as a fixed boundary. A triangular pulse load with a 0.1 s loading duration and a 10 N peak value was applied at the central vertex of the particle system. The material of spherical particles is cement, and the diameters of the sphere are set as 30, 40, and 50 mm. Furthermore, the material properties were tested using a TAJW-2000 electrohydraulic servo testing machine and a rock density tester in the laboratory ( Table 2).  When t = 0.1 s, the stress at the contact point almost disappeared approaches zero. The stress at the center of the sphere is the largest, and the stress wave is diffused throughout the sphere particle. When t = 0.18 s, the two particles are disconnected. Figure 4 shows contact force-deformation curves of three particles with different sizes under the condition of the same material and impact load. The results were obtained using FEM and DEM approaches with different normal contact force models. Table 3 lists the predicted maximum contact force and residual deformation of different particle sizes using different DEM models and FEM approach.
By comparing the simulation curves of the Hertz and FEM models, it can be seen that they intersect at the loading stage in Figure 4A, and the intersection point is near the critical point of loading and unloading of the FEM curve in Figure 4B, and at the unloading stage in Figure 4C. This result suggests that the mechanical response becomes weaker with increasing particle size, which is consistent with the real situation. Since the Hertz theory is an elastic T A B L E 1 Equations used for normal contact force models Models Expression Parameters   calculating effective contact radius. Therefore, effective contact radius is a dynamically changing value in the particle separation phase, which is consistent with actual situation.

F I G U R E 4
Comparison of contact force-deformation curves with different methods for different particle sizes: (A) 3 cm particle diameter; (B) 4 cm particle diameter; and (C) 5 cm particle diameter T A B L E 3 Comparison of the maximum contact force P max and residual deformation δ p 3 cm particle diameter 4 cm particle diameter 5 cm particle diameter Models He model, the effect of particle size on the force-displacement curves from different models is consistent. The expression of the equivalent radius is R R R 1/ = 1/ + 1/ 1 2 . According to the expressions of contact models, the contact deformation decreases under the same load, while the equivalent radius increases with increasing radii R 1 and R 2 . In general, the larger the particle radius, the stronger the deformation resistance of spherical particle. As shown in Table 3, the maximum contact force decreases with increasing particle diameters. For the Cundall model, the values of stiffness and the damping coefficient highly influence the calculated results as compared to particle size.

| IMPROVEMENT AND VALIDATION OF THE NORMAL CONTACT FORCE MODEL
The advantages and disadvantages of each normal contact model were analyzed by simulating sphere particle collision in Section 3. To Based on the above analysis, an improved normal contact force model is proposed and can be expressed as where definitions of the variables can be found in Table 1.  The cement mortar balls are made using an aluminum mould with two separate hollow hemispheres. The balls are polished to reduce the influence of friction between the balls and the inner wall of the cylinder.
In the cylinder wall, 4 mm holes are drilled to connect the piezoelectric sensors. The piezoelectric film sensors, whose effective area includes a central circular film with a diameter of 20 mm and a thickness of 0.25 mm, are placed between particles to measure the contact force. The contact point of two adjacent balls is at the center of the film, which can be regarded as the common tangent plane. The weight of the sensors is much lower than the weight of the ball, so its gravity effect can be ignored.  The error of contact-force peaks in the model measured with respect to that of the experiment is calculated (see Table 5). The  shows that the simulated waveform obtained by the improved model agrees with the measured waveform in the duration and the concentration degree of the main waveform. The average error of the waveform peak is only 10%, which is the smallest error among the models examined in this paper.
Finally, it is worth noting that the improved model was implemented in the DEM software without the need for using complex numerical methods and parameter fitting, which demonstrates the generality and ease of implementation of the improved model.

ACKNOWLEDGMENTS
This study was supported by the National Natural Science Foundation of China (Nos. 51874118, 51778211).