Quantitative assessment and optimization of parallel contact model for flexible paddy straw: a definitive screening and central composite design approach using discrete element method

To simulate the bending behaviour of paddy straw at varied moisture contents after crop harvesting, we created a flexible paddy straw specimen model based on the Hertz–Mindlin with parallel contact bonding model using the discrete element model (DEM) approach. The research presented in this study aims to investigate a new approach called Definitive Screening Design (DSD) for parameterizing and screening the most significant parameters of the DEM model. This investigation will specifically focus on the three-point bending test as a means of parameterization, and the shear plate test will be used for validation purposes. In addition, the most influential DEM parameters were optimized using another Design of Experiments approach called Central Composite Design. The findings from the DSD indicated that parameters such as bonded disk scale, normal stiffness, and shear stiffness have the highest impact on the bending force, while the coefficient of static friction (Straw-Steel) has the least effect. The three bonding parameters were respectively calibrated with the loading rate (0.42, 0.5, and 0.58 mm s−1) and a good agreement between actual and simulated shear force at moisture content M1—35 ± 3.4%, M2—24 ± 2.2% and M3—17 ± 2.6%. Modelled stem helps simulate the straw with low error and increases the accuracy of the simulation. The validated model, with an average relative error of 5.43, 7.63, and 8.86 per cent, produced reasonable agreement between measured and simulated shear force value and loading rate.


Discrete element modelling (DEM)
Recent advancements in computer technology have facilitated the utilization of DEM for simulating granular materials in industrial and laboratory settings 21 .This methodology entails several sequential steps condensed into a single time step, typically involving the following four stages: (1) Initial preparation, which includes calculating particle field attributes, size, and arrangement; (2) Defining the physical characteristics of particles (e.g., mass, motion behaviour, particle interactions, and failure criteria) and applying boundary conditions; (3) Computations for determining particle forces, velocities, accelerations, and positions; (4) Post-processing to generate simulated particle images.These stages collectively constitute the DEM process.

Modeling of the flexible paddy straw
The paddy straw was developed by adding the multi-spheres in meta particles.With a parallel contact model, Hertz-Mindlin was used to develop flexible paddy straws.Bonding properties were selected such that the developed model can represent the paddy straw's actual behaviour.The sphere's diameter defines the straw diameter selected according to the average diameter of paddy straw in the physical experiment.

Flexible straw model
A solid sphere (diameter: 5 mm, straw length: 80 mm) was obtained by considering the shape of the paddy straw, accounting for the overlap between two spheres and having a given contact radius.Uniformly filled the straw unit with 60 spheres using Altair EDEM 2021 software (DEM Solutions Ltd., Edinburgh, UK).
The Hertz-Mindlin (no-slip) contact model was used in the straw model as default in the EDEM programme.Due to its inherent benefits in terms of accuracy and effectiveness, as shown in Fig. 1, this model is preferred.The tangential force component of this particular model was based on research 22 , whereas the normal force component was calculated using the Hertzian contact theory.It is noteworthy that damping elements are present in the normal and tangential forces, with the related recovery and damping coefficients.The total normal force is made up of both the normal force component and the damping normal force component in the context of the Hertz-Mindlin (no-slip) contact model.
where K n , C n and v rel n are the spring constant, damping coefficient and relative velocity in the normal direction, respectively.To model flexible paddy straw within the DEM framework, several sphere units were organized and (1) www.nature.com/scientificreports/interconnected using the Hertz-Mindlin with parallel contact V2 bonding model and these forces and moments can be computed using Eqs.( 2)-( 4) 23 . ).
In addition, the flexible straw would be broken when the maximum normal or tangential bonding stresses (σ max , τ max ) exceeded their critical values (σ critical , τ critical ).

