# Finite Element Simulation and Performance Test of Loading and Mixing Characteristics of Self-Propelled Total Mixed Ration Mixer

Authors: Fuyang Tian; Yuhua Chen; Zhanhua Song; Yinfa Yan; Li Fade; Zhonghua Wang; Benhai Xiong
Journal: Journal of Engineering
 (2020)

---

## Abstract

Simulation analysis and parameter optimization are performed for the loading and mixing devices of a self-propelled total mixed ration mixer. To reveal the three-dimensional movement of silage material under the action of the loading cutter roller, the latter is modeled using SolidWorks software. ANSYS/LS-DYNA software is used to simulate the process of silage cutting, which is modeled using smoothed particle hydrodynamics coupled with the finite element method. The cutting force and power consumption are simulated, and the behavior of the equivalent strain of the silage is determined. The results showed that silage was broken up mainly by extrusion and shear force due to the loading cutter roller. The power consumption according to the simulation is consistent with the value from an empirical formula, confirming the validity of the proposed modeling method. To study the mixing performance and obtain the optimum parameters of the mixing device, the Hertz–Mindlin model is used for the interaction between material particles and mixing device. A three-factor, five-level method is used to optimize the mixing performance. Material-mixing time, loading rate, and auger speed are chosen as experimental factors and mixed uniformity as an evaluation index. It is found that auger speed and material mixing time have significant effects on mixing uniformity. These results provide reference values allowing the analysis of the crushing of silage and selection of the optimum parameters for mixing performance.


---

## Body

## 1. Introduction

A total mixed ration (TMR) is a nutrient-balanced diet based on the physiological needs of ruminants (cattle, sheep, deer, etc.) at different growth stages. The formula for a TMR is designed by animal nutrition experts and includes roughage chopped to an appropriate length, concentrates, and various additives that are fully mixed in certain proportions [1–3]. Outside China, research and development work on ration mixers has been carried out for many years. There has been a trend towards automation and serial operation of both vertical and horizontal structural types of the mixer, while research has focused mainly on processing performance and on methods for testing the uniformity of mixing [4–10]. In China, there has been a concentration on the development of fixed or trailed vertical single-auger and horizontal double-auger ration mixers, together with the required auxiliary mechanical equipment (e.g., silage loaders, belt conveyors, and dischargers), while research has mainly concerned improvements to the design of existing fixed and trailed ration mixers [11–17]. However, there have been few reports in the literature on the loading and mixing working characteristics of self-propelled TMR mixers [4, 7, 11]. Chinese dairy farms currently need to be equipped with several single-function pieces of mechanical equipment to carry out the task of loading material from silage pits and then feeding it to cattle as part of their rations. This leads to problems of high transportation cost, excessively complex working procedures, and large losses of nutrients (represented by losses of silage from the silage pit). Therefore, given the present situation of animal husbandry in China, research on the mechanism of the processing of rations by self-propelled TMR mixers is of great significance for accelerating the use of TMR feeding technologies.The cutting of material in a silage pit is a complex process, involving various highly nonlinear phenomena, including large deformations and crushing of loose material. Its simulation is, therefore, not an easy task. Smoothed particle hydrodynamics (SPH) is a meshless Lagrangian algorithm that has obvious advantages over the finite element method (FEM) in solving highly nonlinear problems such as those involving large deformations and free surfaces [18–21], but it is not as good as FEM in terms of calculational accuracy and dealing with boundary conditions. In this paper, therefore, a coupled SPH-FEM algorithm is proposed that provides a new and effective method for studying the interaction between silage material and a loading cutter roller. In this approach, the FEM is used to model and calculate the small-deformation zone in the cutting area of silage material and SPH is used to model and calculate the large-deformation zone where the mesh is prone to distortion. This combination improves calculational efficiency while ensuring the accuracy of the simulation results, thus maximizing the advantages of both algorithms [20, 22]. The interaction between the granular silage materials and the rigid components of the cutting machinery is modeled using the discrete element method (DEM), which is a numerical simulation method for discontinuous medium problems [23–29], based on the dynamic relaxation method combined with Newton’s second law for various constitutive relations.The work in this paper is based mainly on the finite element software ANSYS/LS-DYNA, using an SPH-FEM coupling algorithm to simulate the cutting of silage material by the loading cutter roller. By analyzing the equivalent strain distribution of silage material at different times during the simulation, the mechanism of the interaction between the silage material and the loading cutter roller is revealed. The DEM software (EDEM) is used to simulate the material mixing process of the mixing device, explore the motion of material particle groups in the bin, and analyze the influence of material mixing time, loading rate, and auger speed on the material mixing performance. The results are expected to provide some guidance regarding the optimal combinations of the material-mixing parameters.


## 2. System Design

### 2.1. Structure

The self-propelled TMR mixer consists of a diesel engine unit, a loading device, a conveyor, a mixing device, and a discharging device. A two-dimensional sketch is shown in Figure1. The power for the whole machine is provided by the diesel engine, which drives a hydraulic pump. By operating a valve, the direction of the hydraulic oil flow can be changed, thus changing the direction of motion of the motors driving the road wheels and allowing the whole machine to move forward or backward. The loading device consists principally of a roller and cutters. During the loading operation in the silage pit, rotary cutters bolted to a spiral blade shift the cut silage material to the center of the loading device, and the material is quickly thrown onto the conveyor to prevent blockages. The mixing device consists principally of a vertical spiral auger and a bin and is used to mix the crude fodder and concentrate feeds.Figure 1
Structure of self-propelled TMR mixer: 1, land wheel; 2, roller; 3, cutter; 4, conveyor; 5, lifting hydraulic cylinder; 6, discharging hydraulic cylinder; 7, moving blade; 8, auger; 9, bin; 10, diesel engine unit; 11, discharging device; 12, road wheels; 13, cab.The self-propelled TMR mixer was designed to realize integrated feeding in the pasture with high efficiency. Self-propelled TMR mixers like the one considered here have played an important role in speeding up the utilization of crop straw as feed. This constitutes one component of the implementation of a policy of “grain to feed” with the aim of alleviating the problems arising when people and livestock compete for grain.


### 2.2. Operating Procedure

The self-propelled TMR mixer realizes the combined functions of loading and conveying of material in the silage pit, mixing of roughage and concentrate, and transporting and discharging of rations. 3D rendering of the whole machine is shown in Figure2.Figure 2
3D rendering of the self-propelled TMR mixer.When the machine is working, a variety of feeds are loaded, depending on the dietary formulation. The lifting action of the loading device on the material in the silage pit and the motion of the conveyor are realized by the reciprocating motion of the piston rod of the double-acting lifting hydraulic cylinder (Figure1 (5)). The cutting of the material in the silage pit is carried out by the rotating loading cutter roller (Figure 1 (2)), and the cut material is then transported to the bin (Figure 1 (9)) of the mixing device by the conveyor (Figure 1 (4)). Driven by a hydraulic motor, a spiral blade radially distributed along the auger (Figure 1 (8)) rotates with the central shaft, and the material moves along the circumferential and axial directions of the spiral. When this material has piled up to a certain height, it falls away freely. As the spiral blade passes through the material, the latter is subjected to inertial and frictional forces, which can cause three-dimensional mixing in the axial, circumferential, and radial directions. This mixing process is mainly shear mixing, together with some diffusive and convective mixing. As a result, the various materials in the mixing bin are thoroughly mixed by the spiral auger. After the machine has been driven into the cattle sheds, the processed uniform rations are delivered directly to the cattle through the discharging device according to the individual needs of the cattle, and the whole task from processing to supply of rations is thus completed.


## 2.1. Structure

The self-propelled TMR mixer consists of a diesel engine unit, a loading device, a conveyor, a mixing device, and a discharging device. A two-dimensional sketch is shown in Figure1. The power for the whole machine is provided by the diesel engine, which drives a hydraulic pump. By operating a valve, the direction of the hydraulic oil flow can be changed, thus changing the direction of motion of the motors driving the road wheels and allowing the whole machine to move forward or backward. The loading device consists principally of a roller and cutters. During the loading operation in the silage pit, rotary cutters bolted to a spiral blade shift the cut silage material to the center of the loading device, and the material is quickly thrown onto the conveyor to prevent blockages. The mixing device consists principally of a vertical spiral auger and a bin and is used to mix the crude fodder and concentrate feeds.Figure 1
Structure of self-propelled TMR mixer: 1, land wheel; 2, roller; 3, cutter; 4, conveyor; 5, lifting hydraulic cylinder; 6, discharging hydraulic cylinder; 7, moving blade; 8, auger; 9, bin; 10, diesel engine unit; 11, discharging device; 12, road wheels; 13, cab.The self-propelled TMR mixer was designed to realize integrated feeding in the pasture with high efficiency. Self-propelled TMR mixers like the one considered here have played an important role in speeding up the utilization of crop straw as feed. This constitutes one component of the implementation of a policy of “grain to feed” with the aim of alleviating the problems arising when people and livestock compete for grain.


## 2.2. Operating Procedure

The self-propelled TMR mixer realizes the combined functions of loading and conveying of material in the silage pit, mixing of roughage and concentrate, and transporting and discharging of rations. 3D rendering of the whole machine is shown in Figure2.Figure 2
3D rendering of the self-propelled TMR mixer.When the machine is working, a variety of feeds are loaded, depending on the dietary formulation. The lifting action of the loading device on the material in the silage pit and the motion of the conveyor are realized by the reciprocating motion of the piston rod of the double-acting lifting hydraulic cylinder (Figure1 (5)). The cutting of the material in the silage pit is carried out by the rotating loading cutter roller (Figure 1 (2)), and the cut material is then transported to the bin (Figure 1 (9)) of the mixing device by the conveyor (Figure 1 (4)). Driven by a hydraulic motor, a spiral blade radially distributed along the auger (Figure 1 (8)) rotates with the central shaft, and the material moves along the circumferential and axial directions of the spiral. When this material has piled up to a certain height, it falls away freely. As the spiral blade passes through the material, the latter is subjected to inertial and frictional forces, which can cause three-dimensional mixing in the axial, circumferential, and radial directions. This mixing process is mainly shear mixing, together with some diffusive and convective mixing. As a result, the various materials in the mixing bin are thoroughly mixed by the spiral auger. After the machine has been driven into the cattle sheds, the processed uniform rations are delivered directly to the cattle through the discharging device according to the individual needs of the cattle, and the whole task from processing to supply of rations is thus completed.


## 3. Materials and Methods

### 3.1. Simulation of Silage Cutting

#### 3.1.1. SPH-FEM Coupling Model

The coupling diagram of SPH particles and FEM meshes is shown in Figure3. On the left are the SPH particles in the large-deformation area and on the right are the FEM meshes in the small-deformation area [30, 31]. The SPH particles at the coupling interface are in one-to-one correspondence with the FEM mesh nodes, and transmission of mechanical parameters is accomplished through penalty function constraints [32]. For the contact boundary constraints on the SPH particles and FEM meshes, the coupling is defined as fixed disconnected contact of the node–surface contact type from ANSYS/LS-DYNA 971 software. The failure criterion for the fixed disconnected contact is [33, 34](1)fnfn,failm1+fsfs,failm2≥1,where fn is the normal force (N); fs is the shear force (N); fn,fail is the normal failure force (N); fs,fail is the shear failure force (N); m1 is the normal force index; and m2 is the shear force index.Figure 3
Schematic diagram of SPH particles coupling with FEM meshes.


#### 3.1.2. Simplification and Assumption of Simulation Model

(1)
The roller and the spiral blade are welded together, while the cutters and the spiral blade are bolted. In the simplified model, parts such as bolts and nuts are ignored and the material of the loading cutter roller is taken to be isotropic and linear elastic to reduce the simulation time [35–37].
(2)
It is assumed that throughout the cutting of the silage material, the rotational speed and the vertical forward speed of the loading cutter roller are constant, and the distance between the roller shaft and the silage material is defined in the same plane; that is, the cutting depth of the silage material is taken to be constant.


#### 3.1.3. Finite Element Model of the Silage and Loading Cutter Roller

