Modeling and experiments of chewing mechanical properties of pellet feed using discrete element method

In order to explore the mechanical properties and breaking behavior of pellet feed during chewing, the experiments of texture mechanics, followed by the modeling and simulation of pellet feed based on the discrete element method were carried out in this research. Five wet basis moisture contents (8%, 10%, 12%, 14% and 16%, respectively), two kinds of loading directions (L and D direction, respectively) of pig pellet feed were selected as variables. First, mechanical parameters including hardness, elasticity, tackiness and chewiness were measured by a texture analyzer. The results of compression tests showed that the hardness of pig pellet feed was 5.44-32.43 N, the elastic index was 0.04-0.94 mm, the tackiness was 0.076.63 N, the chewiness was 5.52-27.39 mJ. Moreover, the hardness, tackiness and chewiness of pig pellet feed decreased significantly with the increase of moisture content but the elasticity showed an adverse varying trend. The hardness, elasticity, tackiness and chewiness along D direction outweighed that of L direction in numerical data at the same moisture content. Then, the chewing and breaking processes of pellet feed were simulated based on the discrete element method (DEM) combined with bonding particle model, in which the whole pellet feed were considered as agglomerations of micro-particles and broke when the stress between micro-particles had exceeded the maximum limit. Multi-parameter optimization experiments were carried out using quadratic orthogonal rotation design, in which stiffness coefficient (X1), critical stress (X2) and bonding radius (X3) were the influencing factors, hardness (Y1), elasticity (Y2), tackiness (Y3) and chewiness (Y4) were evaluating indicators. Based on the regression analysis of the Design-Expert 8.0.6 software and response surface analysis method, the relationship between the three influencing factors and evaluating indicators was established. The similarity between experimental and simulated results in feed morphology and mechanical index proved that the modeling method for pellet feed based on DEM was effective and accurate. This work can provide a reference for the feed forming process and the optimization design of the related feed machinery. Meanwhile, the DEM model provided a new method for evaluating the texture and palatability of pellet feed.


Introduction 
The feed is the material basis of livestock and aquaculture as the main food source. China has become the largest feed producer of the world and the total feed output has reached 228 million tons by the end of 2018 [1] . Pellet feed is the most important feed form and has a wide range of applicability. Compared with mash feed, pellet feed has virtues of low animal pickiness, high feed return rate, anti-harmful bacteria, as well as convenience in storage, management and transportation, which assure it a broad application market and research prospects [2,3] . The production process of pellet feed is as follows: first, raw materials are crushed and mixed, followed by moisture and heat modulation. Next, a mechanical extrusion is performed during which the treated materials are forced to polymerize through the die hole. Therefore, the physical essence of pellet feed is the aggregation of micro-particles.
Chewing and crushing characteristics of pellet feed will directly affect animals' intake, digestion and absorption of nutrition and energy. Therefore, exploring the chewing process of pellet feed is helpful in analyzing the animals' demand for nutrition and taste, as well as evaluating feed palatability, providing a reference for feed processing technology and equipment study. Above all, it is of great theoretical and practical significance to analyze the texture of feed and study its chewing and breaking mechanism.
In recent years, research has been conducted to investigate the texture test of materials [4][5][6] , such as crop stems [7] , rapeseed [8] , biscuits [9] , dough [10] , etc. The results show that the moisture content and loading mode are important factors affecting the texture properties of materials [4,9,10] . They provided a theoretical basis for model establishment as well as research methods to explore the texture mechanical properties of pellet feed [5][6][7][8] .
Chewing is a complex dynamic process and numerical simulation is helpful in studying the relationship between the physical properties of the feed and chewing process. Sun et al. [11] used a two-dimensional PFC software to establish a feed chewing and swallowing model, followed by simulating the first bite and swallowing process of three shapes of feed (rectangular, circular and hollow circular, respectively). It showed that the discrete element method is feasible to simulate the food chewing process but only limited within a two-dimensional plane, the mechanical model of three-dimensional space has not been established. Sun et al. [12] used finite element software ANSYS to simulate the damage effect of three different indenters on the same feed material and analyzed the stress nephogram on the surface of the material after chewing.
However, the results only referred to the deformation of the material and could not reproduce and characterize the fracture process of the material, thus the simulation deviates from the real chewing process.
DEM which is based on the principle of molecular dynamics is a method to analyze discrete particle materials proposed by Professor Cundall in 1971 [13] . It has been widely used in the study of powder dynamics in complex physical fields, material mechanical properties with more complex structure and multiphase mixed material medium. At present, it is also used in compaction and crushing of fertilizer [13] , soil [14] , wheat straw [15] , rock [16] , etc. It is known that feed is formed by extruding and agglomerating powder particles with a discrete structure, and chewing is the dynamic process of feed breaking and moving. Therefore, it is more practical to establish a model of feed agglomeration by using DEM. However, most studies only use repose angle as an evaluation index [17,18] for parameter calibration in DEM, with drawbacks of only represent the flow characteristics of DEM but failed to reflect the physical and mechanical characteristics of the material itself.
In order to solve this problem, this paper chose the pellet feed with different moisture contents as experimental material. First, the compression mechanical tests were carried out by the texture analyzer and the mechanical properties were analyzed. Then, the chewing mechanical model of the pellet feed was constructed based on the discrete element bonding mechanics theory. Finally, the mechanical and morphological changes of the pellet feed were analyzed by simulating the animal chewing process. The results will improve the palatability of pellet feed for animals and feeding rates, thus promote the healthy growth of animals.