Definitive screening design (DSD)
As suggested by Favre and Chaves in 2021, the DSD is a revolutionary statistical screening technique that makes it possible to examine the impact of four or more process factors on one or more response variables while needing fewer runs than conventional screening designs.Consequently, the input factors employed in DEM software simulations underwent a parametric evaluation to explore their impact on responses and identify the most influential parameters for subsequent calibration.This study utilized a Design of Screening Experiments (DSD) approach to conduct parametric assessments on twelve distinct input parameters.These parameters encompassed various characteristics such as the shear modulus of straw, restitution coefficient between straw and straw, static friction coefficient between straw and straw, rolling friction coefficient between straw and straw, restitution coefficient between straw and steel, static friction coefficient between straw and steel, rolling friction coefficient between straw and steel, normal stiffness, shear stiffness, normal strength, shear strength, and bonded disk scale.
Table 1 outlines the operational ranges of these parameters, which are selected based on a comprehensive literature review and prior research experiences.The straw's Poisson's ratio was predetermined at 0.3, as suggested by 24 , and was therefore considered a constant, not subject to variation in the DSD.However, the minimum run size for DSD rose by four runs as a result of the addition of two extra variables, leading to an expanded DSD with an actual total of 29 runs (Table 2) 25 .

DEM model of three-point bending test
To better understand the flexibility of straw, the three-point bending test DEM model was simulated using the Altair EDEM 2021 software, and a setup for the test was included in the straw model (Fig. 2a).Using Kelvin-Voigt models, these two spheres were connected by various linear, parallel spring-damper systems (Fig. 2b).The (2) www.nature.com/scientificreports/bending stiffness of the joints is examined through the stiffness (kb) and placement of the springs.At the joints, dampers (cb) prevent oscillations and guarantee simulation stability, while the straw is subjected to a normal force at the support point.To account for this, an extra set of Kelvin-Voigt systems (kr and cr), specifically in the radial direction (component 2 in Fig. 2b), is included per segment and joint.An additional spring (kt) and damper (ct), both of which are coupled to the two discs within a specific segment (Fig. 2a), are responsible for the tensile stiffness of each segment.As shown in Fig. 2b, the parallel contact version 2 and hysteretic spring contact model are used in the 3-point bending simulation.The flexibility and plasticity properties of the material are considered in this model.
A static block factory generated and fixed multi-sphere particles with a diameter of 5 mm over a three-point bending support of 80 × 5 × 5 mm.The particles followed a normal distribution with a mean and standard deviation of 1 and 0.1 (Fig. 3a).A particle's velocity and kinetic energy were allowed to decrease until they were less than 0.05 m/s and 0.01 J/s, respectively, before being allowed to settle over the two supports.The compressor geometry was given a downward linear translation velocity with a variable downward speed for the probe.

Physical three-point bending test
The test was organized to govern the maximum force to bend the straw up to a certain deformation.The straw length was 80 mm and rested on steel supports (Fig. 3b).A texture analyzer (TA-XT Plus, Stable Micro Systems Ltd., Surrey, UK) with a 490 N load cell capacity was used to conduct the TPA.The test was conducted using the A/3PB probe.Three different levels of moisture content and load rates were used to examine the samples.Three loading rates were used in a bending test to assess the force at the failure region of the straw.

DEM calibration for straw using CCD
After identifying the most dominant factors through DSD in JMP software, which can handle multiple factors simultaneously, the optimization (calibration) phase was carried out using another DOE method, Central Composite Design (CCD) by Design Expert 13 Software.Actual bending force at three moisture contents was obtained during the physical test and kept as the targeted value for optimizing straw DEM properties.

Model validation
A Shear Plate Test (SPT) validated the established DEM model.The dimensions of the shear plate used for the experiment are represented in Fig. 4. Comparisons were made between the shear force detected during the physical test and the DEM estimates for the shear force detected through the simulation test.www.nature.com/scientificreports/Shear test apparatus A force-time deformation graph was developed during the experiment, and textural parameters such as shear strength (SS) (Eq.6) were found with the support of software and the apparatus mentioned in the above section.The shear plate setup consisted of three rectangular plates (The outer plate is fixed, and the center plate is movable) of aluminium material (90 × 75 × 10 mm) having holes of 2-7 mm diameter.The aluminium plates were spaced 5 mm apart 26 .The force was measured using TPA when attached to a Warner Bratzler (HDP/BSW) probe www.nature.com/scientificreports/ to push the sliding plate of the shear plate setup.The data was recorded in TPA, and the force-time deformation graph was plotted on the screen during the experiment.A similar shear plate test was used to determine the shear strength properties of straw 27 .