A three-dimensional solid model of silage and cutter roller is established in SolidWorks software (version 2016) and imported into ANSYS for preprocessing. The silage material is modeled as a cuboid of overall dimensions 2000 mm × 150 mm × 1000 mm using the ∗MAT_173 (∗MAT_MOHR_COULOMB) material model from ANSYS/LS-DYNA 971, which is suitable for granular materials. The classical expression for the shear stress strength is [38](2)τmax=C+σntanϕ,where τmax is the maximum shear stress (MPa), C is the cohesion (MPa), σn is the normal stress (MPa), and φ is the internal friction angle (degrees).The ∗MAT_020 (∗MAT_RIGID) rigid-body material model from ANSYS/LS-DYNA 971 is used for the loading cutter roller, with the material density of the roller taken as 7830 kg/m3, the elastic modulus as 2 × 105 MPa, and Poisson’s ratio as 0.3. The design width of the cutter roller is 1600 mm, the number of cutters is 48, and the thickness of each cutter is 4 mm. The finite element model of the cutter in ANSYS software is shown in Figure 4.Figure 4
Finite element model of the cutter.According to the material properties in a local silage pit and data from the literature [39], the silage density is taken as 420 kg m−3, the friction angle as 34°, the shear modulus as 21 MPa, and the cohesion as 6.2 × 10−3 MPa, while the other parameters of the silage material are set according to the default values of ∗MAT_173 in ANSYS/LS-DYNA 971. The whole simulation model consists of three parts: the loading cutter roller, the inner silage (SPH), and the outer silage (FEM). There needs to be a space between the inner and outer silage to avoid an initial penetration problem during the simulation. Both the inner and outer silage are divided into hexahedral meshes by the sweeping method, with the mesh size for the outer silage being larger than that for the inner silage. To improve the accuracy of the simulation results, the free meshing method is used for the loading cutter roller, and local mesh refinement is carried out at the cutting edge. The element type of the cutter roller is set as 3D Solid164, and the inner silage is transformed into SPH particles by the solid nodes method. The simplified model of the silage and loading cutter roller is shown in Figure 5.Figure 5
Simplified finite element model of the silage and loading cutter roller: 1, loading cutter roller; 2, inner silage (SPH); 3, outer silage (FEM).All nodes at the bottom of the finite element model in Figure5 are defined as node group 1, all SPH particles of the inner silage are defined as node group 2, a layer of SPH particles near the outer silage is defined as node group 3, and the outer side of the outer silage is defined as surface group 1. The keyword ∗BOUNDARY_SPC_SET is added to constrain all degrees of freedom of node group 1. ∗BOUNDARY_NON_REFLECTING is added to define surface group 1 as a nonreflective boundary to eliminate the reflection effect of stress waves. ∗ERODING_NODES_TO_SURFACE is added to define the contact type between node group 2 and the loading cutter roller as node–surface erosion contact (which is used to represent failure of the solid elements). The cutter roller is set as the main contact and node group 2 as the slave contact. ∗NODES_TO_SURFACE_CONSTRAINED_OFFSET is added to define the contact type between node group 3 and the outer silage. The outer silage is set as the main contact and node group 3 as the slave contact. ∗BOUNDARY_SPH_SYMMETRY_PLANE is added to set the symmetry plane of node group 2 as YOZ. The rotational speed of the loading cutter roller is taken as 240 rev·min−1 and the vertical forward velocity as 8 m·s−1 (for the design of the self-propelled TMR mixer, the hydraulic motor BM4-390 is selected to drive the rotation of the loading cutter roller, and its working speed is 240 rev·min−1. Moreover, the radial distance between the cutter tip and the center of the roller is 320 mm. Therefore, the forward velocity of the cutters can be calculated from the formula V=2πrn/60000 to be 8 m·s−1). To ensure that the loading cutter roller can cut out the entire silage model zone during the simulation time, this is set as 200 ms. All motion constraints except translation along the Y axis and rotation around X axis are applied to the loading cutter roller.


### 3.2. Numerical Simulation of Material Mixing

#### 3.2.1. Structural Parameters of Mixing Device

The mixing device of the self-propelled TMR mixer is composed of a bin and a spiral auger; the three-dimensional model structures of which are shown in Figure6. The bin structure adopted in this paper is in the form of a truncated cone, with the volume given by(3)V=13πH×R12+R22+R1R2,where R1 and R2 are the bottom and top radii of the bin and H is the bin height. On substituting the values of 1260 mm, 1870 mm, and 2300 mm for R1, R2, and H, respectively, the volume can be obtained.Figure 6
Mixing device: (a) bin; (b) auger.
(a)
(b)The volume of the mixing bin is found to be 17.92 m3, and when the volume occupied by the spiral auger is removed (the volume of the spiral auger was calculated using the 3D design software SolidWorks), the effective mixing volume of the bin is obtained, which is about 16 m3. The angle between the side wall and bottom wall of such bins is usually in the range 105°–120° [17], and in this paper it is taken as 106°, which meets the requirement for smooth sliding of the mixed material along the walls of the bin.The mixing of material in the bin is a relatively complex process, mainly involving the three modes of shear, diffusion, and convection; the angle, diameter, and pitch of the spiral auger blade are the main structural parameters that affect the mixing of the material [11, 16]. The critical speed of a vertical spiral auger (rev·min−1) is given by(4)nc=42.28tanα+φsdtanμf,where d is the diameter of the spiral blade (m), α is the angle of the blade (degrees), φs is the friction angle between the material and the blade (which in this paper is taken as 21.8°) [40], and μf is the repose angle of the material (taken here as 55°). The angle α is given by(5)tanα=Sπd,where S is the pitch of the spiral auger blade (taken here as 440 mm).Takingd to be the minimum diameter of the spiral auger blade (768 mm), it is found that α = 10.34° and nc = 32 rev·min−1, while taking d to be the maximum blade diameter (2470 mm), α = 3.25° and nc = 15 rev·min−1.The theoretical critical speed range of the auger, which is the lowest speed at which a material particle on the spiral auger blade can be lifted, is 15–32 rev·min−1. Because the spiral auger has a conical structure, the critical speed of the material particles varies with position on the blade, and the theoretical critical speed may deviate from the actual mixing working speed. Consequently, the auger working speed for optimal mixing of material needs to be determined experimentally.


#### 3.2.2. Contact Force Model

The DEM is used to simulate the motion of groups of material particle. Individual particle motions will inevitably cause collisions and thus forces between particles [24, 27]. The contact model is an important component of the DEM, and in this paper, the Hertz–Mindlin contact model is adopted, as shown in Figure 7. This model assumes that the discrete elements are rigid bodies and that there are point contacts between elements with soft contact characteristics (i.e., a certain amount of overlap is allowed between the rigid discrete elements at contact points). The contact force is calculated according to the normal overlap and tangential displacement between particles, and the resultant force on any discrete element is determined by its interaction with the discrete elements with which it is in contact.Figure 7
Model of contact force:Fn is the normal force, Ft is the tangential force, kn and kt are the normal and tangential stiffnesses, Cn and Ct are the normal and tangential damping coefficients, and μ is the friction coefficient.The contact mechanics analysis is based on the following equations:(6)Fn=43E∗R∗δn3/2,Ft=−ktδt,Fn,d=−256βknm∗vn,rel,Ft,d=−256βktm∗vt,rel,where Fn and Ft are the normal and tangential forces (N), E∗ is the equivalent elastic modulus (Pa), R∗ is the equivalent contact radius (m), δn and δt are the normal and tangential overlaps (m), kn and kt are the normal and tangential stiffnesses (N·m−1), Fn,d and Ft,d are the normal and tangential damping forces (N), β is the damping ratio, m∗ is the equivalent mass (kg), and vn,rel and vt,rel are the normal and tangential relative velocities (m·s−1).


#### 3.2.3. Establishment of EDEM Particle Model

SolidWorks 3D software is used to model the rotating auger and the fixed bin of the mixing device. To facilitate the simulation, the components of the device that are irrelevant to material contact during the mixing motion are removed and imported into the EDEM software in the IGES format, as shown in Figure8.Figure 8
Simulation model of the mixing device.Many limitations and difficulties regarding the application of DEM for real industrial-scale applications exist. Limitations on the particle size are mainly due to computational limitations as a result of the higher number of particles to be considered resulting in a drastic increase of necessary calculations. In addition, the computational effort increases dramatically with the decrease of the particle size due to the smaller time step, which results in an increased number of necessary iterations. In the DEM simulation, the time step,ΔTstep, is defined as the time between each iteration. A simulation is stable only if the time step employed is lower than a critical time step, ΔTcritical, which is generally defined as a fraction of the natural frequency of an equivalent mass-spring system. For the commercial EDEM software, the critical time step is defined as the critical Rayleigh time step, TR, as follows [41]:(7)ΔTstep<ΔTcritical=TR=πRρG1/20.1631υ+0.8766−1,where TR is a function of the particle properties such as radius and density, R is the particle radius (m), ρ is the particle density (kg·m−3), G is the shear modulus (MPa), and v is Poisson’s ratio.The shear modulus,G, has a major influence in the critical time step. To arbitrarily reduce the value of the shear modulus as this dramatically increases the critical time step, the simulation time is reduced. However, it has been shown that this approach could lead to errors in the simulation results [42]. Therefore, three similarity principles, namely, geometric, mechanical, and dynamic, under which the numerical model can exactly reproduce the mechanical behavior of a physical model, are established. All the conditions that ensure an exact equivalence between the physical and scaled models in general dynamic cases can be summarized as [43](8)RmRp=hρmρp=h−2umup=hFmum,RmFpup,Rp=QmtQpt=h2,where h is the scale factor, Rm and Rp are the radii of an arbitrary particle in both scaled and physical models (m), ρm and ρp are the densities of scaled and physical models (kg·m−3), um and up are the displacements of scaled and physical models (m), Fm (um, Rm) and Fp (up, Rp) are the interaction forces from neighboring particles of scaled and physical models (N), and Qm (t) and Qp (t) are any external forces applied to the particle of scaled and physical models (N).It is well known that there are three basic quantities for mechanical systems in the SI system: length (L), mass (M), and time (T). After the three basic quantities are chosen, their corresponding conversion or scale factors between the physical and the scaled models can be arbitrarily set in principle.Density as a basic parameter of materials that can be directly obtained, and mass density (ρ) replaces mass as a basic quantity, but length and time remain as the other two basic quantities. The corresponding scale factors for these three quantities, λL, λT, and λρ are specifically chosen as [44](9)λL=h1,λT=h,λρ=1,i.e., the length scale factor is the same as the (spatial) scale factor, <i>h</i>, as expected, and the temporal (or time) scale is the same as the spatial factor, while the density is kept the same for both systems.Because the volume of the bin studied in this paper is large (16 m3) and the size of the material particles is relatively small, in order to avoid the “false crash” state of the CPU [25, 45], therefore, based on similarity principles mentioned above, in the simulation model, the sizes of the bin and the auger are both reduced to 0.25 of their sizes in a realistic model; similarly, the motion timescale of the material particles in the bin is also scaled to 0.25 of the original.On the basis of data in the literature [11, 13], a ratio of concentrate to roughage of 35 : 65 (based on the ratio of dry matter mass) is used to determine the compositions of the experimental rations. The total mass of the mixture produced by the EDEM particle factory is set as 1000 kg and the mass ratio of crop straw to globular corn grain to silage block as 10 : 35 : 55. All the particles in the mixture are assumed to be spherical, and from measurements of an actual mixture, the average particle radii for the three components are obtained as 2 mm, 4 mm, and 3 mm for crop straw, globular corn grain, and silage block, respectively. The Hertz–Mindlin contact model is used for the interaction between material particles and mixing device, and each crop particle model is defined as rigid particles in this paper. The three types of material particles in the discrete element model are shown in Figures 9(a)–9(c). To clearly show the motion of particles in the bin, the crop straw, globular corn grain, and silage block particles are colored yellow, red, and green, respectively, in EDEM, as shown in Figure 9(d).Figure 9
Discrete element model of material particles: (a) crop straw; (b) globular corn grain; (c) silage block; (d) partial view of material at a given time.
(a)
(b)
(c)
(d)The accuracy of DEM relies on the rationality of simulation parameters which should be calibrated before the simulation [41]. The simulation test was based on the results of physical test for the angle of repose, and the discrete element parameters were calibrated by measuring the angle of repose [46, 47]. The GEMM (generic EDEM material model) database contains thousands of material models representing a wide range of materials such as rocks, soils, and ores. It is particularly suited for engineers in the mining, construction, and agriculture industries who design such heavy equipment handling materials. The key contact property parameters including coefficient of restitution, coefficient of static friction, and coefficient of rolling friction can be obtained from the GEMM database [48, 49]. To increase the efficiency of particle generation, three different virtual factories of the static BOX type are created for the three types of material particles. After a start-up time of spiral auger rotation of 0.5 s, the total simulation time is set as 100 s. Therefore, in this paper, according to the GEMM database and data in the literature [50, 51], the material mechanical properties and the material contact property parameters are shown in Tables 1 and 2, respectively.Table 1
Material mechanical properties.
Material
Poisson’s ratio
Shear modulus (Pa)
Density (kg m−3)
Crop straw
0.30
6 × 108
57
Globular corn grain
0.40
1.37 × 108
1256
Silage block
0.30
2.1 × 107
420
Mixing device
0.29
7.7 × 1010
7850Table 2
Material contact property parameters used in the simulation.
Interaction material
Restitution
Static friction
Rolling friction
Crop straw-crop straw
0.33
0.38
0.01
Crop straw-globular corn grain
0.44
0.46
0.11
Crop straw-silage block
0.30
0.76
0.17
Crop straw-mixing device
0.34
0.42
0.01
Globular corn grain-globular corn grain
0.25
0.70
0.05
Globular corn grain-silage block
0.39
0.54
0.15
Globular corn grain-mixing device
0.30
0.50
0.07
Silage block-silage block
0.37
1.29
0.21
Silage block-mixing device
0.41
1.08
0.13


