Calibration and Experiments on the Parameters of the Bonding Particle Model of Eucommia ulmoides Oliver Samara Based on the Discrete Element Method

: The calibration of the discrete element model of EUO samara was conducted in response to the lack of discrete element simulation models and parameters in the development of mechanical hulling technology and equipment of EUO samara. The EUO samara was modeled based on the Hertz–Mindlin with bonding model, and its relevant parameters were measured by physical experiments. Among them, maximum shear force was used as the evaluation index, virtual calibration experiments were conducted for the bonding parameters by the single-factor experiments, and the two-level factorial experiment, the steepest climb experiment, and the Box–Behnken response surface experiment were also used. The results showed that the relative error between the simulated and measured maximum shear force is 0.93%; the model and parameter calibration results established by this research can be used for discrete element simulation research, which is of guiding signiﬁcance for the research and development of hulling equipment of EUO samara.


Introduction
Eucommia ulmoides Oliver (EUO) is a traditional medicinal resource in China; it is widely grown in many provinces in China, including Shaanxi, Sichuan, Yunnan, Guangxi, Zhejiang, and Guizhou [1].As a part of EUO, the EUO samara shell is rich in EUO rubber, and its kernel is rich in EUO seed oil, which is widely used in the fields of national health, industrial raw materials, and medical devices and materials, and it has a high development value [2][3][4].At present, the demand for EUO samara is urgent, but its industrialization is low.Generally, this is because EUO samara generally grows on arbor trees of 15-20 m in height, and most of these EUO trees are located in mountainous areas, making harvesting difficult.More importantly, the low degree of mechanization in the hulling and separation of EUO samara restricts its industrialization [5].In order to improve the mechanization of hulling and separation of EUO samara, it is crucial to optimize the relevant equipment.The optimization process usually includes bench and simulation experiments.However, for agricultural equipment, optimization by simulation tests is easier to achieve and more economical [6].
Recently, the application of the discrete element method (DEM) and its simulation software in the field of agriculture has provided a new way to study the contact characteristics between agricultural materials and machinery components, which has contributed to the development of agricultural machinery [7,8].In general, parameter calibration and discrete element modeling of the target material are essential aspects of discrete element simulation.The agricultural material parameters mainly consist of geometric parameters, contact property parameters, and mechanical property parameters [9].Most of the material properties' parameters can be directly measured or calculated by bench experiments, such as geometry, moisture content, density, Poisson's ratio, modulus of elasticity, etc. [10].However, contact property parameters, such as recovery of the coefficient and friction of the coefficient, and some mechanical property parameters are generally difficult to measure directly and often need to be calibrated by bench experiments combined with simulation experiments [11,12].
Normally, DEM calibration studies center on contact property parameters, and many researchers have completed calibrations of the coefficient of friction and coefficient of restitution for materials such as corn [13][14][15], wheat [16,17], rice [18,19], peanuts [20], soybeans [21], soil [22], and urea [23].The calibration of these contact characteristics is based on the Hertz-Mindlin (no slip) model.The material model is considered as a whole in DEM, and the calibrated model can be used for simulation analysis and optimization design of transport, seeding, and other processes [24,25].However, in processes that require the destruction of material such as hulling, a discrete element model that can be broken is required; this is when a bonding particle model (BPM) based on the Hertz-Mindlin with bonding principle is required [26].The BPM uses a dense accumulation of circular particles of arbitrary size that are bonded together at the point of contact and simulated with DEM [27].This model has been used to simulate rocks as a direct particle model approach to simulate inelastic deformation and fracture of rocks [28].The BPM of rocks directly simulates the behavior of rocks like a complexly shaped particle cemented granular material, where both the particles and the cement are distortable and may be destroyed [29].In principle, the principles of BPM can explain all of the aspects of the relevant behavioral mechanisms.The system is very accurate in predicting sudden behavior changes between particles and corresponds well to real application scenarios.It offers the possibility to study the microscopic mechanisms that generate these complex macroscopic behaviors and to predict them in DEM [30].
The particle-based part of the force-displacement behavior at each contact is described by six parameters, and a detailed description of the contact force-displacement calculation is provided in Figure 1a [31].According to the theory of beams, the formulae for the maximum tensile and maximum shear stresses acting on the parallel bonded edges are as follows [32]: If the maximum tensile stress exceeds the tensile strength (σ max ≥ σ c ) or the maximum shear stress exceeds the shear strength (τ max ≥ τ c ), the parallel bonds break and are removed from the model along with their accompanying forces, bending moments, and stiffness.Based on these principles, the model can also be used to reveal the mechanism of action when the material and mechanical components come into contact, which can guide the development and design of the device.Zhang et al. [33] developed a discrete element model of water chestnut through reverse engineering techniques and calibrated the bonding parameters of the model, using the results to guide the design of water chestnut peeling and dicing equipment.Li et al. [10] calibrated the bonding parameters of corncobs by combining physical bending experiments and simulated bending experiments using a segmental modeling approach for corncobs.Liao et al. [34] addressed the simulation of the chopping aspect of mechanical harvesting of fodder rape by conducting a simulation experiment on the bending damage of fodder rape stalks and conducting a calibration of the bonding parameters of fodder rape to guide the optimization of chopping equipment.EUO samara has a high exploitation value, and the relative relationship between the hulling equipment and EUO samara in the exploitation process should not be ignored.To date, few studies have also been reported on the discrete element modeling and relevant parameter calibration of EUO Samara.Therefore, in this paper, geometric parameters and contact property parameters of EUO samara were obtained by physical experiments, the shape profile of EUO samara was obtained by three-dimensional (3D) scanning technology, and the discrete element model of EUO samara was established by importing DEM discrete element simulation software.Moreover, the maximum shear force of the blade on EUO samara was obtained by a Texture Analyzer.Using this as a reference, a virtual calibration experiment of the bonding parameters of the discrete element model was conducted to obtain the bonding parameters of EUO samara that can be used in the discrete element simulation.

