Calibration of discrete element parameters for simulating wheat crushing

Abstract The mechanical properties of wheat grains were measured and analyzed using discrete element software, which provided crucial data for their processing in a mill. The foundational Hertz–Mindlin model was used as a theoretical framework to evaluate the accumulation angle of wheat grains. The Plackett–Burman, steepest ascent, and Box–Behnken methods were utilized in a series of experiments, with the accumulation angle serving as the dependent variable. Targeting the actual angle of repose in the response surface, the optimal parameters were derived using the regression equations. These included a static‐friction coefficient of 0.3 between individual wheat grains, rolling‐friction coefficient of 0.04 for wheat–wheat interactions, static‐friction coefficient of 0.554 for wheat–tooth roller interactions, collision recovery coefficient of 0.5 for wheat–wheat collisions, collision recovery coefficient of 0.45 for wheat–tooth roller collisions, and rolling‐friction coefficient of 0.05 for wheat–roller interactions. Relying on the bonding contact model of Hertz–Mindlin, virtual uniaxial compression tests were conducted to calibrate the wheat grain bonding parameters. A regression equation for the critical load was subsequently generated using the critical load of the wheat grain bonding model as the response variable. The optimal parameters were calculated and incorporated into the EDEM model for computation, which resulted in a relative error of 1.6% between the calculated and observed values, confirming the accuracy and feasibility of the calibration method, suggesting that the calibrated parameters were accurate.

the milling mechanism within a mill.The EDEM software allows for the visualization of the entire wheat-crushing process.The discrete element method (DEM) introduced by Cundall in 1979 (Cundall & Strack, 1979;Reeves, 2013;Zhong et al., 2016) has been widely applied in diverse fields, including the agriculture, food, chemical, and pharmaceutical industries.When using simulation to study the wheat-milling process, the physical and contact parameters of wheat are critical to understanding its mechanical properties.
Consequently, researching these parameters prior to analyzing the milling mechanism is of utmost importance.Currently, these parameters are acquired primarily through direct measurement and virtual calibration.Wang et al. (2012) measured the recovery coefficients of wheat and rapeseed at various moisture levels based on bounce tests.Liu et al. (2016) utilized a Vernier caliper to measure the long and short axes of wheat grains, treating them as regular ellipsoids.The parameters of the wheat's discrete element model were calibrated using five-ball combinations.However, the wheat's bonding parameters were not examined.Zhang et al. (2020) employed three-dimensional scanning and reverse engineering methods to capture the contours of grains of rice and created balls with different radii for discrete element models to calibrate the staticfriction and dynamic-friction coefficients between grains.However, the accuracy of three-dimensional laser scanning for smaller objects is limited, making it incapable of fully capturing the threedimensional contour of the object.Ying et al. (2018) developed a discrete element model of a wheat grain based on X-ray technology, and calibrated the discrete element parameters of wheat.However, several studies (Chung & Ooi, 2006;Favier et al., 1999;Zeebroeck et al., 2003) have indicated that increasing the particle accuracy significantly extends the simulation time.Furthermore, numerous studies (Härtl, 2008;Markauskas et al., 2010;Pasha et al., 2016) have suggested that an accurate particle shape model does not significantly enhance the simulation accuracy.A rough, non-spherical model constructed with a limited number of spherical units usually offers better model accuracy.
In this study, building upon global research and using wheat particles as the subject, a discrete element model of wheat grains was developed.A cylindrical lifting experiment was conducted to determine the angle of repose of wheat particles.The Plackett-Burman, steepest ascent, and Box-Behnken experimental designs were employed to calibrate the parameters of the discrete element simulation.A comparison between the simulation values and actual experimental results was made, leading to the acquisition of parameters for the discrete element calibration of wheat particles through simulation.Incorporating Hertz-Mindlin's bond theory, and in conjunction with a compressive failure test on wheat particles, a response surface analysis was used to determine bond parameters, such as the normal contact stiffness, tangential contact stiffness, critical normal stress, and critical tangential stress.This will contribute to future research on the mechanism of wheat milling.