#### 3.2.4. Factors and Evaluation Index

The uniformity of mixing of the material particles is the most important index for assessing the performance of the mixing device [1, 7, 15]. In this paper, the grid method is used to establish a model for evaluating this index. The simulation area is divided into 6 × 6 × 6 = 216 grids by the binning method in the EDEM software, as shown in Figure 10.Figure 10
Simulation model of the grid bin.If the grid size of the simulation model is too small, there will be an insufficient number of particles in a single grid to provide an accurate reflection of their dispersion, while if the grid size is too large, the results will not have statistical significance. Therefore, each grid cell should contain a sufficient number of material particles to ensure the validity of the grid cell data statistics, and it is necessary to perform a preliminary screening of the defined grid and remove those grid cells with less than 1000 material particles. The ratio of sample particles to the total number of particles in each grid cell is the sampling ratio, and the ratio of sample particles to the total number of particles in all materials is the overall ratio. The standard deviation of the sample is calculated as(10)S=∑i=1nXi−X¯2n−1,and the mixing uniformity of material particles is calculated as(11)M=1−SX¯×100%,where n is the number of grid cells, Xi is the sampling ratio of the sample particles, and X¯ is the overall ratio of the sample particles.According to the literature [11, 13, 14, 16], the factors affecting the performance of the mixing device include the material mixing time, the loading rate, and the auger speed. From the calculation of the critical speed of the auger described above and from practical experience with ration mixers, the following ranges of values of these factors are adopted in the simulations: a mixing time of 30–150 s, a loading rate of 50–90%, and an auger speed of 8–40 rev·min−1. The three-factor, five-level simulations are designed to analyze the mixing of material particles, with the uniformity of mixing of globular corn grains being taken as the index for evaluating mixing performance. The three factors and five levels are shown in Table 3.Table 3
Factors and levels of mixing simulation.
Level
Factor
Mixing time (s)
Loading rate (%)
Auger speed (rev·min−1)
1
30
50
8
2
60
60
16
3
90
70
24
4
120
80
32
5
150
90
40


## 3.1. Simulation of Silage Cutting

### 3.1.1. SPH-FEM Coupling Model

The coupling diagram of SPH particles and FEM meshes is shown in Figure3. On the left are the SPH particles in the large-deformation area and on the right are the FEM meshes in the small-deformation area [30, 31]. The SPH particles at the coupling interface are in one-to-one correspondence with the FEM mesh nodes, and transmission of mechanical parameters is accomplished through penalty function constraints [32]. For the contact boundary constraints on the SPH particles and FEM meshes, the coupling is defined as fixed disconnected contact of the node–surface contact type from ANSYS/LS-DYNA 971 software. The failure criterion for the fixed disconnected contact is [33, 34](1)fnfn,failm1+fsfs,failm2≥1,where fn is the normal force (N); fs is the shear force (N); fn,fail is the normal failure force (N); fs,fail is the shear failure force (N); m1 is the normal force index; and m2 is the shear force index.Figure 3
Schematic diagram of SPH particles coupling with FEM meshes.


### 3.1.2. Simplification and Assumption of Simulation Model

(1)
The roller and the spiral blade are welded together, while the cutters and the spiral blade are bolted. In the simplified model, parts such as bolts and nuts are ignored and the material of the loading cutter roller is taken to be isotropic and linear elastic to reduce the simulation time [35–37].
(2)
It is assumed that throughout the cutting of the silage material, the rotational speed and the vertical forward speed of the loading cutter roller are constant, and the distance between the roller shaft and the silage material is defined in the same plane; that is, the cutting depth of the silage material is taken to be constant.


### 3.1.3. Finite Element Model of the Silage and Loading Cutter Roller

A three-dimensional solid model of silage and cutter roller is established in SolidWorks software (version 2016) and imported into ANSYS for preprocessing. The silage material is modeled as a cuboid of overall dimensions 2000 mm × 150 mm × 1000 mm using the ∗MAT_173 (∗MAT_MOHR_COULOMB) material model from ANSYS/LS-DYNA 971, which is suitable for granular materials. The classical expression for the shear stress strength is [38](2)τmax=C+σntanϕ,where τmax is the maximum shear stress (MPa), C is the cohesion (MPa), σn is the normal stress (MPa), and φ is the internal friction angle (degrees).The ∗MAT_020 (∗MAT_RIGID) rigid-body material model from ANSYS/LS-DYNA 971 is used for the loading cutter roller, with the material density of the roller taken as 7830 kg/m3, the elastic modulus as 2 × 105 MPa, and Poisson’s ratio as 0.3. The design width of the cutter roller is 1600 mm, the number of cutters is 48, and the thickness of each cutter is 4 mm. The finite element model of the cutter in ANSYS software is shown in Figure 4.Figure 4
Finite element model of the cutter.According to the material properties in a local silage pit and data from the literature [39], the silage density is taken as 420 kg m−3, the friction angle as 34°, the shear modulus as 21 MPa, and the cohesion as 6.2 × 10−3 MPa, while the other parameters of the silage material are set according to the default values of ∗MAT_173 in ANSYS/LS-DYNA 971. The whole simulation model consists of three parts: the loading cutter roller, the inner silage (SPH), and the outer silage (FEM). There needs to be a space between the inner and outer silage to avoid an initial penetration problem during the simulation. Both the inner and outer silage are divided into hexahedral meshes by the sweeping method, with the mesh size for the outer silage being larger than that for the inner silage. To improve the accuracy of the simulation results, the free meshing method is used for the loading cutter roller, and local mesh refinement is carried out at the cutting edge. The element type of the cutter roller is set as 3D Solid164, and the inner silage is transformed into SPH particles by the solid nodes method. The simplified model of the silage and loading cutter roller is shown in Figure 5.Figure 5
Simplified finite element model of the silage and loading cutter roller: 1, loading cutter roller; 2, inner silage (SPH); 3, outer silage (FEM).All nodes at the bottom of the finite element model in Figure5 are defined as node group 1, all SPH particles of the inner silage are defined as node group 2, a layer of SPH particles near the outer silage is defined as node group 3, and the outer side of the outer silage is defined as surface group 1. The keyword ∗BOUNDARY_SPC_SET is added to constrain all degrees of freedom of node group 1. ∗BOUNDARY_NON_REFLECTING is added to define surface group 1 as a nonreflective boundary to eliminate the reflection effect of stress waves. ∗ERODING_NODES_TO_SURFACE is added to define the contact type between node group 2 and the loading cutter roller as node–surface erosion contact (which is used to represent failure of the solid elements). The cutter roller is set as the main contact and node group 2 as the slave contact. ∗NODES_TO_SURFACE_CONSTRAINED_OFFSET is added to define the contact type between node group 3 and the outer silage. The outer silage is set as the main contact and node group 3 as the slave contact. ∗BOUNDARY_SPH_SYMMETRY_PLANE is added to set the symmetry plane of node group 2 as YOZ. The rotational speed of the loading cutter roller is taken as 240 rev·min−1 and the vertical forward velocity as 8 m·s−1 (for the design of the self-propelled TMR mixer, the hydraulic motor BM4-390 is selected to drive the rotation of the loading cutter roller, and its working speed is 240 rev·min−1. Moreover, the radial distance between the cutter tip and the center of the roller is 320 mm. Therefore, the forward velocity of the cutters can be calculated from the formula V=2πrn/60000 to be 8 m·s−1). To ensure that the loading cutter roller can cut out the entire silage model zone during the simulation time, this is set as 200 ms. All motion constraints except translation along the Y axis and rotation around X axis are applied to the loading cutter roller.


## 3.1.1. SPH-FEM Coupling Model

The coupling diagram of SPH particles and FEM meshes is shown in Figure3. On the left are the SPH particles in the large-deformation area and on the right are the FEM meshes in the small-deformation area [30, 31]. The SPH particles at the coupling interface are in one-to-one correspondence with the FEM mesh nodes, and transmission of mechanical parameters is accomplished through penalty function constraints [32]. For the contact boundary constraints on the SPH particles and FEM meshes, the coupling is defined as fixed disconnected contact of the node–surface contact type from ANSYS/LS-DYNA 971 software. The failure criterion for the fixed disconnected contact is [33, 34](1)fnfn,failm1+fsfs,failm2≥1,where fn is the normal force (N); fs is the shear force (N); fn,fail is the normal failure force (N); fs,fail is the shear failure force (N); m1 is the normal force index; and m2 is the shear force index.Figure 3
Schematic diagram of SPH particles coupling with FEM meshes.


## 3.1.2. Simplification and Assumption of Simulation Model

(1)
The roller and the spiral blade are welded together, while the cutters and the spiral blade are bolted. In the simplified model, parts such as bolts and nuts are ignored and the material of the loading cutter roller is taken to be isotropic and linear elastic to reduce the simulation time [35–37].
(2)
It is assumed that throughout the cutting of the silage material, the rotational speed and the vertical forward speed of the loading cutter roller are constant, and the distance between the roller shaft and the silage material is defined in the same plane; that is, the cutting depth of the silage material is taken to be constant.


## 3.1.3. Finite Element Model of the Silage and Loading Cutter Roller

A three-dimensional solid model of silage and cutter roller is established in SolidWorks software (version 2016) and imported into ANSYS for preprocessing. The silage material is modeled as a cuboid of overall dimensions 2000 mm × 150 mm × 1000 mm using the ∗MAT_173 (∗MAT_MOHR_COULOMB) material model from ANSYS/LS-DYNA 971, which is suitable for granular materials. The classical expression for the shear stress strength is [38](2)τmax=C+σntanϕ,where τmax is the maximum shear stress (MPa), C is the cohesion (MPa), σn is the normal stress (MPa), and φ is the internal friction angle (degrees).The ∗MAT_020 (∗MAT_RIGID) rigid-body material model from ANSYS/LS-DYNA 971 is used for the loading cutter roller, with the material density of the roller taken as 7830 kg/m3, the elastic modulus as 2 × 105 MPa, and Poisson’s ratio as 0.3. The design width of the cutter roller is 1600 mm, the number of cutters is 48, and the thickness of each cutter is 4 mm. The finite element model of the cutter in ANSYS software is shown in Figure 4.Figure 4
Finite element model of the cutter.According to the material properties in a local silage pit and data from the literature [39], the silage density is taken as 420 kg m−3, the friction angle as 34°, the shear modulus as 21 MPa, and the cohesion as 6.2 × 10−3 MPa, while the other parameters of the silage material are set according to the default values of ∗MAT_173 in ANSYS/LS-DYNA 971. The whole simulation model consists of three parts: the loading cutter roller, the inner silage (SPH), and the outer silage (FEM). There needs to be a space between the inner and outer silage to avoid an initial penetration problem during the simulation. Both the inner and outer silage are divided into hexahedral meshes by the sweeping method, with the mesh size for the outer silage being larger than that for the inner silage. To improve the accuracy of the simulation results, the free meshing method is used for the loading cutter roller, and local mesh refinement is carried out at the cutting edge. The element type of the cutter roller is set as 3D Solid164, and the inner silage is transformed into SPH particles by the solid nodes method. The simplified model of the silage and loading cutter roller is shown in Figure 5.Figure 5
Simplified finite element model of the silage and loading cutter roller: 1, loading cutter roller; 2, inner silage (SPH); 3, outer silage (FEM).All nodes at the bottom of the finite element model in Figure5 are defined as node group 1, all SPH particles of the inner silage are defined as node group 2, a layer of SPH particles near the outer silage is defined as node group 3, and the outer side of the outer silage is defined as surface group 1. The keyword ∗BOUNDARY_SPC_SET is added to constrain all degrees of freedom of node group 1. ∗BOUNDARY_NON_REFLECTING is added to define surface group 1 as a nonreflective boundary to eliminate the reflection effect of stress waves. ∗ERODING_NODES_TO_SURFACE is added to define the contact type between node group 2 and the loading cutter roller as node–surface erosion contact (which is used to represent failure of the solid elements). The cutter roller is set as the main contact and node group 2 as the slave contact. ∗NODES_TO_SURFACE_CONSTRAINED_OFFSET is added to define the contact type between node group 3 and the outer silage. The outer silage is set as the main contact and node group 3 as the slave contact. ∗BOUNDARY_SPH_SYMMETRY_PLANE is added to set the symmetry plane of node group 2 as YOZ. The rotational speed of the loading cutter roller is taken as 240 rev·min−1 and the vertical forward velocity as 8 m·s−1 (for the design of the self-propelled TMR mixer, the hydraulic motor BM4-390 is selected to drive the rotation of the loading cutter roller, and its working speed is 240 rev·min−1. Moreover, the radial distance between the cutter tip and the center of the roller is 320 mm. Therefore, the forward velocity of the cutters can be calculated from the formula V=2πrn/60000 to be 8 m·s−1). To ensure that the loading cutter roller can cut out the entire silage model zone during the simulation time, this is set as 200 ms. All motion constraints except translation along the Y axis and rotation around X axis are applied to the loading cutter roller.


## 3.2. Numerical Simulation of Material Mixing