Shear plate test DEM simulation
The shear plate test's CAD model (Fig. 4c) was converted into software with exact replicas of the real experiment's dimensions.The shear plate experiment was then modelled using the EDEM 2021 program on a computer with an HP Intel Core i9-9900K CPU running at 3.6 GHz and 32 GB of RAM.A similar length of DEM straw model was generated in the block factory inside the hole of the shear plate by adding the spheres as Meta particles to 80 mm length of straw (Fig. 4d).The calibrated model parameters provided in the table were used to run the simulations.The particles were allowed to settle before the simulations started until their velocities and kinetic energies, respectively, fell below 0.05 m s −1 and 0.01 J. Using a four-core CPU solver, each simulation in the validation study was run once with a timestep of 2.19 × 10 −5 .
The operating conditions (i.e., downward velocity) are similar to the real experiment and kept in the simulated experiment.Outer two plates were fixed while the middle plate was given the linear translation motion towards the downside at certain deformation as experienced in the physical experiment.The correlation between the simulated and actual bending force was evaluated using relative error.The following formula was used to compute the relative error.
The following formula from Eq. ( 8) was used to compute the relative error.
where RE Shearforce is the relative error (%), SF O and SF s are the observed and simulated shear force (N).

Results
The effects of various input parameters on the response of the bending test within the DEM model were initially investigated parametrically using the Discrete Element Method (DEM).The parameters obtained from the DEM research were then refined using an additional Design of Experiments (DOE) called the Central Composite Design (CCD), and this calibration procedure was validated using the shear plate test.

Parameterization of DEM parameters using DSD
Using JMP 15 Pro, the data from the DEM simulation analysis was carried out to accomplish this goal; this strategy decomposes the response variable (Y) into two halves, indicated as Y ME and Y 2nd and shown in Eq. ( 9) 28 .
where Y ME is the expected value of Y determined by regressing Y on the main effects and fictitious variables, Y 2nd is produced when Y 2nd = Y − Y ME.The determination coefficient for the model, which compares predicted values to actual values of bending force, is 0.98, whereas the RMSE is 2.151 (Fig. 5).
The results presented in Table 3 demonstrate that the main effects, including CSF (Straw-Straw), CSF (straw steel), CRF (straw steel), normal stiffness, shear stiffness, and BDS, have a highly significant impact on the bending force (p < 0.01).Among these parameters, BDS has the greatest influence on the bending force, with a corresponding F value of 8.42.It was followed by normal stiffness (8.28), shear stiffness (1.99), CSF straw-steel www.nature.com/scientificreports/(1.67), CRF straw-steel (0.94), and CSF straw-straw (0.74).The bending force is positively impacted by all relevant factors, meaning that when their values rise, the bending force (BF) also tends to rise.The BDS of straw has a significant effect on the bending force (BF).Increasing the BDS value from 3 to 5 results in an increase in the bending force.Similarly, the Normal stiffness has a linear relationship with the bending force, where an increase in the Normal stiffness from 1 × 10 8 to 1 × 10 9 leads to an increase in the bending force (Fig. 6).
Regarding shear stiffness, initially increasing the shear stiffness from 1 × 10 8 to 7 × 10 8 increases the bending force.However, further increasing the shear stiffness from 7 × 10 8 to 1 × 10 9 decreases the bending force.The CSF of straw steel also significantly affects the bending force.Initially, as the CSF value increases from 0.3 to 0.55, the bending force also increases.However, beyond the value of 0.55, the bending force tends to remain relatively constant.