Experimental Materials and Geometric Parameters
EUO samara from the EUO arborescent tree was used in the experiment as the material, which has a rough and dark brown surface and is usually flattened with a prismatic or linear structure [35,36].In this study, 100 randomly selected samples of EUO samaras were measured by using vernier calipers, with dimensions ranging from 29.2 mm to 39.6 mm in length, 8.8 mm to 12.5 mm in width, 1-3 mm in thickness, and a fruit shape index of 2.6 to 3.8.The average moisture content was determined to be 7.51% by using the drying method.
The density  is measured by the specific gravity bottle method and is calculated as where A is the mass of the EUO samara, B is the mass of the specific gravity bottle filled with water, C is the total mass of the remaining water and samara, and  is the density of the water [37].
The weighing of the masses in the above equation was carried out using a high-precision electronic balance (0.0001 g) and randomly selected 200 EUO samaras for 10 replicate experiments with a density measurement of 0.67 g/cm 3 .EUO samara has a high exploitation value, and the relative relationship between the hulling equipment and EUO samara in the exploitation process should not be ignored.To date, few studies have also been reported on the discrete element modeling and relevant parameter calibration of EUO Samara.Therefore, in this paper, geometric parameters and contact property parameters of EUO samara were obtained by physical experiments, the shape profile of EUO samara was obtained by three-dimensional (3D) scanning technology, and the discrete element model of EUO samara was established by importing DEM discrete element simulation software.Moreover, the maximum shear force of the blade on EUO samara was obtained by a Texture Analyzer.Using this as a reference, a virtual calibration experiment of the bonding parameters of the discrete element model was conducted to obtain the bonding parameters of EUO samara that can be used in the discrete element simulation.

Experimental Materials and Geometric Parameters
EUO samara from the EUO arborescent tree was used in the experiment as the material, which has a rough and dark brown surface and is usually flattened with a prismatic or linear structure [35,36].In this study, 100 randomly selected samples of EUO samaras were measured by using vernier calipers, with dimensions ranging from 29.2 mm to 39.6 mm in length, 8.8 mm to 12.5 mm in width, 1-3 mm in thickness, and a fruit shape index of 2.6 to 3.8.The average moisture content was determined to be 7.51% by using the drying method.
The density ρ is measured by the specific gravity bottle method and is calculated as where A is the mass of the EUO samara, B is the mass of the specific gravity bottle filled with water, C is the total mass of the remaining water and samara, and ρ 1 is the density of the water [37].
The weighing of the masses in the above equation was carried out using a highprecision electronic balance (0.0001 g) and randomly selected 200 EUO samaras for 10 replicate experiments with a density measurement of 0.67 g/cm 3 .