### 3.2.1. Structural Parameters of Mixing Device

The mixing device of the self-propelled TMR mixer is composed of a bin and a spiral auger; the three-dimensional model structures of which are shown in Figure6. The bin structure adopted in this paper is in the form of a truncated cone, with the volume given by(3)V=13πH×R12+R22+R1R2,where R1 and R2 are the bottom and top radii of the bin and H is the bin height. On substituting the values of 1260 mm, 1870 mm, and 2300 mm for R1, R2, and H, respectively, the volume can be obtained.Figure 6
Mixing device: (a) bin; (b) auger.
(a)
(b)The volume of the mixing bin is found to be 17.92 m3, and when the volume occupied by the spiral auger is removed (the volume of the spiral auger was calculated using the 3D design software SolidWorks), the effective mixing volume of the bin is obtained, which is about 16 m3. The angle between the side wall and bottom wall of such bins is usually in the range 105°–120° [17], and in this paper it is taken as 106°, which meets the requirement for smooth sliding of the mixed material along the walls of the bin.The mixing of material in the bin is a relatively complex process, mainly involving the three modes of shear, diffusion, and convection; the angle, diameter, and pitch of the spiral auger blade are the main structural parameters that affect the mixing of the material [11, 16]. The critical speed of a vertical spiral auger (rev·min−1) is given by(4)nc=42.28tanα+φsdtanμf,where d is the diameter of the spiral blade (m), α is the angle of the blade (degrees), φs is the friction angle between the material and the blade (which in this paper is taken as 21.8°) [40], and μf is the repose angle of the material (taken here as 55°). The angle α is given by(5)tanα=Sπd,where S is the pitch of the spiral auger blade (taken here as 440 mm).Takingd to be the minimum diameter of the spiral auger blade (768 mm), it is found that α = 10.34° and nc = 32 rev·min−1, while taking d to be the maximum blade diameter (2470 mm), α = 3.25° and nc = 15 rev·min−1.The theoretical critical speed range of the auger, which is the lowest speed at which a material particle on the spiral auger blade can be lifted, is 15–32 rev·min−1. Because the spiral auger has a conical structure, the critical speed of the material particles varies with position on the blade, and the theoretical critical speed may deviate from the actual mixing working speed. Consequently, the auger working speed for optimal mixing of material needs to be determined experimentally.


### 3.2.2. Contact Force Model

The DEM is used to simulate the motion of groups of material particle. Individual particle motions will inevitably cause collisions and thus forces between particles [24, 27]. The contact model is an important component of the DEM, and in this paper, the Hertz–Mindlin contact model is adopted, as shown in Figure 7. This model assumes that the discrete elements are rigid bodies and that there are point contacts between elements with soft contact characteristics (i.e., a certain amount of overlap is allowed between the rigid discrete elements at contact points). The contact force is calculated according to the normal overlap and tangential displacement between particles, and the resultant force on any discrete element is determined by its interaction with the discrete elements with which it is in contact.Figure 7
Model of contact force:Fn is the normal force, Ft is the tangential force, kn and kt are the normal and tangential stiffnesses, Cn and Ct are the normal and tangential damping coefficients, and μ is the friction coefficient.The contact mechanics analysis is based on the following equations:(6)Fn=43E∗R∗δn3/2,Ft=−ktδt,Fn,d=−256βknm∗vn,rel,Ft,d=−256βktm∗vt,rel,where Fn and Ft are the normal and tangential forces (N), E∗ is the equivalent elastic modulus (Pa), R∗ is the equivalent contact radius (m), δn and δt are the normal and tangential overlaps (m), kn and kt are the normal and tangential stiffnesses (N·m−1), Fn,d and Ft,d are the normal and tangential damping forces (N), β is the damping ratio, m∗ is the equivalent mass (kg), and vn,rel and vt,rel are the normal and tangential relative velocities (m·s−1).


### 3.2.3. Establishment of EDEM Particle Model

SolidWorks 3D software is used to model the rotating auger and the fixed bin of the mixing device. To facilitate the simulation, the components of the device that are irrelevant to material contact during the mixing motion are removed and imported into the EDEM software in the IGES format, as shown in Figure8.Figure 8
Simulation model of the mixing device.Many limitations and difficulties regarding the application of DEM for real industrial-scale applications exist. Limitations on the particle size are mainly due to computational limitations as a result of the higher number of particles to be considered resulting in a drastic increase of necessary calculations. In addition, the computational effort increases dramatically with the decrease of the particle size due to the smaller time step, which results in an increased number of necessary iterations. In the DEM simulation, the time step,ΔTstep, is defined as the time between each iteration. A simulation is stable only if the time step employed is lower than a critical time step, ΔTcritical, which is generally defined as a fraction of the natural frequency of an equivalent mass-spring system. For the commercial EDEM software, the critical time step is defined as the critical Rayleigh time step, TR, as follows [41]:(7)ΔTstep<ΔTcritical=TR=πRρG1/20.1631υ+0.8766−1,where TR is a function of the particle properties such as radius and density, R is the particle radius (m), ρ is the particle density (kg·m−3), G is the shear modulus (MPa), and v is Poisson’s ratio.The shear modulus,G, has a major influence in the critical time step. To arbitrarily reduce the value of the shear modulus as this dramatically increases the critical time step, the simulation time is reduced. However, it has been shown that this approach could lead to errors in the simulation results [42]. Therefore, three similarity principles, namely, geometric, mechanical, and dynamic, under which the numerical model can exactly reproduce the mechanical behavior of a physical model, are established. All the conditions that ensure an exact equivalence between the physical and scaled models in general dynamic cases can be summarized as [43](8)RmRp=hρmρp=h−2umup=hFmum,RmFpup,Rp=QmtQpt=h2,where h is the scale factor, Rm and Rp are the radii of an arbitrary particle in both scaled and physical models (m), ρm and ρp are the densities of scaled and physical models (kg·m−3), um and up are the displacements of scaled and physical models (m), Fm (um, Rm) and Fp (up, Rp) are the interaction forces from neighboring particles of scaled and physical models (N), and Qm (t) and Qp (t) are any external forces applied to the particle of scaled and physical models (N).It is well known that there are three basic quantities for mechanical systems in the SI system: length (L), mass (M), and time (T). After the three basic quantities are chosen, their corresponding conversion or scale factors between the physical and the scaled models can be arbitrarily set in principle.Density as a basic parameter of materials that can be directly obtained, and mass density (ρ) replaces mass as a basic quantity, but length and time remain as the other two basic quantities. The corresponding scale factors for these three quantities, λL, λT, and λρ are specifically chosen as [44](9)λL=h1,λT=h,λρ=1,i.e., the length scale factor is the same as the (spatial) scale factor, <i>h</i>, as expected, and the temporal (or time) scale is the same as the spatial factor, while the density is kept the same for both systems.Because the volume of the bin studied in this paper is large (16 m3) and the size of the material particles is relatively small, in order to avoid the “false crash” state of the CPU [25, 45], therefore, based on similarity principles mentioned above, in the simulation model, the sizes of the bin and the auger are both reduced to 0.25 of their sizes in a realistic model; similarly, the motion timescale of the material particles in the bin is also scaled to 0.25 of the original.On the basis of data in the literature [11, 13], a ratio of concentrate to roughage of 35 : 65 (based on the ratio of dry matter mass) is used to determine the compositions of the experimental rations. The total mass of the mixture produced by the EDEM particle factory is set as 1000 kg and the mass ratio of crop straw to globular corn grain to silage block as 10 : 35 : 55. All the particles in the mixture are assumed to be spherical, and from measurements of an actual mixture, the average particle radii for the three components are obtained as 2 mm, 4 mm, and 3 mm for crop straw, globular corn grain, and silage block, respectively. The Hertz–Mindlin contact model is used for the interaction between material particles and mixing device, and each crop particle model is defined as rigid particles in this paper. The three types of material particles in the discrete element model are shown in Figures 9(a)–9(c). To clearly show the motion of particles in the bin, the crop straw, globular corn grain, and silage block particles are colored yellow, red, and green, respectively, in EDEM, as shown in Figure 9(d).Figure 9
Discrete element model of material particles: (a) crop straw; (b) globular corn grain; (c) silage block; (d) partial view of material at a given time.
(a)
(b)
(c)
(d)The accuracy of DEM relies on the rationality of simulation parameters which should be calibrated before the simulation [41]. The simulation test was based on the results of physical test for the angle of repose, and the discrete element parameters were calibrated by measuring the angle of repose [46, 47]. The GEMM (generic EDEM material model) database contains thousands of material models representing a wide range of materials such as rocks, soils, and ores. It is particularly suited for engineers in the mining, construction, and agriculture industries who design such heavy equipment handling materials. The key contact property parameters including coefficient of restitution, coefficient of static friction, and coefficient of rolling friction can be obtained from the GEMM database [48, 49]. To increase the efficiency of particle generation, three different virtual factories of the static BOX type are created for the three types of material particles. After a start-up time of spiral auger rotation of 0.5 s, the total simulation time is set as 100 s. Therefore, in this paper, according to the GEMM database and data in the literature [50, 51], the material mechanical properties and the material contact property parameters are shown in Tables 1 and 2, respectively.Table 1
Material mechanical properties.
Material
Poisson’s ratio
Shear modulus (Pa)
Density (kg m−3)
Crop straw
0.30
6 × 108
57
Globular corn grain
0.40
1.37 × 108
1256
Silage block
0.30
2.1 × 107
420
Mixing device
0.29
7.7 × 1010
7850Table 2
Material contact property parameters used in the simulation.
Interaction material
Restitution
Static friction
Rolling friction
Crop straw-crop straw
0.33
0.38
0.01
Crop straw-globular corn grain
0.44
0.46
0.11
Crop straw-silage block
0.30
0.76
0.17
Crop straw-mixing device
0.34
0.42
0.01
Globular corn grain-globular corn grain
0.25
0.70
0.05
Globular corn grain-silage block
0.39
0.54
0.15
Globular corn grain-mixing device
0.30
0.50
0.07
Silage block-silage block
0.37
1.29
0.21
Silage block-mixing device
0.41
1.08
0.13


### 3.2.4. Factors and Evaluation Index

The uniformity of mixing of the material particles is the most important index for assessing the performance of the mixing device [1, 7, 15]. In this paper, the grid method is used to establish a model for evaluating this index. The simulation area is divided into 6 × 6 × 6 = 216 grids by the binning method in the EDEM software, as shown in Figure 10.Figure 10
Simulation model of the grid bin.If the grid size of the simulation model is too small, there will be an insufficient number of particles in a single grid to provide an accurate reflection of their dispersion, while if the grid size is too large, the results will not have statistical significance. Therefore, each grid cell should contain a sufficient number of material particles to ensure the validity of the grid cell data statistics, and it is necessary to perform a preliminary screening of the defined grid and remove those grid cells with less than 1000 material particles. The ratio of sample particles to the total number of particles in each grid cell is the sampling ratio, and the ratio of sample particles to the total number of particles in all materials is the overall ratio. The standard deviation of the sample is calculated as(10)S=∑i=1nXi−X¯2n−1,and the mixing uniformity of material particles is calculated as(11)M=1−SX¯×100%,where n is the number of grid cells, Xi is the sampling ratio of the sample particles, and X¯ is the overall ratio of the sample particles.According to the literature [11, 13, 14, 16], the factors affecting the performance of the mixing device include the material mixing time, the loading rate, and the auger speed. From the calculation of the critical speed of the auger described above and from practical experience with ration mixers, the following ranges of values of these factors are adopted in the simulations: a mixing time of 30–150 s, a loading rate of 50–90%, and an auger speed of 8–40 rev·min−1. The three-factor, five-level simulations are designed to analyze the mixing of material particles, with the uniformity of mixing of globular corn grains being taken as the index for evaluating mixing performance. The three factors and five levels are shown in Table 3.Table 3
Factors and levels of mixing simulation.
Level
Factor
Mixing time (s)
Loading rate (%)
Auger speed (rev·min−1)
1
30
50
8
2
60
60
16
3
90
70
24
4
120
80
32
5
150
90
40


## 3.2.1. Structural Parameters of Mixing Device