Sample treatment and measure methods
The raw material was pig feed produced by Beijing Xinsanfeng feed factory. The pellet feed with a uniform diameter and a typical cylinder structure was extruded from the die hole of the die roller pellet mill. The geometric sizes of 50 randomly selected pellets were measured by the digital vernier caliper (accuracy is 0.01 mm, Zhangjiakou Jinfeng hardware tool manufacturing Co., Ltd.) and the average is 3.5 mm in diameter and 6.8 mm in length. Considering the influence of length diameter ratio on the strength of pellet feed, the pellet feed was rubbed with a 100-mesh fine sandpaper until they have a uniform length and two flat ends that allowed straight standing on the plate.

Test instrument
The moisture content of pellet feed was measured by Pl2002 electronic balance (Shanghai Mettler Toledo Instrument Co., Ltd.) and DHG-9240a electrothermal constant temperature blast drying oven (Shanghai Jinghong Experimental Equipment Co., Ltd.). The texture mechanics were determined by the TMS pilot texture instrument (FTC company, United States) with a maximum load of 100 N and a sensor accuracy of 0.01 N.

Preparation of test samples
Generally, the moisture content of pellet feed is 11%±1% when it leaves the factory. In order to study the effect of moisture content on the chewing mechanical properties of the pellet feed, the pellet feed was treated with water by a drying and water refilling method [19] . The drying was carried out at 105°C according to ASAE [20] . Samples at 8% and 10% were prepared by drying a known mass of pellet feed at an initial MC to the pre-calculated weight in a hot air oven. Samples with MC of 12%, 14% and 16% were prepared by an adjustment method specified as follows: First, calculate the mass of distilled water that needs to be added using Equation (1), then spray the distilled water evenly onto the pellet feed. Finally, put the watered pellet feed in a sealed bag for 24 h to make the moisture even. The moisture content of pellet feed can be adjusted by the above procedure to 8%, 10%, 12%, 14% and 16%, respectively.
() (100 ) where, Q is mass of distilled water to be added, g; w imass of feed, g; m i is initial moisture content of feed, %; m f is moisture content of feed after adjustment, %.

Chewing mechanical properties of pellet feed
Texture reflects hardness, elasticity and chewiness of food, thus it is an important index for food quality evaluation. Texture profile analyzer (TPA) determines the chewing behavior of teeth by continuously repeating the mechanical compression of samples. It is mainly used to characterize the physical properties of many kinds of food, such as rapeseed [21] , meat [22] , seafood [23] , maize [24] and etc. In this part, parameters of texture mechanical properties including hardness, brittleness, viscosity, cohesion, elasticity, stickiness, resistance and recovery will be analyzed by simulating the chewing movement of teeth. The palatability can be evaluated using the texture test curve of the material given by the microcomputer. As the pellet feed had a shape of a typical cylinder, the axial and radial pressure loading of pellet feed was carried out respectively ( Figure  1), to ensure the reliability and accuracy of the tests. In the test, a rigid flat indenter with a bottom diameter of 30 mm was used. The pellet feed was placed under the center of the indenter. TPA was connected to the computer and controlled by the texture loader software. The specific steps were described as follows: first, turned on the texture analyzer and preheated for 5 min.
Then chose the ''TPA'' module and secondary compression mode. Set parameters: speed of the indenter 2 mm/s, minimum starting force 0.2 N, compression displacement 30% and sampling frequency 0.02 s.
The force, displacement and deformation were recorded simultaneously. 10 samples were tested for each moisture content and the average values were calculated.