Mechanical Property Parameters
The mechanical property parameters include Poisson's ratio, modulus of elasticity, and intrinsic mechanical parameters.
Poisson's ratio µ, also called the transverse deformation coefficient, is the ratio of transverse positive strain to axial positive strain when the EUO samara is subjected to unidirectional tension or compression, and the µ is the elastic constant reflecting the transverse deformation of the EUO samara.The Poisson's ratio of juniper winged fruits was measured experimentally using the definition method, which was calculated as where ε x is the samara transverse strain and ε y is the samara axial strain; ∆L is the absolute samara transverse deformation and L is the samara transverse original length; ∆H is the absolute samara axial deformation and H is the samara axial original length [38].The Poisson's ratio was determined by the uniaxial compression experiment of the Texture Analyzer (TA.XTC-1803061210), as shown in Figure 2a, with a flat indenter at the end of the Texture Analyzer and a loading rate of 0.5 mm/s.The experiment was repeated 20 times, and the Poisson's ratio of the EUO samara was obtained as 0.403.

Mechanical Property Parameters
The mechanical property parameters include Poisson's ratio, modulus of elasticity, and intrinsic mechanical parameters.
Poisson's ratio , also called the transverse deformation coefficient, is the ratio of transverse positive strain to axial positive strain when the EUO samara is subjected to unidirectional tension or compression, and the  is the elastic constant reflecting the transverse deformation of the EUO samara.The Poisson's ratio of juniper winged fruits was measured experimentally using the definition method, which was calculated as where  is the samara transverse strain and  is the samara axial strain; ∆ is the absolute samara transverse deformation and  is the samara transverse original length; ∆ is the absolute samara axial deformation and  is the samara axial original length [38].The Poisson's ratio was determined by the uniaxial compression experiment of the Texture Analyzer (TA.XTC-1803061210), as shown in Figure 2a, with a flat indenter at the end of the Texture Analyzer and a loading rate of 0.5 mm/s.The experiment was repeated 20 times, and the Poisson's ratio of the EUO samara was obtained as 0.403.The modulus of elasticity E is a measure of the magnitude of the ability of the EUO samara to resist elastic deformation and was determined by a uniaxial compression experiment with a universal testing machine (DDL1003050301), with a flat indenter selected at the end of the universal testing machine and a loading rate of 0.5 mm/s.The experiment was repeated 20 times, and the modulus of elasticity of the EUO samara was obtained as 19.98 MPa.The modulus of elasticity E is a measure of the magnitude of the ability of the EUO samara to resist elastic deformation and was determined by a uniaxial compression experiment with a universal testing machine (DDL1003050301), with a flat indenter selected at the end of the universal testing machine and a loading rate of 0.5 mm/s.The experiment was repeated 20 times, and the modulus of elasticity of the EUO samara was obtained as 19.98 MPa.

Coefficient of Restitution
The coefficient of restitution e characterizes the ability of a material to recover after a collision and is only relevant to the material itself.Its value is the ratio of the normal relative separation velocity to the normal relative approach velocity of the two objects at the point of contact before and after the collision.This study involves the contact between the EUO samara and the contact between the EUO samara and stainless steel, and the parameter is measured with the aid of a high-speed camera (I-SPEED TR 4 GB MONO05020701) and coordinate paper [39].
Drop experiment was used to measure the coefficient of restitution e 1 between the EUO samara and stainless steel, as shown in Figure 2c, the EUO samara was lifted to =30 mm so that it was released in free-fall motion and collided with the stainless steel plate radially, then the rebound height H 1 was recorded after the rebound.The formula for calculating the coefficient of restitution e 1 between the EUO samara and stainless steel were The experiment was repeated 20 times, and the coefficient of restitution e 1 was measured to be 0.335 between the EUO samara and stainless steel.
The drop experiment was used to measure the coefficient of restitution e 2 between the EUO samaras.A binder was used to bond the EUO samaras into seeds on a stainless steel plate for the drop experiments, which were repeated 20 times as above, and the coefficient of restitution e 2 between the EUO samaras were measured to be 0.373.