The mixing device of the self-propelled TMR mixer is composed of a bin and a spiral auger; the three-dimensional model structures of which are shown in Figure6. The bin structure adopted in this paper is in the form of a truncated cone, with the volume given by(3)V=13πH×R12+R22+R1R2,where R1 and R2 are the bottom and top radii of the bin and H is the bin height. On substituting the values of 1260 mm, 1870 mm, and 2300 mm for R1, R2, and H, respectively, the volume can be obtained.Figure 6
Mixing device: (a) bin; (b) auger.
(a)
(b)The volume of the mixing bin is found to be 17.92 m3, and when the volume occupied by the spiral auger is removed (the volume of the spiral auger was calculated using the 3D design software SolidWorks), the effective mixing volume of the bin is obtained, which is about 16 m3. The angle between the side wall and bottom wall of such bins is usually in the range 105°–120° [17], and in this paper it is taken as 106°, which meets the requirement for smooth sliding of the mixed material along the walls of the bin.The mixing of material in the bin is a relatively complex process, mainly involving the three modes of shear, diffusion, and convection; the angle, diameter, and pitch of the spiral auger blade are the main structural parameters that affect the mixing of the material [11, 16]. The critical speed of a vertical spiral auger (rev·min−1) is given by(4)nc=42.28tanα+φsdtanμf,where d is the diameter of the spiral blade (m), α is the angle of the blade (degrees), φs is the friction angle between the material and the blade (which in this paper is taken as 21.8°) [40], and μf is the repose angle of the material (taken here as 55°). The angle α is given by(5)tanα=Sπd,where S is the pitch of the spiral auger blade (taken here as 440 mm).Takingd to be the minimum diameter of the spiral auger blade (768 mm), it is found that α = 10.34° and nc = 32 rev·min−1, while taking d to be the maximum blade diameter (2470 mm), α = 3.25° and nc = 15 rev·min−1.The theoretical critical speed range of the auger, which is the lowest speed at which a material particle on the spiral auger blade can be lifted, is 15–32 rev·min−1. Because the spiral auger has a conical structure, the critical speed of the material particles varies with position on the blade, and the theoretical critical speed may deviate from the actual mixing working speed. Consequently, the auger working speed for optimal mixing of material needs to be determined experimentally.


## 3.2.2. Contact Force Model

The DEM is used to simulate the motion of groups of material particle. Individual particle motions will inevitably cause collisions and thus forces between particles [24, 27]. The contact model is an important component of the DEM, and in this paper, the Hertz–Mindlin contact model is adopted, as shown in Figure 7. This model assumes that the discrete elements are rigid bodies and that there are point contacts between elements with soft contact characteristics (i.e., a certain amount of overlap is allowed between the rigid discrete elements at contact points). The contact force is calculated according to the normal overlap and tangential displacement between particles, and the resultant force on any discrete element is determined by its interaction with the discrete elements with which it is in contact.Figure 7
Model of contact force:Fn is the normal force, Ft is the tangential force, kn and kt are the normal and tangential stiffnesses, Cn and Ct are the normal and tangential damping coefficients, and μ is the friction coefficient.The contact mechanics analysis is based on the following equations:(6)Fn=43E∗R∗δn3/2,Ft=−ktδt,Fn,d=−256βknm∗vn,rel,Ft,d=−256βktm∗vt,rel,where Fn and Ft are the normal and tangential forces (N), E∗ is the equivalent elastic modulus (Pa), R∗ is the equivalent contact radius (m), δn and δt are the normal and tangential overlaps (m), kn and kt are the normal and tangential stiffnesses (N·m−1), Fn,d and Ft,d are the normal and tangential damping forces (N), β is the damping ratio, m∗ is the equivalent mass (kg), and vn,rel and vt,rel are the normal and tangential relative velocities (m·s−1).


## 3.2.3. Establishment of EDEM Particle Model

SolidWorks 3D software is used to model the rotating auger and the fixed bin of the mixing device. To facilitate the simulation, the components of the device that are irrelevant to material contact during the mixing motion are removed and imported into the EDEM software in the IGES format, as shown in Figure8.Figure 8
Simulation model of the mixing device.Many limitations and difficulties regarding the application of DEM for real industrial-scale applications exist. Limitations on the particle size are mainly due to computational limitations as a result of the higher number of particles to be considered resulting in a drastic increase of necessary calculations. In addition, the computational effort increases dramatically with the decrease of the particle size due to the smaller time step, which results in an increased number of necessary iterations. In the DEM simulation, the time step,ΔTstep, is defined as the time between each iteration. A simulation is stable only if the time step employed is lower than a critical time step, ΔTcritical, which is generally defined as a fraction of the natural frequency of an equivalent mass-spring system. For the commercial EDEM software, the critical time step is defined as the critical Rayleigh time step, TR, as follows [41]:(7)ΔTstep<ΔTcritical=TR=πRρG1/20.1631υ+0.8766−1,where TR is a function of the particle properties such as radius and density, R is the particle radius (m), ρ is the particle density (kg·m−3), G is the shear modulus (MPa), and v is Poisson’s ratio.The shear modulus,G, has a major influence in the critical time step. To arbitrarily reduce the value of the shear modulus as this dramatically increases the critical time step, the simulation time is reduced. However, it has been shown that this approach could lead to errors in the simulation results [42]. Therefore, three similarity principles, namely, geometric, mechanical, and dynamic, under which the numerical model can exactly reproduce the mechanical behavior of a physical model, are established. All the conditions that ensure an exact equivalence between the physical and scaled models in general dynamic cases can be summarized as [43](8)RmRp=hρmρp=h−2umup=hFmum,RmFpup,Rp=QmtQpt=h2,where h is the scale factor, Rm and Rp are the radii of an arbitrary particle in both scaled and physical models (m), ρm and ρp are the densities of scaled and physical models (kg·m−3), um and up are the displacements of scaled and physical models (m), Fm (um, Rm) and Fp (up, Rp) are the interaction forces from neighboring particles of scaled and physical models (N), and Qm (t) and Qp (t) are any external forces applied to the particle of scaled and physical models (N).It is well known that there are three basic quantities for mechanical systems in the SI system: length (L), mass (M), and time (T). After the three basic quantities are chosen, their corresponding conversion or scale factors between the physical and the scaled models can be arbitrarily set in principle.Density as a basic parameter of materials that can be directly obtained, and mass density (ρ) replaces mass as a basic quantity, but length and time remain as the other two basic quantities. The corresponding scale factors for these three quantities, λL, λT, and λρ are specifically chosen as [44](9)λL=h1,λT=h,λρ=1,i.e., the length scale factor is the same as the (spatial) scale factor, <i>h</i>, as expected, and the temporal (or time) scale is the same as the spatial factor, while the density is kept the same for both systems.Because the volume of the bin studied in this paper is large (16 m3) and the size of the material particles is relatively small, in order to avoid the “false crash” state of the CPU [25, 45], therefore, based on similarity principles mentioned above, in the simulation model, the sizes of the bin and the auger are both reduced to 0.25 of their sizes in a realistic model; similarly, the motion timescale of the material particles in the bin is also scaled to 0.25 of the original.On the basis of data in the literature [11, 13], a ratio of concentrate to roughage of 35 : 65 (based on the ratio of dry matter mass) is used to determine the compositions of the experimental rations. The total mass of the mixture produced by the EDEM particle factory is set as 1000 kg and the mass ratio of crop straw to globular corn grain to silage block as 10 : 35 : 55. All the particles in the mixture are assumed to be spherical, and from measurements of an actual mixture, the average particle radii for the three components are obtained as 2 mm, 4 mm, and 3 mm for crop straw, globular corn grain, and silage block, respectively. The Hertz–Mindlin contact model is used for the interaction between material particles and mixing device, and each crop particle model is defined as rigid particles in this paper. The three types of material particles in the discrete element model are shown in Figures 9(a)–9(c). To clearly show the motion of particles in the bin, the crop straw, globular corn grain, and silage block particles are colored yellow, red, and green, respectively, in EDEM, as shown in Figure 9(d).Figure 9
Discrete element model of material particles: (a) crop straw; (b) globular corn grain; (c) silage block; (d) partial view of material at a given time.
(a)
(b)
(c)
(d)The accuracy of DEM relies on the rationality of simulation parameters which should be calibrated before the simulation [41]. The simulation test was based on the results of physical test for the angle of repose, and the discrete element parameters were calibrated by measuring the angle of repose [46, 47]. The GEMM (generic EDEM material model) database contains thousands of material models representing a wide range of materials such as rocks, soils, and ores. It is particularly suited for engineers in the mining, construction, and agriculture industries who design such heavy equipment handling materials. The key contact property parameters including coefficient of restitution, coefficient of static friction, and coefficient of rolling friction can be obtained from the GEMM database [48, 49]. To increase the efficiency of particle generation, three different virtual factories of the static BOX type are created for the three types of material particles. After a start-up time of spiral auger rotation of 0.5 s, the total simulation time is set as 100 s. Therefore, in this paper, according to the GEMM database and data in the literature [50, 51], the material mechanical properties and the material contact property parameters are shown in Tables 1 and 2, respectively.Table 1
Material mechanical properties.
Material
Poisson’s ratio
Shear modulus (Pa)
Density (kg m−3)
Crop straw
0.30
6 × 108
57
Globular corn grain
0.40
1.37 × 108
1256
Silage block
0.30
2.1 × 107
420
Mixing device
0.29
7.7 × 1010
7850Table 2
Material contact property parameters used in the simulation.
Interaction material
Restitution
Static friction
Rolling friction
Crop straw-crop straw
0.33
0.38
0.01
Crop straw-globular corn grain
0.44
0.46
0.11
Crop straw-silage block
0.30
0.76
0.17
Crop straw-mixing device
0.34
0.42
0.01
Globular corn grain-globular corn grain
0.25
0.70
0.05
Globular corn grain-silage block
0.39
0.54
0.15
Globular corn grain-mixing device
0.30
0.50
0.07
Silage block-silage block
0.37
1.29
0.21
Silage block-mixing device
0.41
1.08
0.13


## 3.2.4. Factors and Evaluation Index

The uniformity of mixing of the material particles is the most important index for assessing the performance of the mixing device [1, 7, 15]. In this paper, the grid method is used to establish a model for evaluating this index. The simulation area is divided into 6 × 6 × 6 = 216 grids by the binning method in the EDEM software, as shown in Figure 10.Figure 10
Simulation model of the grid bin.If the grid size of the simulation model is too small, there will be an insufficient number of particles in a single grid to provide an accurate reflection of their dispersion, while if the grid size is too large, the results will not have statistical significance. Therefore, each grid cell should contain a sufficient number of material particles to ensure the validity of the grid cell data statistics, and it is necessary to perform a preliminary screening of the defined grid and remove those grid cells with less than 1000 material particles. The ratio of sample particles to the total number of particles in each grid cell is the sampling ratio, and the ratio of sample particles to the total number of particles in all materials is the overall ratio. The standard deviation of the sample is calculated as(10)S=∑i=1nXi−X¯2n−1,and the mixing uniformity of material particles is calculated as(11)M=1−SX¯×100%,where n is the number of grid cells, Xi is the sampling ratio of the sample particles, and X¯ is the overall ratio of the sample particles.According to the literature [11, 13, 14, 16], the factors affecting the performance of the mixing device include the material mixing time, the loading rate, and the auger speed. From the calculation of the critical speed of the auger described above and from practical experience with ration mixers, the following ranges of values of these factors are adopted in the simulations: a mixing time of 30–150 s, a loading rate of 50–90%, and an auger speed of 8–40 rev·min−1. The three-factor, five-level simulations are designed to analyze the mixing of material particles, with the uniformity of mixing of globular corn grains being taken as the index for evaluating mixing performance. The three factors and five levels are shown in Table 3.Table 3
Factors and levels of mixing simulation.
Level
Factor
Mixing time (s)
Loading rate (%)
Auger speed (rev·min−1)
1
30
50
8
2
60
60
16
3
90
70
24
4
120
80
32
5
150
90
40


## 4. Results and Discussions

### 4.1. Loading Characteristics

#### 4.1.1. Cutting Process of the Silage

With the material properties, contact type, and boundary constraints on the finite element model of the silage and loading cutter roller as described above, the keyword file is exported to the LS-DYNA Solver after the control parameters have been set, and the postprocessing analysis of the simulation result file D3Plot is carried out in LS-PrePost. Figure11 shows the effect of silage cutting with the loading cutter roller at four different times (20 ms, 60 ms, 100 ms, and 200 ms).Figure 11
Equivalent plastic strain nephograms of silage at different times: (a) 20 ms; (b) 60 ms; (c) 100 ms; (d) 200 ms. (e) Enlarged partial top view of silage cutting at 200 ms.
(a)
(b)
(c)
(d)
(e)To facilitate observation of the strain on the silage contact surface, the loading cutter roller is hidden as shown in Figures11(a) and 11(b). As can be seen from Figure 11(a), at a running time of 20 ms, the cutter roller starts to come in contact with the silage and the latter is subjected to the cutting force of the outer edges of the cutters. At this time, the maximum equivalent strain value is 0.15. From Figure 11(b), at 60 ms, it can be seen that as the cutter roller rotates downward, the contact area between the silage and the surface of the cutter roller gradually increases, and the shear force on the silage from the cutters increases accordingly. As can be seen from Figure 11(c), at a running time of 100 ms, the cutters have completely cut into the silage, the internal damaged area of the silage has been enlarged, and the deformations of the silage cutting layer have clearly increased compared with the previous time. When the deformation becomes sufficiently large, the silage is destroyed, with its material particles being torn away and displaced along the tangent direction of the cutters. These particles diffuse to both sides of the spiral blade of the cutter roller under the centrifugal force. The particles bulge and are subject to a maximum equivalent strain value of 0.34. As can be seen from Figures 11(d) and 11(e), at a running time of 200 ms, the roller has cut out the silage model area, and the spattered silage particles continue to undergo throwing motion due to inertia and the disturbance force. From an overview of this silage cutting process, it can be seen that, under the double action of internal extrusion and shearing by the spiral blade of the loading cutter roller, the deformation and displacement of the silage material are concentrated in the contact area between the silage and the roller.