| Experiment model
This study considered three varieties of wheat: Zhengmai379, Xin-mai26, and Zhoumai36.These varieties are among the most common and representative in China.The moisture content of the wheat grains was determined to be 15.87% (wet basis reference) using a halogen moisture analyzer, which met the necessary condition for wheat processing.A total of 100 wheat particles were randomly and independently selected (Xu et al., 2023), and their dimensions (length, width, and thickness) were measured using a caliper (with a resolution of 0.01 mm).The average dimensions were determined to be 6.15, 3.06, and 2.9 mm for the length, width, and thickness, respectively.The measured results are presented in Figure 1 and Table 1.
To simplify the calculations and reduce the computation time, the equivalent diameter of non-spherical particles was employed, despite the wheat grains being triaxially non-uniform particles.The volumetric diameter was calculated using the formula outlined by Xu et al. (2023).
As observed from the three-axis size measurements of the wheat kernels listed in Table 1, the mean value for the length of the wheat kernels was 6.15 mm, which was used for the length of the ellipsoid.
The mean values for the widths and thicknesses of the wheat kernels were replaced by a circle with a diameter of 3.00 mm (radius R = 1.50 mm).Wheat particles were modeled in the Solidworks software, and the established wheat model was imported into the EDEM software.Based on this, the true shape of wheat grains was approximated, and a simulation model of the wheat grain was constructed using a five-ball combination, as shown in Figure 2. (1)

Practical application
This study primarily focused on the calibration of parameters for simulating wheat grain processing in discrete element software.The discrete element method is a crucial approach in simulating wheat processing, wherein the precision of the wheat's mechanical parameters significantly influences the outcome of the simulations.This study mainly investigated the mechanical property parameters of wheat, laying the foundation for studying wheat milling in flour production.Thus, the research reported in this paper holds substantial importance.
In this study, the cylinder lifting method was employed to measure the angle of repose of wheat.Considering the size of the wheat particles, a cylinder with an inner diameter of 27 mm and a height of 190 mm was used.For the measurement, the cylinder was placed on the horizontal fixed plate of a grinding roller, filled with 30 g of wheat particles, and then slowly elevated at a speed of 0.01 m/s, allowing the particles to accumulate on the grinding roller's flat surface.This experiment was repeated 10 times to obtain the angle of repose of the wheat grains.The experimental setup is illustrated in Figure 3.
A computer image processing technology was used to quantify the angle of repose more accurately.The accumulated particle image was processed in Matlab software, with the image's boundary pixels acquired via grayscale conversion, threshold segmentation, and boundary searching.Finally, the accumulated angle of the particles was obtained by linear fitting.The angle of repose was 31.74°, as shown in Figure 4.

| Simulation parameters
The approximate ranges of the simulation parameters in this study were preliminarily determined based on the relevant parameters of the wheat particles and materials of the grinding roller in the discrete element simulation (Akkoyun & Arslan, 2020;Boac et al., 2009;Frączek et al., 2007;Jia et al., 2014), and are listed in Table 2.

| Simulation model of cylinder lifting
The simulation model agreed with the experimental model.A particle factory was established at the top of the cylinder to generate a particle mass of 30 g, allowing the particles to fall freely.After all the particles reached the specified position and remained stationary, the cylinder was raised upward at a speed of 0.01 m/s until all the wheat particles fell off the cylinder and formed a stable pile, as shown in Figure 5.
In this study, the built-in Hertz-Mindlin (no-slip) contact model in EDEM was utilized for computations.Given the significant volume of calculations, all the simulation parameters were chosen to have a Rayleigh timestep of 20%.Moreover, the size of the grid in the simulation was set to be three times the size of the smallest spherical element.

| Plackett-Burman test
Leveraging the angle of inclination formed by wheat seeds as the response value, the Plackett-Burman test was employed to filter out factors with significant influence on the test and eliminate factors with little impact on the experimental results.This approach reduced the number of experiments and ensured the accuracy of the experimental test results.Six factors affecting the test indices in the cylinder lifting test were selected for examination, and the Plackett-Burman test factor table is presented in Table 3.According to the Plackett-Burman test principle, 12 groups of tests were designed, and the angle of repose for each test group was documented.The design and results of the simulated test are illustrated in Table 4.
The results of the Plackett-Burman test were analyzed to determine the dominance of each factor, as shown in Figure 6.This analysis provided insight into the impact of each factor on the angle of repose of the wheat particles.Among the factors, the static-friction coefficient for wheat-wheat collisions, rolling-friction coefficient for wheat-wheat interaction, and static-friction coefficient for wheat-grinding roller interaction were dominant, with a confidence interval of 95%.The influence of other factors on the angle of accumulation was insignificant.

| Steepest climbing test
To determine the optimal ranges for the dominant factors from the Plackett-Burman test and provide a foundation for subsequent testing, the steepest ascent test was necessary.The Plackett-Burman test identified three dominant factors.Six gradient points were considered for these dominant factors, while median values were used Size and sketch of a wheat particle.
TA B L E 1 Statistical analysis results of wheat triaxial dimensions.for the non-dominant factors.A total of six tests were designed, and the angle of repose was recorded for each.Error values were calculated according to Formula (2).The design and results of the steepest ascent test are presented in Table 5.
As listed in Table 5, as the static-friction coefficient of the wheat-wheat interaction, rolling-friction coefficient of the wheat-wheat interaction, and static-friction coefficient of the wheat-grinding roller interaction increased, the relative error initially decreased and then increased.The smallest relative error was observed in the third group test.Therefore, the optimal value ranges for the test factors could be selected around the third group test, that is, between the second group and fourth group tests.The parameters from the third group test were selected as the center points, with the parameters from the second and fourth group tests selected as high and low levels, respectively, for the Box-Behnken test.