Coefficient of Static Friction
The coefficient of static friction was determined by the inclined plane method.As shown in Figure 2b, a digital inclinometer was used to measure the inclination angle of the inclined plane during the experiment.
To measure the coefficient of static friction between the EUO samara and stainless steel, we used a binder to bond the stainless steel plate to the inclined surface, then put the EUO samara on the stainless steel plate, slowly turned the angle adjuster to increase the inclined angle of the inclined surface until the EUO samara was observed to slide on the stainless steel plate, then we stopped turning, and recorded the angle α displayed by the digital display inclination measuring instrument [40].The formula for calculating the coefficient of static friction γ 1 between the EUO samara and stainless steel is The experiment was repeated 20 times, and the coefficient of static friction between the EUO samara and stainless steel was measured at 0.442.
To measure the coefficient of static friction γ 2 between the EUO samaras, the stainless steel plate was simply replaced with one bonded with EUO samaras stock.The experiment was repeated 20 times, and the coefficient of static friction between the EUO samaras was measured to be 0.676.

Coefficient of Rolling Friction
A multi-functional friction and wear test system (MFT-500003050903) was used for the measurement of the coefficient of rolling friction (as shown in Figure 2d).The coefficient of rolling friction was measured by pair grinding of the materials to be measured, and the experiment was repeated 20 times, and the coefficient of rolling friction γ 3 was measured to be 0.254 between the EUO samara and the stainless steel.
The coefficient of rolling friction γ 4 between the EUO samaras was obtained as 0.05 with reference to similar materials [41].

Bonding Parameters and Calibration Experiments 2.4.1. Bonding Parameters
Shear experiments were carried out on the EUO samara by using a Texture Analyzer to obtain a realistic reference value reflecting the bonding parameters, and a square blade was selected at the end of the Texture Analyzer for the shear experiments.The end of the Texture Analyzer was run at a speed of 0.5 mm/s during the experiments.Due to the brittleness of the EUO samara, the force curve of the shearing process reached its peak when the EUO samara fractured, and there was no fixed rule for the location of the peak between different individuals, so only the maximum value of the shearing force was selected as the evaluation index in this study.The experiment was repeated 20 times, and the maximum shear force F a of the blade on the EUO samara was 43.2 N.This was used as a physical reference value for the bonding parameters for the virtual calibration experiment.

The Single-Factor Experiments
According to the definition of the Hertz-Mindlin bonding model, the breakage of the bonding bond between particles is related to the particle contact radius x 1 , the normal stiffness per unit area x 2 , the shear stiffness per unit area x 3 , the critical normal stress x 4 , the critical shear stress x 5 and the bonded disk radius x 6 .Regarding reference to the simulation parameters of materials with similar morphological and mechanical properties, a preliminary range of discrete element simulation bonding parameters for EUO samara was set as shown in Table 1 [10,33,34], and the influence of each factor on the F a was analyzed by single-factor experiments and the parameter range was determined.Single-factor simulation experiments were carried out in DEM discrete element software with F a as the evaluation index; factors x 1 and x 6 were taken as 5 points on average, factors x 2 -x 5 had a larger range, and in order to take into account different orders of magnitude, points were taken for each order of magnitude endpoint and intermediate position, as shown in Table 2, the fixed factors were taken as intermediate levels, x 1 is 0.8 mm, x 2 is 4.0 × 10 7 N•m −3 , x 3 is 4.0 × 10 7 N•m −3 , x 4 is 4.0 × 10 6 Pa, x 5 is 4.0 × 10 6 Pa and x 6 is 0.8 mm.The two-level factorial experiment can quickly screen out the factors that have a significant impact on the test index during many influencing factors and determine the range of each parameter based on the results of the single-factor experiments.A total of 16 groups of two-level factorial experiments are designed.

The Steepest Climb Experiments
The steepest climb experiment is carried out on the factors screened for significance in the two-level factorial experiments to quickly determine the neighborhood of the optimum value, with the non-significant factors taken at a fixed level in the single-factor experiments.

The Box-Behnken Response Surface Experiments
To obtain the optimal simulation parameters, the response surface experiments were designed according to the Box-Behnken principle.The significant influencing factors x 1 , x 2 , x 3 and x 6 screened by the two-level factorial experiments, and the factors were narrowed down according to the results of the steepest climb experiments.Moreover, the levels in experiments 4 and 6 are the upper and lower limits for the response surface experiments, with four replications set at the central level, the response value of F a of EUO samara was used to obtain the regression equation with the significant influencing factor as the variable.Meanwhile, the relative error e between the measured shear force and the simulated shear force is calculated as follows: where e is the relative error, F a is the maximum measured shear force and F b is the maximum shear force in the simulation experiments.