#### 4.1.2. Variation of Cutting Force and Power Consumption

The variation of the cutting force during silage cutting by the loading cutter roller is shown in Figure12. As can be seen, the cutting force gradually increases from zero and then falls again to zero when the roller has cut out the silage area. With increasing contact area between the roller and the silage, the cutting force shows a sharp upward trend when the elastic deformation stage of the silage material has been overcome. When the bonding between the silage particles fails, the cutting force reaches a peak value of about 1.24 kN. Once the cutters have become completely immersed in the silage, the cutting force fluctuates up and down within a certain range. Since the loading cutter roller has partially cut out the silage material, the bulk density of the silage material also decreases with the continued cutting, and its deformation range begins to decrease. The silage material then starts to soften and the cutting force shows a decreasing trend. The force fluctuations during the whole cutting process from cutting into the silage to cutting out the silage are relatively small and the cutting process is relatively stable.Figure 12
Variation of cutting force with time.The total energy consumption of silage cutting includes the kinetic energy consumed to maintain the rotation of the loading cutter roller and the internal energy consumed by the interaction between the roller and the silage. As shown in Figure13, in the initial stage of silage cutting, as the roller and the silage come into contact with each other, silage deformation occurs and the power consumption increases rapidly. Once the roller has cut out the silage area, it idles; at the same time, the binding force between silage particles decreases because the silage has been broken, and the power consumption decreases. It can be seen from Figure 13 that the maximum power consumption over the whole silage cutting process is 3.87 kW, while the power consumption of the roller fluctuates around 2.69 kW during stable silage cutting. The ANSYS/LS-DYNA simulation could thus be used to estimate the power consumption of silage cutting, thereby providing a practical reference for optimization of the cutter roller structure.Figure 13
Variation of power consumption with time.


#### 4.1.3. Verification of Cutting Power Consumption

According to an empirical formula from an agricultural machinery design manual [52], the power consumption (kW) required for cutting the silage by the loading cutter roller is given by(12)Nr=Qr102000Lw+Hη,where Qr is the silage productivity (kg·min−1), L is the axial length of the roller (mm), w is the resistance coefficient of the movement of the silage along the shell of the loading device (=1.2 in this paper), H is the height through which the silage is lifted (=5000 mm), and η is a correction factor (=1.4). Qr is given by(13)Qr=π24D−2λ2−d2φtnγC×10−10,where D is the outer diameter of the spiral blade of the roller (=500 mm), λ is the clearance between the spiral blade and the shell (=8 mm), d is the inner diameter of the roller (=400 mm), φ is the material loading rate of the loading device (=0.4), t is the pitch of the spiral blade (=285 mm), n is the rotational speed of the loading cutter roller (=240 rev·min−1), γ is the bulk density of the silage material (=420 kg·m−3), and C is the inclined conveying coefficient of the spiral blade (=0.75).By substituting the values of the respective quantities into equations (12) and (13), we can get Nr = 2.78 kW. This value from the empirical formula is basically consistent with the value of 2.69 kW from the ANSYS/LS-DYNA simulation, and therefore, the finite element cutting model can be used to predict the power consumption of silage cutting by the loading cutter roller, which verifies the feasibility of simulation using the SPH-FEM coupling algorithm.


### 4.2. Mixing Characteristics

#### 4.2.1. Paths of Material Particles in the Axial Direction

To analyze the motion of the material particles, the path of a given particle is extracted, as shown in Figure14. The particle interacts with the group of particles surrounding it and tumbles along the tangential direction of the spiral auger blade as well as being lifted along the axial direction of the auger. When the particle reaches a certain height, the friction with the side wall of the bin is overcome and the particle falls under its own gravity. The three types of material particles collide with each other and no particle is able to remain stationary at a given location. Thus, the particles move up and down from the bottom of the bin and the top of the spiral auger to various positions. Under the rotating action of the spiral auger, the particles undergo cyclical movement. Therefore, the mixing process involves shear, diffusion, and convection of particles of different materials.Figure 14
Path of a given particle.


#### 4.2.2. Radial Material Flow at Different Auger Speeds

The flow of material in the bin of the mixing device is simulated for a material mixing time of 100 s and a loading rate of 70% at different auger speeds of 15 rev·min−1, 25 rev·min−1, and 35 rev·min−1. The clipping function in the postprocessing module of the EDEM software is used to generate two-dimensional cross-sections along the direction perpendicular to the auger shaft to provide a visual representation of the mixing of material particles in the bin, as shown in Figure 15.Figure 15
Radial material flow in the bin at different auger speeds: (a) 15 rev·min−1; (b) 25 rev·min−1; and (c) 35 rev·min−1.
(a)
(b)
(c)As can be seen from Figure15(a), when the auger speed is 15 rev·min−1, the overall intensity of mixing of the particles is weak, the three types of particles fill the bin space, and particles of different densities and sizes easily deposit and aggregate under the action of gravity. The particle motion is mainly caused by collisions with the wall of the bin and with other particles. When particles move axially, the effects of extrusion and of friction with surrounding particles and with the wall of the bin are greater, reducing the intensity of particle diffusion and convective motion; relative motion between particles during the circumferential mixing process is also more difficult, and as a consequence, interparticle mixing is insufficient. Thus, a low auger speed is unfavorable for uniform mixing of particles.For auger speeds of 25 rev·min−1 and 35 rev·min−1, the form in which the mixture material displayed is changed in the postprocessing interface of EDEM from the default particle state to vector stream (two-dimensional cross-sections at the position of the discharging device near the bottom wall of the bin are obtained). Figure 15(b) shows that at an auger speed of 25 rev·min−1, the density distributions of the three types of particles are substantially the same and these particles are uniformly distributed in the space of the bin. As the mixing process continues, convective mixing of the particles is enhanced. Pushed by the spiral auger blade, the particles are continuously lifted and tumbled, diffusing and mixing with surrounding particle groups. The mixing between the particles thereby tends to a dynamic equilibrium.As shown in Figure15(c), at an auger speed of 35 rev·min−1, the overall intensity of mixing of the particles in the bin is greater. The linear velocities of the particles around the spiral auger blade are significantly higher than in other regions, resulting in a “boiling” effect. The number of particles around the spiral auger shaft decreases, reducing the retention time of the particles in the bin. The high rotational speed of the auger tends to cause segregation and gradation of material particles, reducing the possibility of contact among the three types of particles. Material particles readily accumulate at the edge of the bin wall and the distribution ceases to be uniform. Thus, it is important to take care when increasing the auger speed of the mixing device, and since if it is too high, the uniformity of distribution of the material particles will be reduced.


#### 4.2.3. Influence of Each Factor on Mixing Uniformity

In the single-factor simulations, the other two factors are held constant at their central-level values (i.e., those for level 3), and at the end of each simulation, the mixing uniformities of the 15 sets of simulation examples are noted. The influence of each factor on mixing uniformity is shown in Figure16.Figure 16
Influence of each simulation factor on mixing uniformity.With the material loading rate and the auger speed held constant, the influence of the material mixing time on mixing uniformity is shown by the blue line in Figure16. The mixing uniformity first increases and then decreases with increasing mixing time. In the initial stage of mixing, the particles are lifted vertically along the shaft of the spiral auger while rolling in the circumferential direction. They are then gradually mixed by shear, convection, and diffusion, leading to a uniform distribution of the three types of material particles. As the mixing time continues to increase, however, owing to differences in physical properties between the three types of particles, they become prone to segregation, resulting in a reduction in mixing uniformity to 70.7%.With the mixing time and the auger speed held constant, the influence of the material loading rate on mixing uniformity is shown by the green line in Figure16. At a low material loading rate (50%), the mixing uniformity is 77.8%. As the loading rate increases, the shear and convective motions generated by the particle extrusion layer become stronger and the probability of particles coming into contact with each other increases, which enhances the cross-mixing of particles and increases the mixing uniformity to 90.7%. When the material loading rate is increased further, the remaining space in the bin becomes smaller, which hinders the relative motion of particles in the circumferential mixing process. Owing to the increase of the material mass in the bin, the extrusion resistance of the particles during axial mixing also increases, the intensities of particle convection and diffusion decrease, and the mixing uniformity drops to 74.9%.With the material mixing time and loading rate held constant, the influence of the auger speed on mixing uniformity is shown by the red line in Figure16. As the auger speed increases from a low value, particles in different regions of the bin are subjected to a combination of shear, convection, and diffusion, and a circumferential and axial three-dimensional mixing motion is generated, leading to a highly uniform particle distribution, with a mixing uniformity of 91.2% when the auger speed is at its central level (24 rev·min−1). A further increase in the auger speed enhances the pushing effect of the spiral auger blade on the particles and increases the centrifugal force on the particles, causing particles with the same physical properties to congregate. This leads to separation and gradation of the three different types of particles, and the mixing uniformity falls sharply to 67.1%.The simulation results show that the effects of the three experimental factors on the mixing uniformity are similar in that the uniformity first increases and then decreases with increasing value of each factor. It can be seen from Figure16, however, that the effect of the auger speed on the mixing uniformity is more significant than that of the mixing time, which in turn has a more significant effect than the material loading rate.


#### 4.2.4. Performance Test of the Prototype

To verify the reliability of the theoretical analysis and numerical simulations, in December 2018, performance tests of silage loading and material mixing with a prototype self-propelled TMR mixer (Figure17) were conducted at the Zhongli Dairy Breeding Cooperative of Feicheng City, Shandong Province. The test material was uniformly formulated in accordance with Chinese standard JB/T 11438-2013 “Total Mixed Ration Mixer” [53], with the mass fractions of a concentrate (including maize flour 20%, soybean meal 10%, and corn grain 5% as tracers) and roughage (including silage 55% and straw 10%) being 35% and 65%, respectively. The testing instruments were principally an electronic scale (model TCS-01, accuracy 0.5 g), a stopwatch (accuracy 1 s), and sampling bags. According to the test method specified in the standard, the mixing uniformity in the bin of the mixing device was measured for a material mixing time of 8 min and a loading rate of 75% at different auger speeds of 22 rev·min−1, 24 rev·min−1, and 26 rev·min−1. After the performance test of the practical prototype, the mixing uniformity of material was calculated using the formula (11) [11–13]. The mixing uniformity for the three times were 90.9%, 92.3%, and 90.1%, respectively, with an average mixing uniformity of 91.1%, which is a suitable value for practical application in the animal husbandry industry.Figure 17
Field test of the prototype of the self-propelled TMR mixer.Based on simulation and performance test results, with the practical production, an auger speed of 22–26 rev·min−1, a material mixing time of 6–8 min (the simulation mixing time of the material in Table 3 should be scaled to the actual mixing time of the material), and a loading rate of 70–75% give a mixing uniformity greater than 90%, meeting the requirements for material mixing stipulated by the national standard [53]. In addition, the optimal auger speed of 22–26 rev·min−1 obtained from the simulation is within the theoretical critical speed range of 15–32 rev·min−1, which indirectly confirms the reliability of the DEM in the analysis of material-mixing characteristics. These results indicate that the proposed method could provide reference values for working parameters when optimizing the mixing performance of self-propelled TMR mixers.


## 4.1. Loading Characteristics

### 4.1.1. Cutting Process of the Silage

With the material properties, contact type, and boundary constraints on the finite element model of the silage and loading cutter roller as described above, the keyword file is exported to the LS-DYNA Solver after the control parameters have been set, and the postprocessing analysis of the simulation result file D3Plot is carried out in LS-PrePost. Figure11 shows the effect of silage cutting with the loading cutter roller at four different times (20 ms, 60 ms, 100 ms, and 200 ms).Figure 11
Equivalent plastic strain nephograms of silage at different times: (a) 20 ms; (b) 60 ms; (c) 100 ms; (d) 200 ms. (e) Enlarged partial top view of silage cutting at 200 ms.
(a)
(b)
(c)
(d)
(e)To facilitate observation of the strain on the silage contact surface, the loading cutter roller is hidden as shown in Figures11(a) and 11(b). As can be seen from Figure 11(a), at a running time of 20 ms, the cutter roller starts to come in contact with the silage and the latter is subjected to the cutting force of the outer edges of the cutters. At this time, the maximum equivalent strain value is 0.15. From Figure 11(b), at 60 ms, it can be seen that as the cutter roller rotates downward, the contact area between the silage and the surface of the cutter roller gradually increases, and the shear force on the silage from the cutters increases accordingly. As can be seen from Figure 11(c), at a running time of 100 ms, the cutters have completely cut into the silage, the internal damaged area of the silage has been enlarged, and the deformations of the silage cutting layer have clearly increased compared with the previous time. When the deformation becomes sufficiently large, the silage is destroyed, with its material particles being torn away and displaced along the tangent direction of the cutters. These particles diffuse to both sides of the spiral blade of the cutter roller under the centrifugal force. The particles bulge and are subject to a maximum equivalent strain value of 0.34. As can be seen from Figures 11(d) and 11(e), at a running time of 200 ms, the roller has cut out the silage model area, and the spattered silage particles continue to undergo throwing motion due to inertia and the disturbance force. From an overview of this silage cutting process, it can be seen that, under the double action of internal extrusion and shearing by the spiral blade of the loading cutter roller, the deformation and displacement of the silage material are concentrated in the contact area between the silage and the roller.