| Box-Behnken test and regression model
The Box-Behnken test was designed using the results of the Plackett-Burman test and steepest ascent test, and 15 experiments were conducted using the Box-Behnken design principle.The angle of repose of each test was recorded.The Box-Behnken design and results are listed in Table 6.
(2) After analyzing the experimental results, a quadratic regression model was established to illustrate the influences of the parameters for the contact between the wheat and grinding roller on the angle of repose.Three substantial factors were also identified.This regression model is depicted in Formula (3).
The results of a variance analysis of the regression equation are presented in Table 7. Specifically, a value of p < .05for the fitted model signified that the regression equation was a good fit.The static-friction coefficient of the wheat-wheat interaction and rolling-friction coefficient of wheat-wheat interaction had values of less than 0.01, indicating that these two factors significantly affected the angle of repose of the wheat particles.Meanwhile, the quadratic terms of the static-friction coefficient of the wheatgrinding roller interaction and the static-friction coefficient of the wheat-grinding roller interaction had values of less than 0.05 and greater than 0.01, respectively.This indicated that these two factors significantly impacted the angle of repose of the wheat particles.

| Response surface analysis
The response surface was generated to analyze the effects of the static-friction coefficient of the wheat-wheat interaction, rollingfriction coefficient of the wheat-wheat interaction, and staticfriction coefficient of the wheat-grinding roller interaction on the angle of repose, and the results are shown in Figure 7.
The regression equation was used to calculate the optimal parameters with respect to the actual angle of repose: the optimal static-friction coefficient for the wheat-wheat interaction was 0.3, the optimal rolling-friction coefficient for the wheat-wheat interaction was 0.04, and the optimal static-friction coefficient for the wheat-grinding roller interaction was 0.554.These optimal parameters were then incorporated into the EDEM software to verify the simulation test.After conducting five verification experiments, the error rate between the simulated angle of repose and actual experimental angle of repose was found to be only 0.9%.This result corroborated the reliability of the parameter calibration, thus confirming its suitability for studying the wheat-milling process and defining the simulation parameters.
When the force between the small particles exceeded the maximum normal contact stress or maximum tangential contact stress of the bonding contact model, the bond broke.This separation of small particles signaled the fracture of the wheat particles.The calculation for this scenario is expressed in Formula (4).
The maximum normal and tangential stress values are defined in Formula (5), where A denotes the contact area, A = R B

2
; R B denotes the bonding radius; J denotes the moment of inertia of the particle, J = 0.5 R B 4 ; S n and S t denote the normal and tangential stiffness values of particles, respectively; δ t denotes the step size; v n and v t denote the normal and tangential velocities of particles, respectively; and w n and w t denote the normal and tangential angular velocities of particles, respectively.
According to the principles underlying the bond model, the rupture of the bonds between particles was related to parameters, such as the normal contact stiffness, tangential contact stiffness, critical normal stress, critical tangential stress, and bond radius of the wheat particles.To simulate the compressive failure of the wheat more accurately, we generated 512 particles by adding small particles with a diameter of 0.2 mm to the original wheat model.
The upper compression plate was then pressed downward at a rate of 6.6 mm/min.The generated wheat particles and compressed wheat are shown in Figure 8. Specifically, Figure 8d represents the TMS-PRO system, which can test the mechanical characteristics of grain, obtaining the destructive force, destructive energy, force-displacement curves, elastic modulus, and more.The forcedisplacement curve of the wheat grain compression process is shown in Figure 9.As shown in Figure 9, the wheat grains exhibited yield points during the compression process.Prior to reaching the (4) T n = −w n S n J t ; TA B L E 2 Parameters required in DEM simulation.

Simulation parameters Value
Poisson's ratio of wheat Density of grinder roll/(kg m −3 ) 7860 Wheat-wheat restitution coefficient 0.1-0.9 Wheat-acrylic static friction coefficient 0.1-0.6 Wheat-wheat rolling friction coefficient 0-0.1 Wheat-grinder roll restitution coefficient 0.1-0.8 Wheat-grinder roll static friction coefficient 0.3-0.8 Wheat  yield point, the pressure incrementally increased with the compression displacement, and dropped sharply after the yield point was reached.