EUO Samara Model
The shape of the EUO samara is mostly irregularly surface (Figure 1b), and it is difficult to model its actual shape features by conventional modeling means.In order to obtain an accurate profile of the EUO samara and improve the accuracy of the simulation, this study used the 3D scanner (HandySCAN BLACKTM-Elite0321111) to scan the profile of the EUO samara to obtain accurate scan data of the EUO samara, imported into Geomagic Wrap to process the 3D scan data, remove abnormal scan points, fill in the gaps and adjust the coordinates to the center point (Figure 1c), then generated the stl model file and imported into DEM to build a discrete element model of the EUO samara with bonded bonds, as shown in Figure 1d.
The EUO samara discrete element modeling process is as follows: The geometry of the box was established in DEM software (version EDEM2020, ALTAIR, Troy, MI, USA, 2019), and the processed EUO samara model was imported into its internal area and set as virtual type (Figure 3a).A large number of particles were generated inside the box to fill the model.After filling, the model was set as physical type.At the same time, change the box type to virtual, and the excess particles will leave the calculation domain and disappear automatically under the action of gravity (Figure 3b).At this point, the internal particle coordinates of the EUO samara model were derived, and the derived coordinates were used to generate a discrete element model of the EUO samara with bonding bonds by particle replacement (Figure 3c,d).In the discrete element model, the smaller the particle radius is, the closer the simulation result is to the real situation, but the amount of calculation will also increase dramatically.In this case, the particle radius was set as 0.5 mm in combination with the blade thickness in actual shear and the shape parameter of EUO samara.Finally, the number of particles in the EUO samara model was 444, and a total of 2423 bonding bonds were generated.The average number of bonds in each particle was more than 5.5, and the bonding was sufficient.

Shear Model
The blade and base models were created in SolidWorks and saved in step mode, with all dimensions consistent with the actual experiments (Figure 3e).After the particle replacement was completed, they were imported into DEM, while the coordinates were adjusted to ensure that the shear angle was consistent with the actual experiments, and the blade loading speed was set at 0.5 mm/s, with the direction vertical downwards, for the shear simulation experiments (Figure 3f).total of 2423 bonding bonds were generated.The average number of bonds in each particle was more than 5.5, and the bonding was sufficient.

Shear Model
The blade and base models were created in SolidWorks and saved in step mode, with all dimensions consistent with the actual experiments (Figure 3e).After the particle replacement was completed, they were imported into DEM, while the coordinates were adjusted to ensure that the shear angle was consistent with the actual experiments, and the blade loading speed was set at 0.5 mm/s, with the direction vertical downwards, for the shear simulation experiments (Figure 3f).

Calibration Results Validation
To test whether the calibrated parameters can be used in a simulation study of the hulling of EUO samaras, simulations were carried out using the optimal parameters solved by regression equations and compared with measured values to verify the accuracy of the discrete element model and parameter calibration of EUO samaras.
The above simulation experiments with a simulation time step of 5% of the Rayleigh time step, namely, 8.08037 × 10 −7 s.The data storage time range was 0.01 s, and the total simulation time was 26 s, the grid size was 1 mm.The DEM simulation was run on a computer containing an eight-core processor (Inter ® Core TM i5-10400F CPU @ 2.90 GHz) with 16 GB of memory (RAM).

The Single-Factor Experiment Results and Analysis
The results of the single-factor experiments are shown in Figure 4.The degree of influence of each factor on  is different, among which the simulation results do not change significantly when the critical normal stress  (Figure 4d) and critical shear stress  (Figure 4d) are taken to different values, indicating that their influence on the shear force is small, and the results of these two groups of experiments are between 40-50 N,

Calibration Results Validation
To test whether the calibrated parameters can be used in a simulation study of the hulling of EUO samaras, simulations were carried out using the optimal parameters solved by regression equations and compared with measured values to verify the accuracy of the discrete element model and parameter calibration of EUO samaras.
The above simulation experiments with a simulation time step of 5% of the Rayleigh time step, namely, 8.08037 × 10 −7 s.The data storage time range was 0.01 s, and the total simulation time was 26 s, the grid size was 1 mm.The DEM simulation was run on a computer containing an eight-core processor (Inter ® Core TM i5-10400F CPU @ 2.90 GHz) with 16 GB of memory (RAM).