Calibration of DEM parameters of paddy straw model using CCD
Based on the findings outlined in Table 3, it can be seen that only four DEM factors have a significant influence on the BF.These factors are CSF (straw-steel), Normal stiffness, Shear stiffness, and BDS.The shear modulus of the straw was kept constant at 1 × 10 6 Pa because, beyond this limit, the straw may become unstable during the simulation 29 .
Table 4 gives information on four factors: their equivalent levels for CSF (straw-steel), Normal stiffness, Shear stiffness, and Bonded disc scale.There are three levels at which each of these parameters is defined.The Central Composite Design (CCD) matrix and the information gathered from the related trials are shown in Table 5.The model fit for the bending test using CCD and the data points line up nicely along a straight line, supporting the model's normality assumption.
The selected model's adjusted and predicted determination coefficient values were 0.99 and 0.98, respectively, and their difference was less than 0.20.This confirms the adequacy of the selected model to describe the effect of independent variables.The showing ANOVA for bending force using CCD.The R 2 value of 0.99 indicates good agreement between independent parameters and bending force.Table 6 indicates that normal stiffness, shear The interaction terms BC and BD had a significant effect on the bending force at 5% and 1% level of significance, respectively.In order to obtain the reduced polynomial model, insignificant terms were ignored from the quadratic model was represented in the Eq.(10).
The physical three-point bending test results showed that the maximum bending force was 2.6 N, 3 N, and 5.5 N at a loading rate of 0.42 mm s −1 applied on straw specimens representing moisture levels of M 1 , M 2, and M 3 , respectively (Table 7).The moisture content was mimicked in the DEM simulation by adjusting the bonding parameters while keeping the other parameters constant.
The maximum and minimum values of the bending force (BF) were found to be 6.8 and 1.15 N, respectively.The BF exhibited an increasing trend with an increase in shear stiffness up to 7.52 × 10 8 N m −3 , after which it  www.nature.com/scientificreports/decreased within the range of variables.A similar trend was observed for BDS.The BF linearly increased from 2.8 to 6.8 N with an increase in BDS (Fig. 7a) in the range of 3-5 at a shear stiffness of 1 × 10 9 N m −3 .Similarly, it increased from 4.4 N to 6.8 N at BDS 5 with an increase in shear stiffness ranging from 1 × 10 7 to 1 × 10 9 N m −3 (Fig. 7b).The BF increased with an increase in normal stiffness up to 6.8 × 10 8 N m −3 after which it decreased.The BF exhibited a relatively constant trend with an increase in BDS ranging from 3 to 5. The BF followed a concave-shaped curve, first increasing and then slightly decreasing within the normal stiffness range of 1 × 10 7 to 1 × 10 9 N m −3 at the shear stiffness range of 1 × 10 7 to 1 × 10 9 N m −3 (Fig. 7c).A similar trend was observed for shear stiffness while keeping the maximum value of Normal stiffness.The calibrated parameters of straw stem at 35 ± 3.4% moisture content and other selected values of the straw model were presented in Table 8.

Model validation
The method outlined in Section "DEM calibration for straw using CCD" was used to validate the constructed DEM straw model (Table 8).The simulated shear force at different loads, such as 0.42, 0.5, and 0.58 mm s −1 with the same deflection depth, is depicted in Fig. 8, and it shows good agreement with the physical shear test.A wellcalibrated model shows a similar trend (Fig. 8) for a maximum shear force at the abovementioned conditions with minimum relative error (Table 9).The maximum physical shear force was measured to be 101.8N, while the maximum simulated shear force was 104.4 N at moisture content M 1 and applied a loading rate of 0.42 mm s −1 .Similar values were observed for other combinations of moisture contents and loads, as depicted in Fig. 9.The highest simulated shear force value of 143 N was obtained at moisture content M3 at a loading rate of 0.58 mm s −1 .
The loading rates of 0.42 mm s −1 , 0.50 mm s −1 , and 0.58 mm s −1 resulted in maximum shear forces with average relative errors of 5.43, 7.63, and 8.86%, respectively (Table 9).The graph (Fig. 9) shows the simulated versus actual shear force at the loading rate of 0.42, 0.5, and 0.58 mm s −1 with an excellent correlation R 2 value of 0.99, 0.98, and 0.96, respectively (Fig. 10).