### 4.1.2. Variation of Cutting Force and Power Consumption

The variation of the cutting force during silage cutting by the loading cutter roller is shown in Figure12. As can be seen, the cutting force gradually increases from zero and then falls again to zero when the roller has cut out the silage area. With increasing contact area between the roller and the silage, the cutting force shows a sharp upward trend when the elastic deformation stage of the silage material has been overcome. When the bonding between the silage particles fails, the cutting force reaches a peak value of about 1.24 kN. Once the cutters have become completely immersed in the silage, the cutting force fluctuates up and down within a certain range. Since the loading cutter roller has partially cut out the silage material, the bulk density of the silage material also decreases with the continued cutting, and its deformation range begins to decrease. The silage material then starts to soften and the cutting force shows a decreasing trend. The force fluctuations during the whole cutting process from cutting into the silage to cutting out the silage are relatively small and the cutting process is relatively stable.Figure 12
Variation of cutting force with time.The total energy consumption of silage cutting includes the kinetic energy consumed to maintain the rotation of the loading cutter roller and the internal energy consumed by the interaction between the roller and the silage. As shown in Figure13, in the initial stage of silage cutting, as the roller and the silage come into contact with each other, silage deformation occurs and the power consumption increases rapidly. Once the roller has cut out the silage area, it idles; at the same time, the binding force between silage particles decreases because the silage has been broken, and the power consumption decreases. It can be seen from Figure 13 that the maximum power consumption over the whole silage cutting process is 3.87 kW, while the power consumption of the roller fluctuates around 2.69 kW during stable silage cutting. The ANSYS/LS-DYNA simulation could thus be used to estimate the power consumption of silage cutting, thereby providing a practical reference for optimization of the cutter roller structure.Figure 13
Variation of power consumption with time.


### 4.1.3. Verification of Cutting Power Consumption

According to an empirical formula from an agricultural machinery design manual [52], the power consumption (kW) required for cutting the silage by the loading cutter roller is given by(12)Nr=Qr102000Lw+Hη,where Qr is the silage productivity (kg·min−1), L is the axial length of the roller (mm), w is the resistance coefficient of the movement of the silage along the shell of the loading device (=1.2 in this paper), H is the height through which the silage is lifted (=5000 mm), and η is a correction factor (=1.4). Qr is given by(13)Qr=π24D−2λ2−d2φtnγC×10−10,where D is the outer diameter of the spiral blade of the roller (=500 mm), λ is the clearance between the spiral blade and the shell (=8 mm), d is the inner diameter of the roller (=400 mm), φ is the material loading rate of the loading device (=0.4), t is the pitch of the spiral blade (=285 mm), n is the rotational speed of the loading cutter roller (=240 rev·min−1), γ is the bulk density of the silage material (=420 kg·m−3), and C is the inclined conveying coefficient of the spiral blade (=0.75).By substituting the values of the respective quantities into equations (12) and (13), we can get Nr = 2.78 kW. This value from the empirical formula is basically consistent with the value of 2.69 kW from the ANSYS/LS-DYNA simulation, and therefore, the finite element cutting model can be used to predict the power consumption of silage cutting by the loading cutter roller, which verifies the feasibility of simulation using the SPH-FEM coupling algorithm.


## 4.1.1. Cutting Process of the Silage

With the material properties, contact type, and boundary constraints on the finite element model of the silage and loading cutter roller as described above, the keyword file is exported to the LS-DYNA Solver after the control parameters have been set, and the postprocessing analysis of the simulation result file D3Plot is carried out in LS-PrePost. Figure11 shows the effect of silage cutting with the loading cutter roller at four different times (20 ms, 60 ms, 100 ms, and 200 ms).Figure 11
Equivalent plastic strain nephograms of silage at different times: (a) 20 ms; (b) 60 ms; (c) 100 ms; (d) 200 ms. (e) Enlarged partial top view of silage cutting at 200 ms.
(a)
(b)
(c)
(d)
(e)To facilitate observation of the strain on the silage contact surface, the loading cutter roller is hidden as shown in Figures11(a) and 11(b). As can be seen from Figure 11(a), at a running time of 20 ms, the cutter roller starts to come in contact with the silage and the latter is subjected to the cutting force of the outer edges of the cutters. At this time, the maximum equivalent strain value is 0.15. From Figure 11(b), at 60 ms, it can be seen that as the cutter roller rotates downward, the contact area between the silage and the surface of the cutter roller gradually increases, and the shear force on the silage from the cutters increases accordingly. As can be seen from Figure 11(c), at a running time of 100 ms, the cutters have completely cut into the silage, the internal damaged area of the silage has been enlarged, and the deformations of the silage cutting layer have clearly increased compared with the previous time. When the deformation becomes sufficiently large, the silage is destroyed, with its material particles being torn away and displaced along the tangent direction of the cutters. These particles diffuse to both sides of the spiral blade of the cutter roller under the centrifugal force. The particles bulge and are subject to a maximum equivalent strain value of 0.34. As can be seen from Figures 11(d) and 11(e), at a running time of 200 ms, the roller has cut out the silage model area, and the spattered silage particles continue to undergo throwing motion due to inertia and the disturbance force. From an overview of this silage cutting process, it can be seen that, under the double action of internal extrusion and shearing by the spiral blade of the loading cutter roller, the deformation and displacement of the silage material are concentrated in the contact area between the silage and the roller.


## 4.1.2. Variation of Cutting Force and Power Consumption

The variation of the cutting force during silage cutting by the loading cutter roller is shown in Figure12. As can be seen, the cutting force gradually increases from zero and then falls again to zero when the roller has cut out the silage area. With increasing contact area between the roller and the silage, the cutting force shows a sharp upward trend when the elastic deformation stage of the silage material has been overcome. When the bonding between the silage particles fails, the cutting force reaches a peak value of about 1.24 kN. Once the cutters have become completely immersed in the silage, the cutting force fluctuates up and down within a certain range. Since the loading cutter roller has partially cut out the silage material, the bulk density of the silage material also decreases with the continued cutting, and its deformation range begins to decrease. The silage material then starts to soften and the cutting force shows a decreasing trend. The force fluctuations during the whole cutting process from cutting into the silage to cutting out the silage are relatively small and the cutting process is relatively stable.Figure 12
Variation of cutting force with time.The total energy consumption of silage cutting includes the kinetic energy consumed to maintain the rotation of the loading cutter roller and the internal energy consumed by the interaction between the roller and the silage. As shown in Figure13, in the initial stage of silage cutting, as the roller and the silage come into contact with each other, silage deformation occurs and the power consumption increases rapidly. Once the roller has cut out the silage area, it idles; at the same time, the binding force between silage particles decreases because the silage has been broken, and the power consumption decreases. It can be seen from Figure 13 that the maximum power consumption over the whole silage cutting process is 3.87 kW, while the power consumption of the roller fluctuates around 2.69 kW during stable silage cutting. The ANSYS/LS-DYNA simulation could thus be used to estimate the power consumption of silage cutting, thereby providing a practical reference for optimization of the cutter roller structure.Figure 13
Variation of power consumption with time.


## 4.1.3. Verification of Cutting Power Consumption

According to an empirical formula from an agricultural machinery design manual [52], the power consumption (kW) required for cutting the silage by the loading cutter roller is given by(12)Nr=Qr102000Lw+Hη,where Qr is the silage productivity (kg·min−1), L is the axial length of the roller (mm), w is the resistance coefficient of the movement of the silage along the shell of the loading device (=1.2 in this paper), H is the height through which the silage is lifted (=5000 mm), and η is a correction factor (=1.4). Qr is given by(13)Qr=π24D−2λ2−d2φtnγC×10−10,where D is the outer diameter of the spiral blade of the roller (=500 mm), λ is the clearance between the spiral blade and the shell (=8 mm), d is the inner diameter of the roller (=400 mm), φ is the material loading rate of the loading device (=0.4), t is the pitch of the spiral blade (=285 mm), n is the rotational speed of the loading cutter roller (=240 rev·min−1), γ is the bulk density of the silage material (=420 kg·m−3), and C is the inclined conveying coefficient of the spiral blade (=0.75).By substituting the values of the respective quantities into equations (12) and (13), we can get Nr = 2.78 kW. This value from the empirical formula is basically consistent with the value of 2.69 kW from the ANSYS/LS-DYNA simulation, and therefore, the finite element cutting model can be used to predict the power consumption of silage cutting by the loading cutter roller, which verifies the feasibility of simulation using the SPH-FEM coupling algorithm.


## 4.2. Mixing Characteristics

### 4.2.1. Paths of Material Particles in the Axial Direction

To analyze the motion of the material particles, the path of a given particle is extracted, as shown in Figure14. The particle interacts with the group of particles surrounding it and tumbles along the tangential direction of the spiral auger blade as well as being lifted along the axial direction of the auger. When the particle reaches a certain height, the friction with the side wall of the bin is overcome and the particle falls under its own gravity. The three types of material particles collide with each other and no particle is able to remain stationary at a given location. Thus, the particles move up and down from the bottom of the bin and the top of the spiral auger to various positions. Under the rotating action of the spiral auger, the particles undergo cyclical movement. Therefore, the mixing process involves shear, diffusion, and convection of particles of different materials.Figure 14
Path of a given particle.


### 4.2.2. Radial Material Flow at Different Auger Speeds

The flow of material in the bin of the mixing device is simulated for a material mixing time of 100 s and a loading rate of 70% at different auger speeds of 15 rev·min−1, 25 rev·min−1, and 35 rev·min−1. The clipping function in the postprocessing module of the EDEM software is used to generate two-dimensional cross-sections along the direction perpendicular to the auger shaft to provide a visual representation of the mixing of material particles in the bin, as shown in Figure 15.Figure 15
Radial material flow in the bin at different auger speeds: (a) 15 rev·min−1; (b) 25 rev·min−1; and (c) 35 rev·min−1.
(a)
(b)
(c)As can be seen from Figure15(a), when the auger speed is 15 rev·min−1, the overall intensity of mixing of the particles is weak, the three types of particles fill the bin space, and particles of different densities and sizes easily deposit and aggregate under the action of gravity. The particle motion is mainly caused by collisions with the wall of the bin and with other particles. When particles move axially, the effects of extrusion and of friction with surrounding particles and with the wall of the bin are greater, reducing the intensity of particle diffusion and convective motion; relative motion between particles during the circumferential mixing process is also more difficult, and as a consequence, interparticle mixing is insufficient. Thus, a low auger speed is unfavorable for uniform mixing of particles.For auger speeds of 25 rev·min−1 and 35 rev·min−1, the form in which the mixture material displayed is changed in the postprocessing interface of EDEM from the default particle state to vector stream (two-dimensional cross-sections at the position of the discharging device near the bottom wall of the bin are obtained). Figure 15(b) shows that at an auger speed of 25 rev·min−1, the density distributions of the three types of particles are substantially the same and these particles are uniformly distributed in the space of the bin. As the mixing process continues, convective mixing of the particles is enhanced. Pushed by the spiral auger blade, the particles are continuously lifted and tumbled, diffusing and mixing with surrounding particle groups. The mixing between the particles thereby tends to a dynamic equilibrium.As shown in Figure15(c), at an auger speed of 35 rev·min−1, the overall intensity of mixing of the particles in the bin is greater. The linear velocities of the particles around the spiral auger blade are significantly higher than in other regions, resulting in a “boiling” effect. The number of particles around the spiral auger shaft decreases, reducing the retention time of the particles in the bin. The high rotational speed of the auger tends to cause segregation and gradation of material particles, reducing the possibility of contact among the three types of particles. Material particles readily accumulate at the edge of the bin wall and the distribution ceases to be uniform. Thus, it is important to take care when increasing the auger speed of the mixing device, and since if it is too high, the uniformity of distribution of the material particles will be reduced.


### 4.2.3. Influence of Each Factor on Mixing Uniformity