Test results and analysis 2.2.1 Analysis of variance
Hardness represents the pressure born per unit area of feed. A certain degree of hardness can guarantee the integrity of feed during the transportation and feeding process. Elasticity refers to the degree of feed recovery after the first compression deformation with no external force. Tackiness reflects the cohesion of feed, which is the energy performance needed to break feed before swallowing. Chewiness refers to the work done in the process of simulating animal mastication and reflects the continuous resistance of feed during mastication. Texture reflects feed quality and palatability, thus has an important influence on feed and process control.

Influence of various factors on mechanical parameters
TPA results were shown in Table 1 (average values were obtained by at least 5 parallel samples). The data statistics software SPSS was used to make a variance analysis on the test results of feed texture mechanical parameters. The dependent variables were hardness, elasticity, tackiness and chewiness. The fixed factors were loading mode and moisture content.  As can be seen from the analysis in Table 2, the loading mode has a significant effect on hardness (p<0.01), elasticity (p<0.01) and chewiness (p<0.01) of pellet feed. The influence of moisture content on hardness (p<0.01), elasticity (p<0.01) and chewiness (p<0.01) of pellet feed were very significant, but the influence of moisture content on tackiness was not significant (p=0.344). The elasticity increased with the increase of moisture content, but hardness, viscosity and chewiness decreased with the increase of moisture content in spite of loading modes. As can be seen in Table 1 [25] and Ma et al. [26] , respectively. The reason may be due to the softening of feed internal tissue and decrease of adhesive force between micro-particles as moisture content increased. The ability to bear the load and resist crushing weakened, resulting in the decrease of hardness, adhesive viscosity and chewiness of pellet feed. When the moisture content is higher than 16%, the feed shows liquid slurry characteristics and its hardness drops sharply, so it is difficult to form pellet in this situation.

3
Construction of pellet feed model and its parameter check 3.1 DEM model DEM was proposed by Cundall and Strack in 1979 and was mainly used in the field of geotechnical mechanics. Based on this, Potyondy and Cundall developed the bonded particle model (BPM) in 2004, as shown in Figure 2. BPM is based on bonding dense/loose arbitrarily sized circular (2D) or sphere (3D) particles together by virtual elastic bonds at their contact points to form a breakable agglomerate.
The interparticle bond can transfer translational and rotational motion, as well as forces and torques between particles [27] . In addition, each bond can resist tensile and shear deformation to some extent. The mechanical properties of materials can be described by force (F i ) and moment (M i ) at the contact point. Newton's second law was used to determine translation and rotation caused by contact force, external force and field force that acting on the particles. The force displacement law was used to calculate the contact force at each contact point due to relative motion. When reached the maximum normal and tangential shear stresses, the bond breaks and the contact action was solved by Hertz Mindlin (no slip) model. M are normal and tangential moments, respectively, N· m; Fi,b is the force between particle A and particle B, N; ni and τi are normal and tangential components respectively; Rb is a bonding radius, mm; Lb is the overlapping amount of particle A and particle B, mm.

Figure 2 BPM contact model
The forces and moments increments within a time step acting on the bond are calculated according to the following formulas [27,28] : where, δt is time step, v n is the normal velocity of particles; v t is the tangential velocity of particles; ω n is the normal angular velocity; ω t is tangential angular velocity and R B is bonding radius. When normal and tangential shear stress exceed the preset threshold, the bond will break according to the following equations: where, δ max is normal critical stress, Pa; τ max is tangential critical stress, Pa.