Discussion
Sensitivity analysis for the DEM straw model identified key parameters, given the challenge of handling numerous parameters in a single experiment 25 .Therefore, parameters with non-significant influence should be screened out first before proceeding with the actual experiment 34 .Parameterization of the DEM model revealed that bonded disc scale (BDS), normal stiffness and shear stiffness were the most influential parameters.BDS and Normal stiffness affect the bending force, whereas the bending force initially increases with increasing the shear stiffness (1 × 10 8 to 7 × 10 8 N m −3 ) and then decreases with further increasing the value of shear stiffness (7 × 10 8 to 1 × 10 9 N m −3 ).BDS affects linearly on the bending force; when the bonded disc scale is increased, the contact area between the straw particles and the other particles in the DEM is increased, resulting in an increase in bending force 7,9 .Higher normal stiffness means that it will be more difficult to compress the material, leading to a higher bending force.This is because when the disk contacts the straw, the straw resists compression, which causes an increase in the bending force.Higher shear stiffness means that it will be more difficult to deform the material, leading to a higher bending force 16 .The optimized (calibrated) values of DEM parameters at three different moisture content.The straw model was validated by the shear plate test setup with three different loading rates.Maximum and minimum bending force and shear force were varied from 2.3 to 8.9 N and 101.8 to 143 N 35 .The relative error between the actual and simulated model is less than 10 per cent under the acceptable limit and represents the real behaviour of straw 12 .

Figure 1 .
Figure 1.Schematic diagram of Hertz Mindlin contact model.

Figure 2 .
Figure 2. A bendable stem of crop (a) Placement of meta particles-based straw on three-point support (b) Representation of spring and mass damper between two particles 15 .

Figure 3 .
Figure 3. Three-point bending test (a) Simulated test in DEM (b) Physical test on TPA platform.

Figure 4 .
Figure 4. Shear strength measurement (a) Physical shear plate setup (b) dimensions of shear plates, mm and (c) labelled CAD model of shear plate setup (d) simulated shear plate setup.

Figure 5 .
Figure 5. Actual versus predicted plot for bending force using DSD.

Figure 7 .
Figure 7. Response surface of bending force (BF) (a) BDS and shear stiffness (b) shear stiffness and normal stiffness (c) bonded disk scale (BDS) and normal stiffness.

Figure 8 .
Figure 8.(a) Simulation of shear plate test in EDEM software (b) comparison of simulated and actual shear force.
www.nature.com/scientificreports/bending test with a coefficient of determination near unity (R 2 = 0.9965).The flexible straw could withstand normal and tangential displacements, which enabled the simulation of crop straw's deformation characteristics.When using the best bonding parameters in conjunction with moisture content, M 1 , M 2, and M 3 were calibrated and validated to be (normal stiffness per unit area: 5.05 × 10 8 , 5.05 × 10 8 and 1 × 10 9 N m −3 , Shear stiffness per unit area: 5.05 × 10 8 , 9.33 × 10 8 and 5.05 × 10 8 N m −3 , BDS 3.34, 3.6 and 3.65).The relative error of the average of maximum actual and simulated bending force at different load rates and moisture content was 10.4, 4.25, and 7.2%; similarly, for shear force, it varies from 5.43, 7.63, and 8.86%, respectively.It depicts the accuracy of the regression model used to predict the bending and shear force.The flexible paddy straw described in this study can be implemented in DEM models to simulate the interactions between straw particles and equipment.Similarly, other crops can benefit by adopting modelling and parameter calibration techniques.

Figure 9 .
Figure 9. Relationship between an actual and simulated shear force of paddy straw stem.

Figure 10 .
Figure 10.Actual versus simulated values of shear force.

Table 1 .
Selected input parameters and their ranges for DSD.

Table 2 .
Experimental design matrix for DSD.

Table 3 .
Estimates of model coefficients with corresponding standard errors (SE) and p values for bending force (BF) using DSD.**Significant at 1% level of significance, *significant at 5% level of significance, NS nonsignificant., and BDS significantly affected bending force at 1% significance level.The created model was statistically significant (p < 0.001), as shown in Table6, indicating that lack of fit is insignificant.As a result, it can be concluded that the constructed model has excellent statistical quality.Except for CSF (straw-steel) with Normal stiffness, Shear stiffness, and BDS, all of the main and interaction effects favourably impacted the response.

Table 4 .
DEM input parameters for the straw calibration using the CCD experiment.

Table 5 .
Experimental design matrix for CCD.

Table 6 .
ANOVA for the bending force using CCD.**Significant at 1% level of significance, *significant at 5% level of significance, NS non-significant.

Table 7 .
Calibrated parameters of the straw specimen using CCD.

Table 9 .
Observed and simulated shear force at different loading rates and moisture content with corresponding relative error.