In the single-factor simulations, the other two factors are held constant at their central-level values (i.e., those for level 3), and at the end of each simulation, the mixing uniformities of the 15 sets of simulation examples are noted. The influence of each factor on mixing uniformity is shown in Figure16.Figure 16
Influence of each simulation factor on mixing uniformity.With the material loading rate and the auger speed held constant, the influence of the material mixing time on mixing uniformity is shown by the blue line in Figure16. The mixing uniformity first increases and then decreases with increasing mixing time. In the initial stage of mixing, the particles are lifted vertically along the shaft of the spiral auger while rolling in the circumferential direction. They are then gradually mixed by shear, convection, and diffusion, leading to a uniform distribution of the three types of material particles. As the mixing time continues to increase, however, owing to differences in physical properties between the three types of particles, they become prone to segregation, resulting in a reduction in mixing uniformity to 70.7%.With the mixing time and the auger speed held constant, the influence of the material loading rate on mixing uniformity is shown by the green line in Figure16. At a low material loading rate (50%), the mixing uniformity is 77.8%. As the loading rate increases, the shear and convective motions generated by the particle extrusion layer become stronger and the probability of particles coming into contact with each other increases, which enhances the cross-mixing of particles and increases the mixing uniformity to 90.7%. When the material loading rate is increased further, the remaining space in the bin becomes smaller, which hinders the relative motion of particles in the circumferential mixing process. Owing to the increase of the material mass in the bin, the extrusion resistance of the particles during axial mixing also increases, the intensities of particle convection and diffusion decrease, and the mixing uniformity drops to 74.9%.With the material mixing time and loading rate held constant, the influence of the auger speed on mixing uniformity is shown by the red line in Figure16. As the auger speed increases from a low value, particles in different regions of the bin are subjected to a combination of shear, convection, and diffusion, and a circumferential and axial three-dimensional mixing motion is generated, leading to a highly uniform particle distribution, with a mixing uniformity of 91.2% when the auger speed is at its central level (24 rev·min−1). A further increase in the auger speed enhances the pushing effect of the spiral auger blade on the particles and increases the centrifugal force on the particles, causing particles with the same physical properties to congregate. This leads to separation and gradation of the three different types of particles, and the mixing uniformity falls sharply to 67.1%.The simulation results show that the effects of the three experimental factors on the mixing uniformity are similar in that the uniformity first increases and then decreases with increasing value of each factor. It can be seen from Figure16, however, that the effect of the auger speed on the mixing uniformity is more significant than that of the mixing time, which in turn has a more significant effect than the material loading rate.


### 4.2.4. Performance Test of the Prototype

To verify the reliability of the theoretical analysis and numerical simulations, in December 2018, performance tests of silage loading and material mixing with a prototype self-propelled TMR mixer (Figure17) were conducted at the Zhongli Dairy Breeding Cooperative of Feicheng City, Shandong Province. The test material was uniformly formulated in accordance with Chinese standard JB/T 11438-2013 “Total Mixed Ration Mixer” [53], with the mass fractions of a concentrate (including maize flour 20%, soybean meal 10%, and corn grain 5% as tracers) and roughage (including silage 55% and straw 10%) being 35% and 65%, respectively. The testing instruments were principally an electronic scale (model TCS-01, accuracy 0.5 g), a stopwatch (accuracy 1 s), and sampling bags. According to the test method specified in the standard, the mixing uniformity in the bin of the mixing device was measured for a material mixing time of 8 min and a loading rate of 75% at different auger speeds of 22 rev·min−1, 24 rev·min−1, and 26 rev·min−1. After the performance test of the practical prototype, the mixing uniformity of material was calculated using the formula (11) [11–13]. The mixing uniformity for the three times were 90.9%, 92.3%, and 90.1%, respectively, with an average mixing uniformity of 91.1%, which is a suitable value for practical application in the animal husbandry industry.Figure 17
Field test of the prototype of the self-propelled TMR mixer.Based on simulation and performance test results, with the practical production, an auger speed of 22–26 rev·min−1, a material mixing time of 6–8 min (the simulation mixing time of the material in Table 3 should be scaled to the actual mixing time of the material), and a loading rate of 70–75% give a mixing uniformity greater than 90%, meeting the requirements for material mixing stipulated by the national standard [53]. In addition, the optimal auger speed of 22–26 rev·min−1 obtained from the simulation is within the theoretical critical speed range of 15–32 rev·min−1, which indirectly confirms the reliability of the DEM in the analysis of material-mixing characteristics. These results indicate that the proposed method could provide reference values for working parameters when optimizing the mixing performance of self-propelled TMR mixers.


## 4.2.1. Paths of Material Particles in the Axial Direction

To analyze the motion of the material particles, the path of a given particle is extracted, as shown in Figure14. The particle interacts with the group of particles surrounding it and tumbles along the tangential direction of the spiral auger blade as well as being lifted along the axial direction of the auger. When the particle reaches a certain height, the friction with the side wall of the bin is overcome and the particle falls under its own gravity. The three types of material particles collide with each other and no particle is able to remain stationary at a given location. Thus, the particles move up and down from the bottom of the bin and the top of the spiral auger to various positions. Under the rotating action of the spiral auger, the particles undergo cyclical movement. Therefore, the mixing process involves shear, diffusion, and convection of particles of different materials.Figure 14
Path of a given particle.


## 4.2.2. Radial Material Flow at Different Auger Speeds

The flow of material in the bin of the mixing device is simulated for a material mixing time of 100 s and a loading rate of 70% at different auger speeds of 15 rev·min−1, 25 rev·min−1, and 35 rev·min−1. The clipping function in the postprocessing module of the EDEM software is used to generate two-dimensional cross-sections along the direction perpendicular to the auger shaft to provide a visual representation of the mixing of material particles in the bin, as shown in Figure 15.Figure 15
Radial material flow in the bin at different auger speeds: (a) 15 rev·min−1; (b) 25 rev·min−1; and (c) 35 rev·min−1.
(a)
(b)
(c)As can be seen from Figure15(a), when the auger speed is 15 rev·min−1, the overall intensity of mixing of the particles is weak, the three types of particles fill the bin space, and particles of different densities and sizes easily deposit and aggregate under the action of gravity. The particle motion is mainly caused by collisions with the wall of the bin and with other particles. When particles move axially, the effects of extrusion and of friction with surrounding particles and with the wall of the bin are greater, reducing the intensity of particle diffusion and convective motion; relative motion between particles during the circumferential mixing process is also more difficult, and as a consequence, interparticle mixing is insufficient. Thus, a low auger speed is unfavorable for uniform mixing of particles.For auger speeds of 25 rev·min−1 and 35 rev·min−1, the form in which the mixture material displayed is changed in the postprocessing interface of EDEM from the default particle state to vector stream (two-dimensional cross-sections at the position of the discharging device near the bottom wall of the bin are obtained). Figure 15(b) shows that at an auger speed of 25 rev·min−1, the density distributions of the three types of particles are substantially the same and these particles are uniformly distributed in the space of the bin. As the mixing process continues, convective mixing of the particles is enhanced. Pushed by the spiral auger blade, the particles are continuously lifted and tumbled, diffusing and mixing with surrounding particle groups. The mixing between the particles thereby tends to a dynamic equilibrium.As shown in Figure15(c), at an auger speed of 35 rev·min−1, the overall intensity of mixing of the particles in the bin is greater. The linear velocities of the particles around the spiral auger blade are significantly higher than in other regions, resulting in a “boiling” effect. The number of particles around the spiral auger shaft decreases, reducing the retention time of the particles in the bin. The high rotational speed of the auger tends to cause segregation and gradation of material particles, reducing the possibility of contact among the three types of particles. Material particles readily accumulate at the edge of the bin wall and the distribution ceases to be uniform. Thus, it is important to take care when increasing the auger speed of the mixing device, and since if it is too high, the uniformity of distribution of the material particles will be reduced.


## 4.2.3. Influence of Each Factor on Mixing Uniformity

In the single-factor simulations, the other two factors are held constant at their central-level values (i.e., those for level 3), and at the end of each simulation, the mixing uniformities of the 15 sets of simulation examples are noted. The influence of each factor on mixing uniformity is shown in Figure16.Figure 16
Influence of each simulation factor on mixing uniformity.With the material loading rate and the auger speed held constant, the influence of the material mixing time on mixing uniformity is shown by the blue line in Figure16. The mixing uniformity first increases and then decreases with increasing mixing time. In the initial stage of mixing, the particles are lifted vertically along the shaft of the spiral auger while rolling in the circumferential direction. They are then gradually mixed by shear, convection, and diffusion, leading to a uniform distribution of the three types of material particles. As the mixing time continues to increase, however, owing to differences in physical properties between the three types of particles, they become prone to segregation, resulting in a reduction in mixing uniformity to 70.7%.With the mixing time and the auger speed held constant, the influence of the material loading rate on mixing uniformity is shown by the green line in Figure16. At a low material loading rate (50%), the mixing uniformity is 77.8%. As the loading rate increases, the shear and convective motions generated by the particle extrusion layer become stronger and the probability of particles coming into contact with each other increases, which enhances the cross-mixing of particles and increases the mixing uniformity to 90.7%. When the material loading rate is increased further, the remaining space in the bin becomes smaller, which hinders the relative motion of particles in the circumferential mixing process. Owing to the increase of the material mass in the bin, the extrusion resistance of the particles during axial mixing also increases, the intensities of particle convection and diffusion decrease, and the mixing uniformity drops to 74.9%.With the material mixing time and loading rate held constant, the influence of the auger speed on mixing uniformity is shown by the red line in Figure16. As the auger speed increases from a low value, particles in different regions of the bin are subjected to a combination of shear, convection, and diffusion, and a circumferential and axial three-dimensional mixing motion is generated, leading to a highly uniform particle distribution, with a mixing uniformity of 91.2% when the auger speed is at its central level (24 rev·min−1). A further increase in the auger speed enhances the pushing effect of the spiral auger blade on the particles and increases the centrifugal force on the particles, causing particles with the same physical properties to congregate. This leads to separation and gradation of the three different types of particles, and the mixing uniformity falls sharply to 67.1%.The simulation results show that the effects of the three experimental factors on the mixing uniformity are similar in that the uniformity first increases and then decreases with increasing value of each factor. It can be seen from Figure16, however, that the effect of the auger speed on the mixing uniformity is more significant than that of the mixing time, which in turn has a more significant effect than the material loading rate.


## 4.2.4. Performance Test of the Prototype

To verify the reliability of the theoretical analysis and numerical simulations, in December 2018, performance tests of silage loading and material mixing with a prototype self-propelled TMR mixer (Figure17) were conducted at the Zhongli Dairy Breeding Cooperative of Feicheng City, Shandong Province. The test material was uniformly formulated in accordance with Chinese standard JB/T 11438-2013 “Total Mixed Ration Mixer” [53], with the mass fractions of a concentrate (including maize flour 20%, soybean meal 10%, and corn grain 5% as tracers) and roughage (including silage 55% and straw 10%) being 35% and 65%, respectively. The testing instruments were principally an electronic scale (model TCS-01, accuracy 0.5 g), a stopwatch (accuracy 1 s), and sampling bags. According to the test method specified in the standard, the mixing uniformity in the bin of the mixing device was measured for a material mixing time of 8 min and a loading rate of 75% at different auger speeds of 22 rev·min−1, 24 rev·min−1, and 26 rev·min−1. After the performance test of the practical prototype, the mixing uniformity of material was calculated using the formula (11) [11–13]. The mixing uniformity for the three times were 90.9%, 92.3%, and 90.1%, respectively, with an average mixing uniformity of 91.1%, which is a suitable value for practical application in the animal husbandry industry.Figure 17
Field test of the prototype of the self-propelled TMR mixer.Based on simulation and performance test results, with the practical production, an auger speed of 22–26 rev·min−1, a material mixing time of 6–8 min (the simulation mixing time of the material in Table 3 should be scaled to the actual mixing time of the material), and a loading rate of 70–75% give a mixing uniformity greater than 90%, meeting the requirements for material mixing stipulated by the national standard [53]. In addition, the optimal auger speed of 22–26 rev·min−1 obtained from the simulation is within the theoretical critical speed range of 15–32 rev·min−1, which indirectly confirms the reliability of the DEM in the analysis of material-mixing characteristics. These results indicate that the proposed method could provide reference values for working parameters when optimizing the mixing performance of self-propelled TMR mixers.


## 5. Conclusions

(1)
For the loading device of a self-propelled TMR mixer, an SPH-FEM coupling algorithm has been used to establish a finite element model of the silage and the loading cutter roller. Numerical simulation of silage cutting has been carried out using ANSYS/LS-DYNA. The distribution of silage strain and the variations of cutting force and power consumption have been investigated, and the results show that silage is broken up mainly by extrusion and shear forces due to the spiral blade of the loading cutter roller. The maximum cutting force is 1.24 kN, and the power consumption fluctuates around 2.69 kW, with a maximum value of 3.87 kW. The simulation results have been verified by those obtained from an established empirical formula.
(2)
For the mixing device of the self-propelled TMR mixer, a Hertz–Mindlin dynamic contact model of three kinds of discrete element material particles has been established, the particle motion path has been simulated, and the material-mixing mechanism has been analyzed. The material moves along the axial direction of the spiral auger while tumbling along the circumferential direction, leading to a combination of shear, diffusive, and convective motion in different regions of the bin. Simulations of mixing performance have shown that the auger speed, mixing time, and loading rate have a descending order of influence on the mixing performance and that an auger speed of 22–26 rev·min−1, a mixing time of 6–8 min, and a loading rate of 70–75% give an optimal mixing effect.



---

*Source: 6875816-2020-03-13.xml*