Parameter setting
The model of feed chewing and breaking was built and the chewing process was simulated based on DEM. The core idea was to replace the pellet feed with micro-particles [28] . Firstly, a cylinder was defined and was moved out of the calculation domain by using particle volume. Then, position information of the cylinder was recorded, and micro-particles with a certain particle size were closely generated around the center of the position. After replacing the cylinder, the particle group was triggered by the bond model.
Thus, micro-particles were bonded to form agglomerate to characterize pellet feed. When the normal stress and tangential stress between micro-particles exceed the critical value, the bond failed and agglomerate began to break.
The pellet feed is made of bonding particle system by extrusion. From the view of physical structure composition, its overall structure is uniform and isotropic. The average particle size of micro-particles was selected as the calculation value in the simulation. The feed with initial moisture content is about 10% thus it is selected as the simulation object. The particle size distribution of the powder before forming pellet was determined by screening analysis and the average size of powder was 821 μm. Then, the discrete element method combined with bonding particle model of pellet feed was established, thus the overall structure and mechanical properties of micro particle groups were close to the actual situation of pellet feed. The physical parameters selected for simulation materials were as close to the actual situation as possible with a density of 860 kg/m 3 , the elastic modulus of 3×10 3 MPa [28,29] . Stiffness coefficient of 1×10 6 N/m, critical stress of 40 kPa and bonding radius of 4.4 mm were used to build pellet feed model. Figure 3 showed a bonded particle model, which was an agglomerate composed of 832 spherical micro-particles with an individual radius of 0.4 mm, forming a cylinder shape with 6.5 mm in height and 6 mm in diameter. In the a. Elevation view b. Vertical view Figure 3 Bonded particle model of pellet feed in the DEM the simulation, a single particle itself cannot be broken down. However, the interaction forces between adjacent particles could be destroyed by the external force and the coordinate information of each particle was recorded for subsequent simulation and calculation.

Simulation analysis and results
The upper and lower teeth used in the chewing process were simplified into two circular plates, which were separately located above and below the pellet feed in Figure 4. The chewing model was basically consistent with the texture experiment. The chewing process was shown in Figures 4b-4d. To be specific: first, feed particles were replaced with bonded micro-particles of 0.4 mm in radius, as shown in Figure 4b. Then, upward moving the lower plate to squeeze large feed particles. The bonds were gradually destroyed with clusters breaking, as shown in Figure 4c. Finally, the lower plate was returned to the initial position ( Figure  4d), completing a chewing cycle. Based on the above parameters, feed chewing and crushing models were set up. The data were collected in the DEM and curves were drawn to obtain the change of cluster force. The maximum instantaneous pressure of pellet feed during crushing was called hardness of pellet feed. Figure 5 shows the typical compression force-displacement response obtained from experimental and simulated tests. The abscissa is the vertical displacement of the top plate and the ordinate is the vertical component of the reaction force between pellet feed and top plate. The peak load (highest ordinate) corresponds to the pellet feed's failure point and represents the value of hardness.

Simulation parameters
In order to characterize the BPM model of the pellet feed, bonding parameters including stiffness coefficient, critical stress and bonding radius were determined [30] .
Other physical parameters selected by simulation were as close as possible to the actual situation. Based on the principle of orthogonal rotation combination test, the stiffness coefficient (X 1 ), critical stress (X 2 ) and bonding radius were used as test variables and factor level table was established as Table 3 (x 1 -x 3 was the true value of each variable and X 1 -X 3 was the coded value of each variable).

Test results and analysis
In discrete element simulation based on BPM, the level coding value of each influencing factor was taken as independent variable and hardness (Y 1 ), elasticity (Y 2 ), tackiness (Y 3 ) and chewiness (Y 4 ) were obtained from simulated results as evaluation indexes. Table  4 listed the design matrix and simulation results with a moisture content of 10%.

