- 1College of Biological and Agricultural Engineering, Jilin University, Changchun, China
- 2Research Center, Hohhot Branch of Chinese Academy of Agricultural Mechanization Sciences, Hohhot, China
The forage crop Caragana korshinskii Kom. is of high quality, and the biomechanical properties of its plant system are of great significance for the development of harvesting equipment and the comprehensive utilisation of crop resources. However, the extant research on the biomechanical properties of Caragana korshinskii Kom. is inadequate to enhance and refine the theoretical techniques for mechanised harvesting. In this study, we established a discrete element model of CKS based on the Hertz-Mindlin bonding contact model. By combining physical experiments and numerical simulations, we calibrated and validated the intrinsic and contact parameters. The Plackett-Burman design test was employed to identify the significant factors influencing bending force, and the optimal parameter combination for these factors was determined through response surface analysis. When the shear stiffness per unit area was 3.56×109 Pa, the bonded disk scale was 0.93 mm, the normal stiffness per unit area was 9.68×109 Pa, the normal strength was 5.62×107 Pa, the shear strength was 4.27×107 Pa, the discrete element numerical simulation results for three-point bending, radial compression, axial tension, and shear fracture exhibited a maximum failure force error of 3.32%, 4.37%, 4.87% and 3.74% in comparison to the physical experiments. In the cutting experiments, a smaller radial angle between the tool edge and the stem resulted in less damage to the cutting section, which was beneficial for the smoothness of the stubble after harvesting and the subsequent growth of the stem. The discrepancy in cutting force between the physical and numerical simulations was 3.89%, and the F-x (force versus displacement) trend was consistent. The multi-angle experimental validation demonstrated that the discrete element model of CKS is an accurate representation of the real biomechanical properties of CKS. The findings offer valuable insights into the mechanisms underlying crop-machine interactions.
1 Introduction
Caragana korshinskii Kom., a member of the Leguminosae family, is a shrub that is widely distributed in Inner Mongolia, Ningxia, Shanxi, Shaanxi, and Gansu provinces of China, as well as in Mongolia (Zhang et al., 2010; Yang et al., 2013). It can be utilised as forage and an ecological protection species, possessing considerable economic and ecological value (Liang et al., 2019; Roy et al., 2021). In order to promote population renewal and sustainable growth, stem cutting operations must be performed every three to five years, in accordance with the plant’s growth characteristics and agronomic requirements (Wang et al., 2013; Gao et al., 2021). The morphology of the stubble section after cutting is a pivotal factor influencing the rejuvenation of new branches (Gao et al., 2021). At present, the biomechanical properties of CKS are not sufficiently researched to support the development of harvesting equipment. This results in the residual stems exhibiting a lack of smoothness following harvesting, which in turn affects their subsequent growth and causes a significant waste of biological resources.
To provide evidence for the design of harvesting machinery, the rejuvenation of plants, and the increase in crop yield based on the mechanical properties of plant stems, many researchers have conducted extensive studies on the biomechanical properties of stems in crops such as barley, corn, wolfberry, tea, hemp, and other crops. These studies have played a significant role in guiding agricultural production activities (Khan et al., 2010; Liu et al., 2012; Aydin and Arslan, 2018; Chen et al., 2019; Ma et al., 2020; Guo et al., 2021; Luo et al., 2022; Xia et al., 2022). However, conventional physical testing techniques for stem biomechanical properties are only capable of examining the macroscopic aspects and are unable to elucidate the failure mechanisms within the internal biological tissues of the stems at the microscopic level. This limitation impedes a comprehensive investigation into the dynamic interaction between the stems and the cutting devices.
The discrete element method has been employed extensively in the construction of plant stem models in light of the accelerated advancement of computers (Schramm and Tekeste, 2022). (Fu et al., 2023) employed a combination of physical measurements and computer simulation to ascertain the parameters of micro-crushing of corn stalk pulverisation, thereby providing technical support for the development of stalk micro-crushing equipment. (Leblicq et al., 2016a) employed the discrete element method (DEM) to numerically simulate the compression and bending processes (Leblicq et al., 2016b) of corn stem, thereby demonstrating the efficacy of DEM in characterising the mechanical properties of the stem. (Huan et al., 2022) investigated the cutting and flinging process of king grass stalks in mechanical harvesting with a discrete element method, elucidating the mechanism of interaction between king grass stalks and harvesting machinery. (Lei et al., 2022) calibrated the discrete element numerical simulation parameters of the tail stem and tail leaves of crushed sugarcane (STL) by combining physical experiments and numerical simulation with an optimised design, and verified the reliability of the calibrated parameters. (Hu et al., 2023) calibrated the parameters of the DEM model of ramie, verified the accuracy of the ramie DEM model, and provided theoretical support for the simulation of ramie stalk peeling. (Jiang et al., 2021) established a model of a cotton stalk using the discrete element method and identified the optimal parameter combinations, which facilitated the simulation of the shredding apparatus and the analysis of the material preparation. (Wang et al., 2020) established the DEM of citrus fruit stalks, which is capable of accurately simulating the bending and shearing actions of citrus fruit stalks. The results of this study are of considerable importance for the optimisation of the end-effector of harvesting robots for citrus fruit stalks. (Liu et al., 2022) established a DEM model of taro using a combination of physical tests and numerical modelling. The model was verified through clamping and pulling experiments, thereby providing a reference for optimising the taro harvester. The aforementioned studies demonstrate the efficacy of the DEM as a methodology for simulating crop stalks. By integrating physical experimentation with computational modelling, it is feasible to delve deeper into the internal micro-mechanical attributes and failure mechanisms of plant stems. In comparison to conventional experimental testing techniques, computer-based numerical simulation offers a number of advantages, including the ability to generate accurate and reliable test results, facilitate repetition, and verify the accuracy of traditional experiments. However, there is a paucity of literature on the DEM model of CKS, and the field of numerical simulation research in this area remains largely uncharted territory.
The objective of this study was to conduct a series of mechanical tests on CKS, including three-point bending, radial compression, axial tension, and shear fracture tests. A discrete element model of CKS was established based on DEM, calibrated and validated. The optimal parameter combination for the CKS model was determined, and the maximum failure force and F-x curve trends were compared between physical experiments and numerical simulations. Cutting experiments were designed to analyse the cutting force of physical and numerical simulation tests to verify the correctness of the discrete element model and to provide a correct and reliable CKS discrete element model for the DEM numerical simulation of the operation process of the Caragana korshinskii Kom. harvesting equipment.
2 Materials and methods
2.1 Experimental materials
2.1.1 Experimental samples
Caragana korshinskii Kom. was used as the experimental material. Samples were gathered utilising the five-point sampling method (Ren et al., 2021) from Fanjiayao Village, Helinger County, Hohhot, Inner Mongolia Autonomous Region, at the experimental field of the Hohhot branch of the Chinese Academy of Agricultural Mechanization Sciences (40°16′N, 111°44′E). Only specimens of Caragana korshinskii Kom. that exhibited no obvious signs of pests or diseases were selected. The stems were pruned to remove any superfluous branches and leaves, leaving only the requisite stems for the production of experimental samples. In consideration of regional characteristics, the cutting operation of CKS was primarily concentrated between 0 mm and 120 mm above the ground (Chen et al., 2020), which was intercepted and the stems were cut into 100 mm and 20 mm. The 100 mm segments were employed in three-point bending, axial tension, shear fracture, and cutting tests, whereas the 20 mm segments were utilised in compression tests, as illustrated in Figure 1A. To ensure the accuracy of the test data, all CKS samples were maintained at a moisture content of between 35.6% and 44.3% (wet basis). The diameter of the CKS utilised in the physical experiments was 6 ± 0.47 mm and the diameter of the discrete element model of the CKS established in this study was set to 6 mm. To guarantee the precision of the experiments, all physical tests were conducted within 48 hours.
Figure 1. Experimental sample and device. (A) Experiment sample preparation; (B) Cutting experiment device. 1. Movable frame 2. Fixture 3. Pulley 4. Circular cutter 5. Guide rails 6. Belt 7. Motor; (C) MTS Exceed E43; (D) Three-point bending experiment; (E) Radial compression experiment; (F) Axial tension experiment; (G) Shear fracture experiment.
The CKS were subjected to four distinct tests: three-point bending (Figure 1D), radial compression (Figure 1E), axial tension (Figure 1F), and shear fracture (Figure 1G). These tests were conducted using an MTS Exceed E43 universal testing machine, as illustrated in Figure 1C. During the course of the tests, the force, displacement, and force-displacement (F-x) curves of the samples were recorded using the computer that was connected to the testing machine. The testing machine is equipped with a force sensor with a maximum capacity of 10 kN, and the testing speed ranges from 0.001 to 500 mm/min. The four distinct test procedures were conducted by utilising different fixtures. The loading speeds for the three-point bending, radial compression, and shear fracture tests were set at 15 mm/min, while the loading speed for the axial tension test was set at 2 mm/min. Each test method was performed on ten samples.
2.1.2 DEM model
2.1.2.1 Contact model
According to previous studies, the biomechanical properties of CKS could be better reflected by choosing the Hertz-Mindlin (no-slip) and standard rolling friction models to study the interaction between CKS particles and steel (Liu et al., 2018), on the basis of these models, the Hertz-Mindlin bonding contact model was used to describe the interactions between CKS particles, with parallel bonds representing the interaction forces between two particles, as illustrated in Figure 2A. Fn and Fn’ are the normal forces of the bond, Fs and Fs’ are the tangential forces of the bond; Mn and Ms are the normal and tangential moments of the bond, respectively. where kn and ks are the normal and shear stiffness per unit area, Vn, Vt, ωn, ωt are the normal and tangential velocities, and the normal and tangential angular velocities, respectively; Lb is the length of the bond; Rb is the bonded disk scale; Rp is the radius of the particles; and Rc is the radius of the contact between the particles. Rc is 1.2 to 2 times of Rp (Zhang et al., 2019), this spacing can help to generate the parallel bond better. The force acting on the particle and the action force were determined using (Equation 1) (Jiménez-Herrera et al., 2018).
Figure 2. The discrete element model of CKS. (A) Hertz-Mindlin contact model with the bonding contact; (B) The size of physical CKS and the DEM Model of CKS.
Where,
Where I, IP, and Ab are the cross-sectional moment of inertia, polar moment of inertia, and cross-sectional area of the bond, respectively. The forces F and moments M (Equation 2) acting in the normal and shear directions vary due to the relative normal and shear displacements of the sphere as well as the normal and shear rotations of the sphere. Both Fs and Fs’ act on the bond and the bond breaks down when these forces exceed a critical value, as shown in (Equation 3)
Where σmax and τmax are the normal and tangential critical stresses, respectively.
2.1.2.2 The DEM model of experimental materials
In DEM numerical simulations, the multi-sphere method had been a commonly employed technique for establishing irregular particle models (Liu et al., 2018; Zeng and Chen, 2019; Liao et al., 2020; Fan et al., 2022). In the context of discrete element numerical simulations of non-spherical particles, the multi-sphere packing model offered a number of advantages, including a relatively fast computational speed and a straightforward approach to contact judgment. However, as the number of filling particle elements increased, the simulation time increased significantly, necessitating the use of a reasonable number of filling particle elements (Kruggel-Emden et al., 2008; Cabiscol et al., 2018). The modeling process of the CKS was as follows: (1) Established the geometric model of the CKS used in the experiments in Unigraphics NX; (2) Imported the model into ICEM CFD software for structured meshing to obtain the coordinates of the centroids of each mesh element; (3) Used the element particle model in EDEM to generate the DEM model of the experimental CKS by matching the mesh centroid coordinates to the particle center coordinates. The diameter and length of the DEM model were maintained in accordance with the specifications of the physical experiments. A total of 1,660 standard spherical particles with a length of 100 mm and 340 standard spherical particles with a length of 20 mm were generated, as illustrated in Figure 2B (The methodology employed to define the particle parameters of the DEM model was outlined in the Supplementary Material.).
2.2 Determination of intrinsic parameters of the DEM model
The intrinsic parameters of the DEM model include density, Poisson’s ratio, and shear modulus.. The density ρ was determined by standard methods, with a value of 1014 kg·m-3 (ASTM D854-10, 2010), and the Poisson’s ratio was set to 0.4 (Madison et al., 2010). The shear modulus G was calculated through a three-point bending test, as illustrated in Figure 1D. In the three-point bending test, the CKS was positioned on the test bench with a sample length of 100 mm and a diameter of 5.93 ± 0.26 mm. The distance between the two support points was 80 mm, and a cylindrical indenter with a diameter of 10 mm was applied at a rate of 15 mm/min. The bending force and displacement (FB-x) curve was recorded by the MTS Exceed E43. The elastic modulus E is calculated using (Equation 4), and the shear modulus G is then obtained using (Equation 5), where y is the deflection (displacement), FB is the bending force, and S is the spanning distance between the support points. (scale length), E is the modulus of elasticity, J is the moment of inertia of the CKS cross-section, dF/dy is the slope of the FB-x curve, and μ is the Poisson’s ratio of the CKS.
The three-point bending experiment was conducted ten times. The average maximum bending force FB was found to be 115.67 N, while the average modulus of elasticity E was calculated to be 6.46 × 109 Pa. The shear modulus G was determined to be 2.31 ×109 Pa.
2.3 Determination of other parameters of the DEM model
The study established the contact parameters of the material and the parallel bond parameters of the DEM model. The methods for determining the material’s contact parameters were relatively well established; detailed measurement methods could be found in the Supplementary Material and would not be discussed here. In the Hertz-Mindlin bonding contact model of the DEM, five parallel bond parameters were used to characterise the interaction forces between two CKS discrete element particles: normal stiffness per unit area, shear stiffness per unit area, normal strength, shear strength, and bonded disk scale. The ranges for these five parameters were based on the findings of previous research (Luo et al., 2019; Liao et al., 2020). In order to ascertain the requisite level values for the parallel bond parameters, preliminary experimental studies were conducted. Table 1 presents the values for the contact and parallel bond parameters.
2.4 Optimized design of discrete element model parameters
Firstly, numerical simulation tests were conducted in the EDEM software, utilising the three-point bending test force FB as an indicator. The Plackett-Burman design test (Zhang et al., 2023) was employed to identify the significant factors affecting the outcome. Secondly, a mathematical model of FB and the significant factors was established using a response surface analysis methodology, and the optimal model parameter combination was established by solving the mathematical model. Finally, numerical simulations of three-point bending, radial compression, axial tension, and shear fracture were conducted using the DEM model with the optimal parameters. The maximum simulated failure forces and the trends of the force-displacement curves were compared with the physical test results to verify the accuracy of the DEM model of CKS.
2.4.1 Plackett-Burman design experiment
In this study, the Plackett-Burman module in Design-Expert software was employed to evaluate the effect of eleven factors on the bending force FB as the response value. Each factor was assigned both a high and a low level. A total of 15 experimental runs were conducted, with three centre points established at the outset. The experimental design and results are presented in Table 2. Additionally, a Pareto chart was established.
2.4.2 Response surface experiment
Based on the results of the Plackett-Burman design experiment, the Box-Behnken module in Design-Expert software was employed to undertake response surface optimisation experiments. The optimal parameter combination for the three most significant factors was identified. A total of 15 experimental runs were conducted, with three centre points established. The experimental design and results are presented in Table 3. Subsequently, the central composite design experiment was employed to ascertain the optimal parameter combination for the remaining two factors. Three centre points were established during the course of the experiments, resulting in a total of nine experimental runs. The experimental design and results are presented in Table 4.
2.5 Cutting experiments
A cutting test bench was employed for both physical and numerical simulation experiments, as illustrated in Figure 1B. The power generated by the motor was conveyed through a belt to operate the circular cutter. The CKS was affixed to a movable frame via a fixture, with a traction motor positioned on the traction side of the movable frame. This configuration enabled the CKS to be transported horizontally along the guide rails to the circular cutter. A high-precision sensor connected the traction motor and the movable clamping device transmited the force data to a computer for analysis. During the test, the circular cutter had a diameter of 400 mm, a rotational speed of 600 rpm, and a material feed rate of 0.6 m/s. The numerical simulation experiments were conducted under identical operating conditions to those of the physical experiments.
3 Results and discussion
3.1 Plackett-Burman design experiment analysis
As illustrated in Figure 3A, when the contribution of the DEM parameters exceeds the t-value limit, a notable impact is observed on the response. Conversely, when the contribution of the parameters is below the t-value limit, the effect on the response is almost negligible (Xia et al., 2019). A positive effect indicates that as the level of the factor in question increases, the value of FB also increases. Conversely, a negative effect represents the opposite phenomenon. The effects on FB, in descending order, are as follows: ks, Rb, kn, σmax, τmax, x3, x2, x1. The parameters that have a significant effect on FB are ks, Rb, and kn. Furthermore, the effects of the three significant factors are positively correlated.
Figure 3. Pareto chart and the response surface analysis. (A) Pareto chart of significant parameters; (B) Effects of ks and kn on FB; (C) Effects of Rb and kn on FB; (D) Effects of Rb and ks on FB; (E) Effects of σmax and τmax.
3.2 Response surface analysis
The Design-Expert software was employed to align the test outcomes of the Box-Behnken design with the fitting operation, thereby establishing a mathematical model of the bending force FB and the influencing factors, as illustrated in (Equation 6). Subsequently, the analysis of variance (ANOVA) was conducted, and the test protocols (factors not enumerated were assumed to be at the level of 0) and the outcomes are presented in Table 3.
Table 5 reveals that the model P-value is 0.0013, which is less than 0.01, indicating that the model is statistically significant. The P-value of the misfit term is 0.4805, which is greater than 0.05, suggesting that the misfit term is not significant. The coefficient of determination R2 is 0.9775, and the coefficient of variation is 6.31%. The adjusted coefficient of determination R2 adj is 0.9370, indicating that the model has a high level of predictive power. The difference between the adjusted coefficient of determination and the prediction coefficient R2 pre is less than 0.2, which suggests that the model does not have a large block effect and is therefore highly predictive. It can be observed that both ks and Rb have a highly significant effect on FB (P < 0.01), while kn has a significant effect on FB (P < 0.05). The remaining terms are not significant.
The response surface of the three significant factors affecting the FB (Figure 3B) demonstrates that the FB required for the numerical simulation tests of three-point bending increases in proportion to the values of kn and ks, in accordance with the theoretical basis presented in (Equation 1). In conjunction with Table 5, it is evident that ks exerts a more pronounced influence on FB than kn. This is due to the fact that when the model is subjected to bending forces, a pair of reversed tangential forces are generated internally, and an increase in ks enhances the bond’s capacity to resist deformation in the tangential direction (Wang et al., 2020). This, in turn, leads to enhanced tangential elastic recovery and stability, which improves inter-particle contact stresses and increases FB in order to break the bond. Figure 3C and Figure 3D demonstrate that an increased Rb also results in a larger required FB, which is consistent with the theoretical basis of (Equations 1, 2). As indicated by the P-values of the three significant factors in Table 5, the order of influence on FB is ks> Rb> kn, which aligns with the order of influence of the three factors as depicted in the Pareto chart (Figure 3A). The parallel bond parameter combinations were obtained by solving the response surface analysis of Design-Expert software. The results yielded the following values: ks is 3.56×109 Pa, Rb is 0.93 mm, kn is 9.68×109 Pa.
It should be noted that the remaining two parameters may fall outside the preset range. Therefore, in order to determine the values of σmax and τmax, the central composite design method was employed, utilising the optimal parameter combinations of ks, Rb, and kn determined in the previous section. The mathematical model of the bending force FB with σmax and τmax is presented in Equation 7.
As the values of σmax and τmax must be treated as positive numbers, the non-compliant test items were excluded from further consideration. Subsequently, ANOVA was conducted, and the resulting test outcomes were presented in Table 6. The model P-value was 0.0229, which was less than 0.01, indicating that the model was significant. The P-value of the misfit term was 0.9964, which was greater than 0.05, demonstrating that the misfit term was not significant. The coefficient of determination R2 was 0.9635, the coefficient of variation was 4.77%, and the adjustment coefficient R2 adj was 0.9028. The difference between the prediction coefficient R2 pre of 0.9179 and the value of 0.9028 was less than 0.2, indicating that the model did not have a large block effect and that it performed better in prediction. The values of σmax and τmax had a significant effect on FB (P < 0.05), while the remaining terms where not significant. The order of influence on FB was σmax > τmax, which was consistent with the order of influence of the two factors corresponding to the Pareto chart (Figure 3A).
Increasing σmax, τmax can raise the thresholds of normal strength and shear strength of the bond, without altering the remaining parameters in (Equation 3). The sole means of achieving this is to augment the FB, and the bond will be rendered unstable when the FB surpasses the normal strength and shear strength of the bond. This is in accordance with the observed trend of increasing σmax, τmax, which subsequently results in an increase of FB, as illustrated in Figure 3E. The five parallel bond parameters of CKS were determined by response surface analysis of Design-Expert software and are as follows: ks is 3.56×109 Pa, Rb is 0.93 mm, kn is 9.68×109 Pa, σmax is 5.62×107 Pa, and τmax is 4.27×107 Pa.
3.3 Validation of the correctness of the DEM model
3.3.1 Analytical verification of shear failure experiments
Figure 4A shows the shear failure physical and numerical simulation experiments of CKS. As illustrated in Figure 4B, the shear fracture process can be delineated into four distinct stages: elastic deformation, linear elastic deformation, shear failure, and fracture. The elastic deformation stage (OA) demonstrates a negligible variation in shear force with increasing displacement. This phenomenon can be attributed to the resilience of the phloem in the stem, where the load applied by the tool primarily compresses the stem without the blade cutting into it, resulting in elastic deformation of the stem (Guo et al., 2021). As the shear force increases, the process progresses to the linear elastic deformation stage (AB). In this phase, the shear force increases exponentially with displacement as the blade penetrates the stem, resulting in predominantly plastic deformation. Subsequently, at point B, the process transitions to the shear failure stage, wherein the rate of increase in shear force decelerates. This is due to the fact that, upon entering the plastic deformation stage, the original biological tissue structure of the stem is destroyed. The application of the tool results in the accumulation of fractured tissues, which in turn causes the shear force to continue increasing. However, the accumulated tissues are unable to attain the shear strength of the original biological structure, resulting in a lower rate of force increase in comparison to the plastic deformation stage. At point C, the maximum shear force is reached, after which there is a rapid decrease in shear force as the stem fractures, thereby marking the conclusion of the test.
Figure 4. Analysis of shear fracture experiments. (A) The shear failure physical and numerical simulation experiments of CKS; (B) The FS-x curves of physical and numerical simulation experiments.
A comparison of the FS-x curves derived from the numerical simulation and physical tests reveals that the rate of change in shear force as it transitions from the elastic deformation stage to the linear elastic deformation stage in the numerical shear fracture test is lower than in the physical shear fracture test. This is due to the fact that the DEM model of CKS simplifies the material by not distinguishing between the phloem and xylem, treating it as a homogeneous material. This results in a more gradual change in shear force in the numerical simulation. The discrepancy between the maximum shear force FS values obtained from the physical and numerical simulation tests is 3.74%, as presented in Table 7.
The linear elastic deformation stage is of particular importance for the study of mechanical properties, including the elastic modulus. The F-x curves derived from the numerical simulation and physical tests conducted during this phase were extracted for fitting analysis. The discrepancy in the slope between the numerical simulation and physical tests is 5.21%, indicating a satisfactory fit. The fitted curves and results are presented in the blue box of Figure 4B, where kP represents the average slope of the physical tests and kS denotes the average slope of the numerical simulation tests. This indicates that the mechanical properties of the DEM model are similar to those of the actual CKS.
3.3.2 Analysis of three-point bending physics experiments
Figure 5A shows the three-point bending physical and numerical simulation experiments of CKS. As illustrated in Figure 5B, the three-point bending test commences with an increase in bending force, resulting in the CKS entering the linear elastic deformation stage. During this phase, the CKS is capable of recuperating its deformation under the prevailing load. In an ideal scenario, the force and displacement should be proportional, in accordance with Hooke’s law. However, the force-displacement curve obtained from the experiment did not demonstrate a proportional relationship. This discrepancy may be attributed to the disparate physical properties of the phloem and xylem within the CKS, which renders it a heterogeneous material. However, the stem displayed characteristics consistent with the elastic stage. Subsequently, the CKS entered the yield stage, exhibiting substantial plastic deformation. The accumulation of various tissue cells resulted in a deceleration in the rate of increase in bending force. As the load continued to increase, reaching the maximum bending force, the stem entered the fracture stage. The bending force decreased markedly, resulting in the fracture of the CKS. At this moment, the bending force ceased to increase significantly with displacement, thereby marking the conclusion of the test. The three-point bending characteristics of the CKS are consistent with the test results of other stalks (Réquié et al., 2018). However, there are notable differences in the degree of force and the values of deflection deformation.
Figure 5. Analysis of three-point bending experiments. (A) The three-point bending physical and numerical simulation experiments of CKS; (B) The FB-x curves of physical and numerical simulation experiments.
The two curves exhibit a general consistency in their respective trends. The numerical simulation yielded a value of FB of 111.83 N, exhibiting a 3.32% discrepancy with the target value of 115.67 N. This outcome substantiates the assertion that the DEM model of CKS effectively encapsulates the correlation between bending force and displacement. The overall trend of the numerical simulation curve is in leading of that of the physical test curve. This discrepancy can be attributed to the initial stage of the physical test, wherein the phloem of CKS, characterised by a high moisture content and soft texture, exerts a comparatively diminished influence on the rate of change in bending force (Guo et al., 2021). In contrast, the DEM model employed in the numerical simulation simplifies the material without distinguishing between the phloem and xylem, resulting in a more gradual change in bending force. In both the physical and numerical simulation tests, upon reaching the maximum bending force and entering the fracture stage, the change in bending force no longer follows a regular pattern of decrease. In the physical test, upon reaching the yield limit, the axial fibres of the CKS xylem enter the fracture stage. The accumulation of fibre cells results in a gradual reduction in bending force until the accumulated fibre cells reach their yield limit, at which point there is a sudden decline in bending force. This fracture stage then repeats the aforementioned process, resulting in a stepwise decrease in bending force. In contrast, the fracture stage of the DEM simulation curve demonstrates a more gradual decline in bending force in comparison to the physical test. This is due to the necessity of ensuring a reasonable computation time in the numerical simulation, which entails the simplification of the model and the reduction of the number of parallel bonds between particles to a level that is considerably lower than that observed in real fibre cells. Consequently, the accumulated particles exhibit diminished resistance to bending loads. During the linear elastic phase of the experimental procedure, the slope error of the F-x curve derived from the physical and numerical simulation tests is 2.15%.
3.3.3 Analytical verification of axial tension and radial compression experiments
Figure 6 illustrates the processes of radial compression, while Figure 7 depicts the processes of radial tension in the physical and numerical simulation tests. To guarantee the dependability of the axial tension experiments, the samples were subjected to a bespoke processing procedure. The ends of the samples were shaped into steps for easy clamping, and the middle part of the samples was made into a cylindrical shape (diameter was 1.10 ± 0.04 mm) with a smooth transition to avoid stress concentration (Li et al., 2024).
Figure 6. Analysis of radial compression experiments. (A) Performance of CKS in compression experiment; (B) The FC-x curves of physical and numerical simulation experiments; (C) Change process of total bond numbers.
Figure 7. Analysis of axial tension experiments. (A) Axial tension experiments of CKS; (B) The FT-x curves of axial tension experiments.
As illustrated in Figure 6C, it is evident that during the displacement range of 0-0.6 mm, the number of parallel bonds exhibits a pronounced decline. At this moment, the CKS is undergoing the expansion phase of cracking. Subsequent to 1.0 mm, the reduction in parallel bonds becomes more gradual, and the CKS gradually approaches the limit of its compressive strength. Upon reaching the maximum pressure, the CKS undergoes complete collapse, accompanied by the rupture of the parallel bonds within the model (Zhao et al., 2023).
The F-x curves for radial compression and axial tension in the physical tests exhibited trends that were generally consistent with those reported in previous studies, including those by (Guo et al., 2021; Zhao et al., 2023; Li et al., 2024). The F-x curves derived from the numerical simulation tests exhibited slight discrepancies compared to those obtained from the physical tests. This can be attributed to the use of uniform particles in the numerical model, which did not accurately represent the heterogeneous mechanical properties of the actual phloem and xylem. However, with regard to the maximum failure force, the discrepancy between the physical and numerical simulation tests was 4.37% for radial compression and 4.87% for axial tension, which is within an acceptable range.
3.3.4 Cutting experiments analysis verification
The F-x curves of the physical cutting process and the numerical simulation cutting process demonstrate a similar trend (Figure 8A). In the initial phase of the cutting process, the cutting force exhibited in the physical test displays greater fluctuations than in the numerical simulation. This is attributed to the more intricate internal structure of the stem, which differs from the discrete element model. As the tool advances further into the material, the cutting force reaches a maximum of 271.64 N, after which it declines rapidly. This occurs as the tool transitions from cutting the harder xylem to the softer phloem. However, the tool’s penetration also causes compression, which in turn gives rise to fluctuations in the cutting force during this stage. Subsequently, a second sharp decline in the cutting force is observed once the majority of the stem’s tissue has been cut, marking the complete severing of the stem and the conclusion of the test.
Figure 8. Analysis of cutting experiments. (A) The F-x curves of cutting experiments; (B) Sectional morphology of stems cut with different blades; (C) The progress of numerical simulation cutting test (purple bonds – the bond staus is normal, green bonds - the bond staus is breaking).
Furthermore, the numerical simulation cutting test also exhibits two sharp drops in cutting force. This phenomenon can be attributed to the exertion of a pushing action by the cutting tool on the parallel bonds during the cutting process. This action results in discontinuous fractures of the parallel bonds, preventing a smooth decrease in cutting force. An examination of the status of the parallel bonds reveals a greater number of broken bonds on the flat side in comparison to the inclined side, resulting in a more even and smooth cut on the flat side (Figure 8C). The results of the physical tests demonstrate that the actual CKS exhibits an uneven cross-section following cutting with a bevelled edge in comparison to a flat edge. This finding is in alignment with the performance observed in the numerical simulation tests (Figure 8B). The results indicate that the utilisation of a flat cutting method for stems is more advantageous for the subsequent growth of the stems and for the achievement of a smoother stubble after cutting.
The cutting performance was characterised by the shear stress τ (Equation 8). In this context, Fmax represents the maximum cutting force, while Sa denotes the cross-sectional area of the cut portion of the cutting specimen. The discrepancy between the numerical simulation results and the physical experiment was 3.89%, indicating that the discrete element model parameters had been calibrated reliably for subsequent studies.
3.4 Discussion
In this study, we found that the parallel bond parameter exerted a significant influence on the maximum bending force, while the contact parameter exerted a comparatively lesser effect on the F-x curve for three-point bending. Furthermore, the order of significance of the parallel bond parameters affecting stem failure force was found to be consistent with that reported in previous research (Zhao et al., 2023). Comparing and analysing physical tests with numerical simulation tests, the maximum failure force errors for the three-point bending, shear fracture, and radial compression tests were found to be small, and the F-x curve trends were observed to be consistent, same as previous studies (Sun et al., 2023; Zhao et al., 2023). However, the maximum failure force error for the axial tension test was smaller than that reported by (Guo et al., 2021), and the F-x curve trend was closer. The elastic deformation stage represented a crucial methodology for investigating the biomechanical properties of materials. Numerical simulation tests of three-point bending and shear fracture fitted equations at this stage with high accuracy and small errors between simulated and actual slopes. This provided a robust foundation for further studies on the elastic modulus of CKS. In the cutting test, both the inclined and flat sides of the tool and the parallel bond of the model were damaged. The inclined edge exhibited a tearing effect, resulting in more significant damage to the section. In contrast, the section after cutting by the flat edge displayed a relatively smooth surface (Li et al., 2024). This same performance was observed in the physical test. A comparison of the maximum cutting force and F-x between the physical and numerical simulation cutting tests reveals a small discrepancy. In comparison to other studies, this paper presented a comprehensive analysis of the biomechanical properties of CKS at both the micro and macro levels, and offered verification of the discrete element model from multiple perspectives. Furthermore, the F-x curves of all numerical simulation tests exhibited a leading trend compared to the physical tests. Additionally, the mechanical properties of xylem and phloem could not be adequately characterised by a single granular material. This limitation was to be addressed in future studies through the modelling of xylem and phloem separately. To date, there have been few studies on CKS. The use of DEM numerical simulation methods facilitated an understanding of the failure modes of the mechanical properties of real CKS at a microscopic level, thereby contributing to the subsequent in-depth study of tool-stem interaction mechanisms. It provided theoretical support for the study of CKS harvesting theory and harvesting machinery.
4 Conclusion
1. A DEM model of CKS based on the Hertz-Mindlin contact model was established, the parameters of which were calibrated. The optimal parameter combinations for each significant factor were then determined by Plackett-Burman design and response surface analysis experiments: shear stiffness per unit area is 3.56×109 Pa, the bonded disk scale is 0.93 mm, the normal stiffness per unit area is 9.68×109 Pa, the normal strength is 5.62×107 Pa, the shear strength is 4.27×107 Pa.
2. The numerical simulation tests for three-point bending, radial compression, axial tension, shear fracture, and cutting tests using the CKS model based on the aforementioned parameter combinations results with errors of 3.32%, 4.37%, 4.87%, 3.74%, and 3.89%, respectively, in comparison to their corresponding physical experiments. Additionally, the trend of F-x of the numerical simulation tests is approximately the same as the curve trend of the physical tests.
3. The cutting tests demonstrated that a smaller radial angle between the tool edge and the stem resulted in a more even and smooth cut surface. The DEM model of CKS established in this study effectively represents the biomechanical properties of real CKS, provides a basis for subsequent in-depth study of crop-machine interaction mechanisms, and guides the design and development of cutting devices for CKS harvesting machinery.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
CQ: Conceptualization, Data curation, Methodology, Writing – original draft. ZS: Investigation, Software, Writing – review & editing. LT: Validation, Visualization, Writing – review & editing. ZG: Supervision, Writing – review & editing. YH: Funding acquisition, Project administration, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the project of the National Key R&D Program of China (2022YFD2001901).
Acknowledgments
Thanks to the Hohhot Branch of Chinese Academy of Agricultural Mechanization Sciences, the samples used for the experiment were collected from their experimental field.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2024.1457243/full#supplementary-material
Glossary
A: The effective area of CKS, m2
Ab: The cross-sectional area of the bonds, m2
E: Elastic modulus, Pa
F-x: Force and deformation curve
Faverage: The average value of FB, N
FC: Compression force, N
FB: Bending force, N
FT: Axial tensile force, N
Fn, Fn’: Normal adhesion forces of the bonds, N
FS: Shear force, N
Fs, Fs’: Tangential adhesion forces of the bonds, N
G: Shear modulus, Pa
H: CKS physical experiment drop height, m
H’: CKS DEM model simulation experiment drop height, m
h: CKS physical experiment post-crash rise height, m
h’: CKS DEM model numerical simulation experiment post-crash rise height, m
I: Moment of inertia of bonds, m4
IP: Polar moment of inertia of bonds, m4
J: Moment of inertia of CKS, m4
kn: Normal stiffness per unit area, N·m−3
ks: Shear stiffness per unit area, N·m−3
L: The effective length of CKS, mm
Lb: The length of bonds, mm
ΔL: The change in the effective length, mm
S: Span of three-point bending experiment, mm
S1: CKS physics experiment incline movement distance, m
S1’: CKS DEM model simulation experiment incline movement distance, m
S2: CKS physics experiment plane motion distance, m
S2’: CKS DEM model simulation experiment plane motion distance, m
Vn: Normal velocities of the particles, m·s−1
Vt: Tangential velocities of the particles, m·s−1
ωn: Normal angular velocities of the particles, rad·s−1
ωt: tangential angular velocities of the particles, rad·s−1
σmax: Critical normal stress, Pa
τmax: Critical shear stress, Pa
μ: Poisson’s ratio
ρ: Density, kg·m−3
θ: The angle of the inclined plane of the static friction physical experiment, °
θ’: The angle of the inclined plane for static friction simulation experiment, °
x1: Coefficient of restitution between CKS and Steel
x2: Coefficient of static friction between CKS and Steel
x3: Coefficient of rolling friction between CKS and Steel
δt: Time step
α: The angle of the inclined plane of the dynamic friction physical experiment, °
α’: The angle of the inclined plane for dynamic friction simulation experiment, °
Vn: Normal velocities of the particles, m·s−1
Vt: Tangential velocities of the particles, m·s−1
CKS: Stem of Caragana korshinskiil
DEM: Discrete element method
References
ASTM D854-10 (2010). Standard test methods for specific gravity of soil solids by water pycnometer. Am. Soc. Test. Mater.
Aydin, I., Arslan, S. (2018). Mechanical properties of cotton shoots for topping. Ind. Crops Prod. 112, 396–401. doi: 10.1016/j.indcrop.2017.12.036
Cabiscol, R., Finke, J. H., Kwade, A. (2018). Calibration and interpretation of DEM parameters for simulations of cylindrical tablets with multi-sphere approach. Powder Technol. 327, 232–245. doi: 10.1016/j.powtec.2017.12.041
Chen, J. H., Zhao, N., Fu, N., Li, D., Wang, L. J., Chen, X. D. (2019). Mechanical properties of hulless barley stem with different moisture contents. Int. J. Food Eng. 15, 1–10. doi: 10.1515/ijfe-2018-0033
Chen, Y. F., Dong, S. P., Xing, D. L., Gao, W., Li, J. X., Li, Y., et al. (2020). Current situation and development suggestions of caragana microphylla stumping harvest. Trans. Chin. Soc Agric. Eng. 10, 30–34. doi: 10.3969/j.issn.2095-1795.2020.07.008
Fan, G. J., Wang, S. Y., Shi, W. J., Gong, Z. F., Gao, M. (2022). Simulation parameter calibration and test of typical pear varieties based on discrete element method. Agronomy 12, 1720. doi: 10.3390/agronomy12071720
Fu, M., Chen, X. Q., Gao, Z. F., Wang, C. M., Xu, B., Hao, Y. L. (2023). Parameters calibration of discrete element model for crushed corn stalks. Inmateh-Agricult. Eng. 69, 399–408. doi: 10.35633/inmateh-69-37
Gao, Y. Y., Kang, F., Kan, J. M., Wang, Y. T., Tong, S. C. (2021). Analysis and experiment of cutting mechanical parameters for Caragana korshinskii (C.K.) Branches. Forests 12, 1359. doi: 10.3390/f12101359
Guo, J., Karkee, M., Yang, Z., Fu, H., Li, J., Jiang, Y. L., et al. (2021). Discrete element modeling and physical experiment research on the biomechanical properties of banana bunch stalk for postharvest machine development. Comput. Electron. Agric. 188, 106308. doi: 10.1016/j.compag.2021.106308
Hu, Y., Xiang, W., Duan, Y. P., Yan, B., Ma, L., Liu, J. J., et al. (2023). Calibration of ramie stalk contact parameters based on the discrete element method. Agriculture-Basel 13 (5). 1070. doi: 10.3390/agriculture13051070
Huan, X. L., Wang, D. C., You, Y., Ma, W. P., Zhu, L., Li, S. B. (2022). Establishment and calibration of discrete element model of king grass stalk based on throwing test. Inmateh-Agricult. Eng. 66, 19–30. doi: 10.35633/inmateh-66-02
Jiang, P., Li, Y. P., Li, J. L., Meng, H. W., Peng, X. B., Zhang, B. C., et al. (2021). Experimental research on the bending and fracture characteristics of cotton stalk. Trans. ASABE. 64, 1771–1779. doi: 10.13031/trans.14589
Jiménez-Herrera, N., Barrios, G. K. P., Tavares, L. M. (2018). Comparison of breakage models in dem in simulating impact on particle beds. Adv. Powder Technol. 29, 692–706. doi: 10.1016/j.apt.2017.12.006
Khan, M. M. R., Chen, Y., Lagüe, C., Landry, H., Peng, Q. J., Zhong, W. (2010). Compressive properties of Hemp (Cannabis sativa L.) stalks. Biosyst. Eng. 106, 315–323. doi: 10.1016/j.biosystemseng.2010.04.004
Kruggel-Emden, H., Rickelt, S., Wirtz, S., Scherer, V. (2008). A study on the validity of the multi-sphere discrete element method. Powder Technol. 188, 153–165. doi: 10.1016/j.powtec.2008.04.037
Leblicq, T., Smeets, B., Ramon, H., Saeys, W. (2016a). A discrete element approach for modelling the compression of crop stems. Comput. Electron. Agric. 123, 80–88. doi: 10.1016/j.compag.2016.03.022
Leblicq, T., Smeets, B., Vanmaercke, S., Ramon, H., Saeys, W. (2016b). A discrete element approach for modelling bendable of crop stems. Comput. Electron. Agric. 124, 149–88. doi: 10.1016/j.compag.2016.03.022
Lei, J. L., Wang, S. W., Lei, D. Y., Liu, Z. C. (2022). Measurement and calibration of discrete element simulation parameters of crushed sugarcane tail leaves. Bioresources 17, 5984–5998. doi: 10.15376/biores.17.4.5984-5998
Li, S. B., Huan, X. L., Wang, T. Y., Hui, Y. T., You, Y., Wang, D. C. (2024). Biomechanical properties and discrete element modeling of PSR stalks during silage harvest. Comput. Electron. Agric. 217, 108644. doi: 10.1016/j.compag.2024.108644
Liang, H. B., Xue, Y. Y., Shi, J. W., Li, Z. S., Liu, G. H., Fu, B. J. (2019). Soil moisture dynamics under Caragana korshinskii shrubs of different ages in Wuzhai County on the Loess Plateau, China. Earth Environ. Sci. Trans. R. Soc. Edinburgh. 109, 387–396. doi: 10.1017/s1755691018000622
Liao, Y. T., Liao, Q. X., Zhou, Y., Wang, Z. T., Jiang, Y. J., Liang, F. (2020). Parameters calibration of discrete element model of fodder rape crop harvest in bolting stage. Trans. Chin. Soc Agric. Mach. 51, 73–82. doi: 10.6041/j.issn.1000-1298.2020.06.008
Liu, F. Y., Zhang, J., Chen, J. (2018). Modeling of flexible wheat straw by discrete element method and its parameter calibration. Int. J. Agric. Biol. Eng. 11, 42–46. doi: 10.25165/j.ijabe.20181103.3381
Liu, Q., Mathanker, S. K., Zhang, Q., Hansen, A. C. (2012). Biomechanical properties of Miscanthus stems. Trans. ASAE. 55, 1125–1131. doi: 10.13031/2013.42231
Liu, W. R., Zhang, G. Z., Zhou, Y., Liu, H. P., Tang, N. R., Kang, Q. X., et al. (2022). Establishment of discrete element flexible model of the tiller taro plant and clamping and pulling experiment. Front. Plant Sci. 13, 1019017. doi: 10.3389/fpls.2021.635440
Luo, W. W., Hu, Z. C., Wu, F., Gu, F. W., Xu, H. B., Chen, Y. Q. (2019). Design and optimization for smashed straw guide device of wheat clean area planter under full straw field. Trans. Chin. Soc Agric. Eng. 35, 1–10. doi: 10.11975/j.issn.1002-6819.2019.18.001
Luo, K., Wu, Z. M., Cao, C. M., Qin, K., Zhang, X. C., An, M. H. (2022). Biomechanical characterization of bionic mechanical harvesting of tea buds. Agriculture-Basel 12 (9), 1361. doi: 10.3390/agriculture12091361
Ma, F. L., Lekai, L. K., Wang, Y. T., Li, P., Zhu, C. W. (2020). Biomechanical properties of wolfberry plant organs. Engenharia Agricola. 40, 162–176. doi: 10.1590/1809-4430-Eng.Agric.v40n2p162-176/2020
Madison, W. (2010). Wood handbook-woods as an engineering material [M]. General Technical Report FPL-GTR-190. WI: U.S. Department of Agriculture, Forest Service, Forest Products Laboratory.
Ren, J. H., Liu, X. L., Yang, W. P., Yang, X. X., Li, W. G., Xia, Q., et al. (2021). Rhizosphere soil properties, microbial community, and enzyme activities: Short-term responses to partial substitution of chemical fertilizer with organic manure. J. Environ. Manage. 299, 1136501. doi: 10.1016/j.jenvman.2021.113650
Réquié, S., Goudenhooft, C., Bourmaud, A., Le Duigou, A., Baley, C. (2018). Exploring the link between flexural behaviour of hemp and flax stems and fibre stiffness. Ind. Crops Prod. 113, 179–186. doi: 10.1016/j.indcrop.2018.01.035
Roy, R., Mostofa, M. G., Wang, J. X., Fornara, D., Sarker, T., Zhang, R. Q. (2021). Revegetation intervention of drought-prone coal-mined spoils using Caragana korshinskii under variable water and nitrogen-phosphorus resources. Agric. Water Manage. 246, 106712. doi: 10.1016/j.agwat.2020.106712
Schramm, M., Tekeste, M. Z. (2022). Wheat straw direct shear simulation using discrete element method of fibrous bonded model. Biosyst. Eng. 213, 1–12. doi: 10.1016/j.biosystemseng.2021.10.010
Sun, K., Yu, J. Q., Zhao, J. W., Liang, L. S., Yu, Y. J. (2023). A DEM-based general modeling method and experimental verification for wheat plants in the mature period. Comput. Electron. Agric. 214, 108283. doi: 10.1016/j.compag.2023.108283
Wang, Y., Zhang, Y., Yang, Y., Zhao, H. M., Yang, C. H., He, Y., et al. (2020). Discrete element modelling of citrus fruit stalks and its verification. Biosyst. Eng. 200, 400–414. doi: 10.1016/j.biosystemseng.2020.10.020
Wang, Y., Zhao, H. L., Zhao, X. Y. (2013). Effects of land use intensity on the restoration capacity of sandy land vegetation and soil moisture in fenced sandy land in desert area. Contemp. Prob. Ecol. 6, 128–136. doi: 10.1134/S1995425513010198
Xia, R., Li, B., Wang, X. W., Li, T. J., Yang, Z. J. (2019). Measurement and calibration of the discrete element parameters of wet bulk coal. Measurement 142, 84–95. doi: 10.1016/j.measurement.2019.04.069
Xia, Y. D., Klinger, J., Bhattacharjee, T., Thompson, V. (2022). The elastoplastic flexural behaviour of corn stalks. Biosyst. Eng. 216, 218–228. doi: 10.1016/j.biosystemseng.2022.02.016
Yang, Q., Zhang, T., Wang, Y., Li, G., Yin, J. J., Han, X. M., et al. (2013). Construction of a suppression subtractive hybridization library of Caragana korshinskii under drought stress and cloning of ckWRKY1Gene. scientia Silvae sinicae. 49, 62–68. doi: 10.11707/j.1001-7488.20130709
Zeng, Z. W., Chen, Y. (2019). Simulation of straw movement by discrete element modelling of straw-sweep-soil interaction. Biosyst. Eng. 180, 25–35. doi: 10.1016/j.biosystemseng.2019.01.009
Zhang, F. W., Song, X. F., Zhang, X. K., Zhang, F. Y., Wei, W. C., Dai, F. (2019). Simulation and experiment on mechanical characteristics of kneading and crushing process of corn straw. Trans. Chin. Soc Agric. Eng. 35, 58–65. doi: 10.11975/j.issn.1002-6819.2019.09.007
Zhang, H. Q., Tang, M., Chen, H., Tian, Z. Q., Xue, Y. Q., Feng, Y. (2010). Communities of arbuscular mycorrhizal fungi and bacteria in the rhizosphere of Caragana korshinkii and Hippophae rhamnoides in Zhifanggou watershed. Plant Soil. 326, 415–424. doi: 10.1007/s11104-009-0022-1
Zhang, Z. G., Xu, H. W., Xue, H. T., Ni, C., Wang, C. (2023). Calibration and experimental of discrete element parameters of panax notoginseng stem. Trans. Chin. Soc Agric. Mach. 54, 61–70 + 91. doi: 10.6041/j.issn.1000-1298.2023.11.006
Keywords: discrete element method, stem of Caragana korshinskii Kom., mechanical properties, parameter calibration, stem cutting
Citation: Qingqiu C, Shengwei Z, Tao L, Gaixia Z and Hongfang Y (2024) Discrete element modelling and mechanical properties and cutting experiments of Caragana korshinskii Kom. stems. Front. Plant Sci. 15:1457243. doi: 10.3389/fpls.2024.1457243
Received: 17 July 2024; Accepted: 07 October 2024;
Published: 31 October 2024.
Edited by:
Mark Blyth, University of East Anglia, United KingdomReviewed by:
Changsu Xu, Southwest University, ChinaLiqing Chen, Anhui Agricultural University, China
Copyright © 2024 Qingqiu, Shengwei, Tao, Gaixia and Hongfang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yuan Hongfang, eWhmMTk4NDgyOEAxNjMuY29t