| Results and analysis of bonding model
To determine the optimal bonding parameters of wheat particles, we designed a Box-Behnken experiment to obtain a regression equation that captured the influences of these factors on the relative error rate.Building upon related research (Huang et al., 2014;Oladunmoye et al., 2014;Zeng et al., 2018), we initially selected the bonding parameters listed in Table 8 for the wheat.
To determine the critical load of the bonding model for wheat grains, we used it as the response value and designed 27 sets of experiments.Table 9 illustrates the experimental design and respective results.An analysis of these results produced a second-order regression model that accounted for the impact of the four parameters of the wheat grain bonding model, as illustrated in Equation ( 6).
Table 10 presents the results of a variance analysis of the regression equation.The coefficient of determination, R2, of the re-

| Model verification
With the actual critical load as the target for the response surface, the regression equation was used to calculate the optimal parameters.The optimal normal stiffness was 3.61 × 10 9 N/m 3 ; optimal tangential stiffness was 1 × 10 7 N/m 3 ; optimal critical normal stress was 1 × 10 7 Pa; and optimal critical tangential stress was 1 × 10 8 Pa.We input these optimal parameters into the EDEM software to validate the simulation test.After conducting five tests, we found that the error rate between the simulated critical load and regression equation was 1.6%.This low error rate confirmed the reliability of the parameter calibration and indicated that these parameters could be used for simulating the wheat-milling process.

| CON CLUS IONS
Considering the challenges observed during the wheat-milling process, this study thoroughly investigated the recovery coefficient, as well as the static-and rolling-friction coefficients for contacts between wheat particles and between the wheat and grinding roller.
Additionally, we explored the influence of various factors on the wheat's angle of repose.A failure model was established to simulate the compression-induced fracturing of wheat particles, and we also examined the impact of bond parameters on the destructive force exerted on wheat.The key findings from this research are as follows.
1. Significant contact parameters influencing the angle of repose of wheat were determined.The optimal static-friction coefficient for the wheat-wheat interaction, rolling-friction coefficient for the wheat-wheat interaction, and static-friction coefficient for the wheat-grinding roller interaction were determined to be 0.3, 0.04, and 0.554, respectively.These optimal parameters were then applied in EDEM for simulation tests.After conducting five experiments, it was verified that the discrepancy between the discrete element simulation angle of repose and the actual test angle of repose was only 0.9%.This validated the reliability of the parameter calibration, which could be used to examine the wheat-milling process and determine the simulation parameters.

E 2
Model of wheat seed for discrete element.(a) Spherical position map of wheat particles.(b) Spherical filling diagram of wheat particles.

F
bonding model was introduced via the API function of the discrete element software.The bonding model could bond the small particles and prevent normal and tangential movement between (3) = −9.2+ 12.4x 2 + 251x 3 + 124.5x 5 + 35.2 x 3 + 38.7x 2 x 5 − 55x 3 x 5 Fitting image of accumulated angle.(a) Original plot of accumulation angle.(b) Bipartite graph of accumulation of angle.(c) Fitting image.
Simulation process diagram.(a) Complete filling of wheat particles.(b) Accumulated angle after cylinder lifting.
The Plackett-Burman test of the factor table.
gression equation was close to one, indicating a high degree of fit for the regression equation.Specifically, a value of p > .05for the fitted model suggested that the regression equation was well fitted.The p values for the normal stiffness, tangential stiffness, quadratic term of the normal stiffness, and interaction term between the normal stiffness and tangential stiffness were <.0001, indicating their significant impacts on the critical load of the wheat grain bonding model.However, the critical normal stress and critical tangential stress, which had p values >.05, appeared to have no substantial influence on the critical load.The results of the response surface analysis illustrated the effects of the normal contact stiffness, tangential contact stiffness, normal critical stress, and tangential critical stress on the crushing force of wheat, and the results are shown in Figure 10.

F
Experiment and simulation of wheat extrusion.(a) Generated wheat particle model.(b) Wheat grain bond map.(c) Simulation calculation of wheat compression.(d) Compression experiment of wheat.F I G U R E 9 Force-displacement curve of compression of wheat grain.TA B L E 8 Test factors.

F
Effects of various factors on the critical load.(a) Effects of y 1 and y 2 on the critical load.(b) Effects of y 1 and y 3 on the critical load.(c) Effects of y 1 and y 4 on the critical load.(d) Effects of y 2 and y 3 on the critical load.(e) Effects of y 2 and y 4 on the critical load.(f) Effects of y 3 and y 4 on the critical load.
The Plackett-Burman trial design and results.
F I G U R E 6 Significance of Plackett-Burman test parameter.
Box-Behnken trial design and results.Regression equation analysis of variance.

Source of variation Degrees of freedom Sum of squares Mean square p-value probability > F
Experimental design and experimental results.
TA B L E 9