Optimization of regression model
The design expert software was used to analyze the response surface map by fixing one of the three factors to level 0 and the influence of the other two factors on pellet feed hardness and chewiness was investigated. The results are shown in Figures 6  and 7, from which the interaction effect between two parameters can be obtained directly. From Figures 6a-6b and Figures 7a-7b, it can be seen that hardness and chewiness of pellet feed showed an increasing trend with the increase of stiffness coefficient (X 1 ), which is associated with the resistance for pellet feed deformation and crushing.
From Figures 6a, 6c, 7a and 7c, it can be seen that hardness and chewiness presented an increasing trend with the increase of critical stress (X 2 ). The bond between adjacent particles will break when the normal force or tangential force exceeds critical stress. This may due to the fact that high critical stress results in an increase of contact force between micro-particles, which makes pellet feed less likely to break. The greater bond stiffness between particles, the stronger particle bond, thus greater load is required when the failure occurs, which means a higher value in hardness and chewiness. Liu et al. [31] had similar conclusions on the uniaxial compressive strength of concrete blocks under different critical stresses. From Figures 6b, 6c, 7b and 7c, it can be seen that hardness and chewiness increased with the increase of the bonding radius (X 3 ). This can be related to the contact surface between particles, which is in a positive correlation with X 3 . It is believed that large contact surface ensures a more efficient transfer of force and moment between adjacent particles, resulting in higher bonding force between particles. As discussed above, high bonding force guaranteed higher value in hardness and chewiness, making pellet feed less likely to be compressed or damaged.
Based on response surface method and multi-objective parameter optimization in the range of -2 ≤ X i ≤ 2 (i = 1, 2, 3), the optimum BPM parameters of pellet feed was: X 1 = 1.65, X 2 = −0.89 and X 3 = −0.35, namely, the stiffness coefficient was 1.66×10 6 N/m, critical stress was 31.1 kPa and bonding radius was 3.7 mm. a.
b. c. Figure 6 Response surface plots of hardness with moisture content of 10% a. b. c. Figure 7 Response surface plots of chewiness with moisture content of 10%

Test verification
Based on the optimum BPM parameters above, the simulated values of hardness and chewiness for pellet feed were (14.79± 0.23) N and (13.41±0.14) mJ, respectively with the moisture content of 10%. In order to estimate the accuracy of the simulation, the simulated results were compared with the experiment results. The error (δ) can be calculated by: where, θ 1 is simulated values; θ 2 is the experimental value. Simulated and experimental results were shown in Tables 6  and 7. Further analysis showed that simulated values were slightly higher than the experimental one because the latter was more complex. In fact, hardness and chewiness are greatly affected by factors such as flatness of the upper/lower end surfaces, internal components, density distribution and local concentrated stress. Once tiny cracks formed, the whole pellet feed was more easily to reach the breaking strength. It can be seen from the above analysis that the BPM model could be used to predict and simulate the chewing and breaking process of pellet feed.
Similarly, DEM simulation and parameter optimization with different moisture contents were carried out. The stiffness coefficient, critical stress and bonding radius were reduced with the increase of moisture content. So as hardness and chewiness, which also decreased as moisture content increased.
The maximum difference between simulated and experimental values was 5.49% for the case of 10% moisture content, indicating that the model fitting effect was good with different moisture contents.

Conclusions
(1) The hardness, elasticity, stickiness and chewiness of pellet feed with different moisture content (8%, 10%, 12%, 14% and 16%, respectively) and loading mode (L and D direction) were measured by compression experimental test. The mechanical characteristic parameters of pellet feed as a function of moisture content and loading mode were analyzed. Specifically, the hardness was 5.44-32.43 N, elasticity was 0.04-0.94 mm, tackiness was 0.07-6.63 N and chewiness was 5.52-27.39 mJ. Moisture content had a significant effect on the mechanical parameters of pellet feed (p<0.01). The elasticity increased with the increase of moisture content, however hardness, tackiness and chewiness decreased with the increase of moisture content. The loading mode had a significant effect on hardness, elasticity, tackiness and chewiness (p<0.01). The values of the texture mechanical parameters in radial loading direction were higher than that in axial loading direction at the same moisture content.
(2) The discrete element method combined with bonding particle model was proposed to study the chewing process of pellet feed.
By building bonds between micro-particles, a three-dimensional model of pellet feed was obtained, which was similar to real pellet feed and could be used to simulate the crushing behavior of feed. Through regression analysis and response surface analysis of test results by design expert software, the optimum parameters were obtained: stiffness coefficient was 1.66×10 6 N/m, critical stress was 31.1 kPa and bonding radius was 3.7 mm. The effectiveness of simulation and regression model was verified by comparing with the experimental values. It provided a new way to establish and characterize the discrete element mechanical model of granular materials like feed and study the breaking mechanism.