The Single-Factor Experiment Results and Analysis
The results of the single-factor experiments are shown in Figure 4.The degree of influence of each factor on F a is different, among which the simulation results do not change significantly when the critical normal stress x 4 (Figure 4d) and critical shear stress x 5 (Figure 4d) are taken to different values, indicating that their influence on the shear force is small, and the results of these two groups of experiments are between 40-50 N, which are close to the F a , indicating that the fixed level of each factor is close to the correct parameter and the parameter range is chosen reasonably.
With the increase in particle contact radius x 1 , the shear force also gradually increases when x 1 is between 0.7-0.8mm, the simulation result is closest to the measured value (Figure 4a).The shear force gradually increases as the normal stiffness per unit area x 2 increases, and the simulation result is closest to the measured value when x 2 is between 2.0 × 10 7 N/m 3 -4.0 × 10 7 N/m 3 (Figure 4b).When the shear stiffness per unit area x 3 increases, the shear force also increases, and the simulation result is closest to the measured values as the x 3 is between 2.0 × 10 7 N/m 3 and 4.0 × 10 7 N/m 3 (Figure 4c).As the bonded disk radius x 6 increases, the shear force also increases, and the simulation result is closest to the measured values when x 6 is between 0.7-0.8mm (Figure 4f).

The Steepest Climb Experiments Results and Analysis
According to the results of the two-level factorial experiments, the factors x 1 , x 2 , x 3 and x 6 , which have a significant effect on the maximum shear, were selected for the steepest climb experiments.The non-significant factors x 4 and x 5 were set as fixed levels in the single-factor experiments.The experiment protocols and results are shown in Table 5.With the increase in each factor level, the shear force of the blades on the EUO samaras increased, among which the maximum values of simulated shear force in experiment five and experiment six were 40.9 N and 61.2 N, respectively.The results of experiment five were close to F a , which was not conducive to the subsequent development of response surface experimental design.Therefore, it was appropriate to extend the range of factors and use the factor levels from experiments four and six as the low and high levels for the subsequent Box-Behnken experiments.

The Box-Behnken Response Surface Experiments Results and Analysis
In order to obtain the optimal combination of simulation parameters for the discrete element model of EUO samara, Box-Behnken experiments were conducted on the significance factors according to the results of the two-level factorial experiments and the steepest climb experiments, and the experimental protocols and results are shown in Table 6.The regression results were fitted to obtain the regression equation with F a of the blade on EUO samara as the response value and x 1 , x 2 , x 3 and x 6 as variables, and the regression equation is shown as follows: The Box-Behnken response surface experiments results were analyzed by variance in Design-Expert software, and the results are shown in Table 7.The model has a good fit because of the high coefficient R 2 of 0.9503.The model is significant, and there are no other major factors affecting the response value of this experiment.p < 0.01 for factors x 1 , x 2 , x 3 and x 6 , with a highly significant effect on the shear force, p > 0.05 for the interaction of factors x 1 , x 2 , x 3 and x 6 and the squared terms of factors x 1 , x 2 , x 3 and x 6 , which have a non-significant effect on shear force.Response surface diagrams of the effect of the interaction of the factors on F a are obtained separately from the regression equation, and the relevant results are shown in Figure 5. because of the high coefficient R 2 of 0.9503.The model is significant, and there are no other major factors affecting the response value of this experiment.p < 0.01 for factors  ,  ,  and  , with a highly significant effect on the shear force, p > 0.05 for the interaction of factors  ,  ,  and  and the squared terms of factors  ,  ,  and  , which have a non-significant effect on shear force.Response surface diagrams of the effect of the interaction of the factors on  are obtained separately from the regression equation, and the relevant results are shown in Figure 5.When x 2 , x 3 and x 6 are fixed, the shear force increases with the increase in x 1 , and the increasing trend is obvious (Figure 5a-c), indicating that x 1 is a comparatively critical factor affecting the magnitude of the shear force.The reason may be that x 1 is a condition to determine whether the contact model is triggered, and a large particle contact radius is more likely to trigger the Hertz-Mindlin with bonding model [42], while the Hertz-Mindlin with bonding model is used as the basis for this study, accordingly a large x 1 will increase the mechanical action range of the particles inside the model, thus increasing the shear strength of the model.
Secondly, when the x 1 and x 6 are fixed, there is also an increasing trend of the shear force with the increase in the x 2 and x 3 (Figure 5a,b,d-f).The reason is that in the shear process, the model will be not only subject to the action of the shear force but also to the action of the squeeze force, F a in the component force also has a certain squeeze force, so when the increase in x 2 and x 3 , the shear force will also increase.
Finally, the trend of the shear force with x 6 is shown in Figure 5c,e,f.It can be seen that when the x 2 , x 3 and x 1 are fixed, the shear force increases with the increase in x 6 .This is because a larger x 6 between particles will increase the number of model bonding bonds and change their spatial distribution, which will enhance the mechanical properties of the whole model and increase the stability of the model, and the shear strength of the corresponding model will increase [26].

Optimal Parameters and Validation of Calibration Results
Based on the F a of the blade on the EUO samara, the simulation parameters were optimally solved, and the quadratic regression equation obtained by fitting was parametrically solved with the aid of Design-Expert software, setting the optimization target to 43.2 N. The best combination of the significance parameters was obtained as follows: the particle contact radius x 1 was 0.79 mm, the normal stiffness per unit area x 2 was 3.20 × 10 7 N•m −3 , the shear stiffness per unit area x 3 was 3.83 × 10 7 N•m −3 , the bonded disk radius x 6 was 0.83 mm, and other factors were fixed values in the single-factor experiments.The optimal combination of parameters was verified in a simulation experiment, and the maximum shear force was 43.6 N. The relative error e between the simulation and the measured maximum shear force was 0.93%; the simulation results are in general agreement with the measured results, which shows that the established discrete element model of EUO samara and the calibrated parameters are accurate and reliable, and can guide the design and optimization of EUO samaras shelling equipment through discrete element simulation.

Conclusions
In this paper, based on the Hertz-Mindlin with bonding model, physical experiments and virtual calibration experiments were combined to calibrate the simulation parameters of the discrete element model of EUO samara, and the best combination of parameters for the EUO samara bonding model was obtained and validated, with the following main conclusions: (1) The moisture content and the density of the EUO samara were measured at 7.51% and 0.67 g/cm 3 , respectively.The bench experiments and physical experiments measured that the Poisson's ratio of the juniper fins was 0.403, and the modulus of elasticity was 19.98 MPa.The coefficient of restitution of EUO samara-to-samara and EUO samarato-stainless steel were 0.373 and 0.335, the coefficient of static friction were 0.676 and 0.442, respectively, and the coefficients of rolling friction were 0.349 and 0.254.Moreover, the maximum shear force of the blades on the EUO samara was 43.2 N. (2) The range of simulation parameters for EUO samara was initially determined with reference to similar materials.The influence trend of each factor on the shear force was analyzed by single-factor experiments, and the parameter range was reduced.The significant factors obtained from the two-level factorial experiments were particle contact radius, the normal stiffness per unit area, the shear stiffness per unit area, and the bonded disk radius.(3) The steepest climb experiments were used to narrow the range of significance parameters further.A quadratic regression model was established between the significance factors and the maximum shear force, while the optimal combination of the relevant parameters was obtained by optimal solution, and the following results were obtained: the particle contact radius was 0.79 mm, the normal stiffness per unit area was 3.20 × 10 7 N•m −3 , the shear stiffness per unit area was 3.83 × 10 7 N•m −3 , and the bonded disk radius was 0.83 mm.(4) The simulation was validated by using the best combination of parameters.The relative error between the simulated and measured values was 0.93%, which showed that the established discrete element model of EUO samara and the calibrated parameters were accurate and reliable.All of these could guide the design and optimization of EUO samara hulling equipment through discrete element simulation.

Figure 5 .
Figure 5. Response surface diagrams: (a) the particle contact radius and the normal stiffness per unit area, (b) the particle contact radius and the shear stiffness per unit area, (c) the particle contact radius and the bonded disk radius, (d) the normal stiffness per unit area and the shear stiffness per unit area, (e) the normal stiffness per unit area and the bonded disk radius, and (f) the shear stiffness per unit area and the bonded disk radius.

Figure 5 .
Figure 5. Response surface diagrams: (a) the particle contact radius and the normal stiffness per unit area, (b) the particle contact radius and the shear stiffness per unit area, (c) the particle contact radius and the bonded disk radius, (d) the normal stiffness per unit area and the shear stiffness per unit area, (e) the normal stiffness per unit area and the bonded disk radius, and (f) the shear stiffness per unit area and the bonded disk radius.

Table 2 .
The Single-factor experiment level.

Table 4 .
Analysis of variance of the two-level factorial experiments.

Table 5 .
The steepest climb experiments protocols and results.

Table 6 .
The Box-Behnken response surface experiments protocols and results.

Table 7 .
Analysis of variance of the Box-Behnken response surface experiments.