3D Finite Element Simulation and Experimental Validation of a Mole Rat’s Digit Inspired Biomimetic Potato Digging Shovel

Featured Application: This work proves the feasibility and effectiveness of biomimetic macroscopic surface modiﬁcation method to reduce the draught force of potato digging shovels based on mole rat’s digits, which is also a reference and experience for the optimization and modiﬁcation of other soil engaging components. Abstract: To reduce the draught force of a traditional planar potato digging shovel (DZ), a biomimetic potato digging shovel (YS), inspired by the mole rat’s digits, is designed, using the biomimetic macroscopic surface modiﬁcation method. The ﬁnite element simulations, soil bin experiments, and ﬁeld experiments for DZ and YS are conducted to explore the factors affecting draught force and to verify the feasibility and effectiveness of the biomimetic potato digging shovel. Results show that the soil–shovel interaction models predict the draught force well, but the simulations for the soil rupture distance ratio need to be further improved. The studied factors all have a great inﬂuence on the draught force of DZ and YS and they follow the order of cutting speed > digging depth > mounting angle. For the single shovels, YS, compared with DZ, increases the draught force at a low mounting angle but decreases the draught force by over 8.41% when the mounting angle is higher than 30 ◦ ; for the grouped shovels, the draught force and fuel consumption of YS, compared with those of DZ, decline by over 13.33% and 9.18%, respectively. The reasons for the reduction in the draught force of YS are to make the soil mass tend to move upward and to change the soil’s state of motion and stress continually; thus, the compaction to the soil is reduced, and the soil becomes easier to be broken.


Introduction
Gansu province is one of the major potato production areas in China. The hills and mountains scattered in Gansu, however, enable only the medium and small potato harvesters to be applicable there. Potato digging shovels are the key power-consumption components in potato harvesting, whose consumption accounts for more than 27% [1]. Therefore, to improve the performance of planar potato digging shovels and, thus, reduce the draught force are of great economic and regional importance in Gansu province.
The existing methods to reduce the drought force include the fluid injection [2,3], vibration [4], electro-osmosis [3], magnetization [5], macroscopic or microcosmic surface modification [3,6], use of different materials or coatings [7] and so on. However, these methods require higher strength and stability from the components or machinery, which is contrary to the demand for simplified and miniaturized agricultural machinery in Gansu province.
Chen et al. [8] and Ren et al. [9] examined the characteristics of soil-burrowing animals' soil engaging organs in the 1990s and they found that the soil engaging surfaces of these animals were characterized by a spatial surface with convex domes, concave hollows, ridges, furrows or other non-smooth construction units, which was beneficial to the decrease in the compaction to soil as well as the adhesion between the soil and interface and finally the reduction in the drought force. Based on the above findings, scholars applied biomimetics method to improve the performance of soil-engaging components and obtained favorable results. For instance, the biomimetic non-smooth moldboard ploughs imitating the surface morphology of dung beetle reduced the drought force by 15-18% and was used in Northeastern China [10]. By equipping some electrodes with specific shape and arrangement on the surface of soil-engaging components, the drought force decreased by 15-40%, and the components were granted lower electro-osmosis voltage and higher security [11,12]. Guo et al. [13] designed some different biomimetic subsoilers based on the macroscopic surface characters of chicken clawed toes, mole crickets clawed toes, home mouse clawed toes, and so on, and they found that the biomimetic subsoilers reduced the drought force by more than 7.1%.
In view of the potential advantages of biomimetics, researchers have also developed various biomimetic potato digging shovels with excellent draught force reduction or antiadhesion performance, such as a biomimetic convex dome shaped shovel [14], mole cricket inspired shovel [15], earthworm inspired shovel [16], wild boar inspired shovel [17], etc. In fact, the macroscopic shape, microscopic morphology and cutting edge of soil engaging components have a great effect on drought force [18]. The macroscopic surface modification is a relatively simple, reliable and low-cost method for the reduction of drought force [19].
Mole rats in Gansu province who forage and nest in the potato field have two strong forelimbs and two outward front paws, which grant them brilliant soil digging abilities. The structure of the mole rat's digit is similar to that of the potato digging shovel, and the mole rat's habitat is the same as the working environment of potato digging shovels. Therefore, the mole rat's digit is a great biomimetic prototype to improve the performance of potato digging shovels by biomimetic macroscopic surface modification method.
The finite element (FE) method excels in the analysis of dynamic problems regarding material failure and large-scale deformation, which has been proven an effective approach for modeling soil-tool interaction [20,21]. ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA) has provided several different yield criteria to describe the mechanical behavior of soil, such as the Mohr-Coulomb criterion [22], Drucker-Prager criterion [23], Cam-Clay criterion [24] and their extended forms. Thereinto, the Drucker-Pager model and its extended forms are used to model frictional materials or those whose yield is associated with hardening [20,25]. The extended Drucker-Pager model has three yield surface forms, including linear form, hyperbolic form and general exponential form, of which the linear form can be calibrated with triaxial test data at different levels of confining pressures.
To our knowledge, the FE simulations of the working process of potato digging shovels imitating the mole rat's digits are scarcely reported in the literature. The objectives of the study are to (1) design a kind of mole-rat-digit-inspired biomimetic potato digging shovel (YS) with a draught force reduction property via the biomimetic macroscopic surface modification method, (2) provide a method to establish soil-shovel interaction models via ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA), (3) compare the performance of planar digging shovel (DZ) and YS under different operation conditions using the finite element models and experimental data, and (4) propose a perspective of draught force reduction mechanism for biomimetic potato digging shovels. This paper provides a reference and experience for the optimization and modification of potato digging shovels and other soil-engaging components.

Characteristics of Mole Rat's Digit
As shown in Figure 1a, the mole rat's paw is characterized by five wide and shovelshaped digits which are highly similar to potato digging shovels. Our previous studies indicate that the No. 4 digit on a mole rat's paw, located the second farthest from the mole rat's body, comes into contact with the soil most, and its ventral side is the major soil-engaging surface [26,27]. Therefore, the ventral side of the No. 4 digit is chosen as the biomimetic prototype to design biomimetic potato digging shovels. The digit of the mole rat is scanned by a non-contact 3D scanner (Vivid 910, Konica Minolta, Inc., Tokyo, Japan) and the data are exported as STL format files to build its 3D model (Figure 1b). The perpendicular bisector plane of the 3D model (Figure 1c) is captured and imported into MATLAB software to extract the coordinates of the digit's ventral side, and the fitting equation of these coordinates is obtained and expressed as Equation (1). Figure 1d shows the extracted coordinates and the fitting curve.
where x and y are the extracted coordinates and 80 ≤ x ≤ 170. The determination coefficient of the equation is 0.99, which indicates that the fitting precision is satisfied.
Appl. Sci. 2022, 12, 1761 3 of 21 conditions using the finite element models and experimental data, and (4) propose a perspective of draught force reduction mechanism for biomimetic potato digging shovels. This paper provides a reference and experience for the optimization and modification of potato digging shovels and other soil-engaging components.

Characteristics of Mole Rat's Digit
As shown in Figure 1a, the mole rat's paw is characterized by five wide and shovelshaped digits which are highly similar to potato digging shovels. Our previous studies indicate that the No. 4 digit on a mole rat's paw, located the second farthest from the mole rat's body, comes into contact with the soil most, and its ventral side is the major soilengaging surface [26,27]. Therefore, the ventral side of the No. 4 digit is chosen as the biomimetic prototype to design biomimetic potato digging shovels. The digit of the mole rat is scanned by a non-contact 3D scanner (Vivid 910, Konica Minolta, Inc., Tokyo, Japan) and the data are exported as STL format files to build its 3D model (Figure 1b). The perpendicular bisector plane of the 3D model (Figure 1c) is captured and imported into MATLAB software to extract the coordinates of the digit's ventral side, and the fitting equation of these coordinates is obtained and expressed as Equation (1). Figure 1d shows the extracted coordinates and the fitting curve. y= − 0.0002 9 ·x 3 +0.0051·x 2 +0.5645·x (1) where x and y are the extracted coordinates and 80 ≤ x ≤ 170. The determination coefficient of the equation is 0.99, which indicates that the fitting precision is satisfied.  Figure 2a shows the widely used 4U-110 potato harvester in Gansu province, China. The key components in the 4U-110 potato harvester include potato digging shovels (a-1 in Figure 2a), conveying and separating device (a-2 in Figure 2a), harvester's frame (a-3 in Figure 2a) and ground wheel (a-4 in Figure 2a). As shown in Figure 2b, the overall dimension of the original planar potato digging shovel (DZ) on the 4U-110 potato harvester is 280 mm × 140 mm × 10 mm (length × width × thickness). A potato digging shovel can be divided into three parts (P-a, P-b and P-c) based on their functions: P-a penetrates and cuts soil, P-b conveys the soil backward, which affects the drought force most, and P-c guarantees the firm attachment with the harvester's frame, while (d) and (e) refer to the boundary of the three parts mentioned above, as shown in Figure 2b,c. The difference between the original plane potato digging shovel (DZ) and the mole rat's digit inspired biomimetic shovel (YS) is that the planar structure of DZ's P-b is replaced by a spatial surface, taking Equation (1) as a collimation line ((f) in Figure 2c), while P-a and P-c remain the same.  The key components in the 4U-110 potato harvester include potato digging shovels (a-1 in Figure 2a), conveying and separating device (a-2 in Figure 2a), harvester's frame (a-3 in Figure 2a) and ground wheel (a-4 in Figure 2a). As shown in Figure 2b, the overall dimension of the original planar potato digging shovel (DZ) on the 4U-110 potato harvester is 280 mm × 140 mm × 10 mm (length × width × thickness). A potato digging shovel can be divided into three parts (P-a, P-b and P-c) based on their functions: P-a penetrates and cuts soil, P-b conveys the soil backward, which affects the drought force most, and P-c guarantees the firm attachment with the harvester's frame, while (d) and (e) refer to the boundary of the three parts mentioned above, as shown in Figure 2b,c. The difference between the original plane potato digging shovel (DZ) and the mole rat's digit inspired biomimetic shovel (YS) is that the planar structure of DZ's P-b is replaced by a spatial surface, taking Equation (1) as a collimation line ((f) in Figure 2c), while P-a and P-c remain the same.

Soil Material Model
In this paper, soil is considered an elastic-plastic continuum that reveals material hardening and the soil is modeled by the linear form of the extended Drucker-Prager model. The yield function [20,25,28] in ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA) is expressed as Equation (2).
where F is the yield function, t is the deviatoric stress given by Equation (3), p is the normal stress given by Equation (4), β and d are the internal friction angle and cohesion of the linear form of extended Drucker-Prager model, and β has a relationship with the internal friction angle φ in the Mohr-Coulomb criterion given by Equation (7).
where σ 1 , σ 2 and σ 3 are compressive stress in a triaxial test, r is the third invariant of deviatoric stress, and q is the mises equivalent stress. k is defined as the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression, which controls the dependence of the yield surface on the value of the intermediate principal stress. To ensure that the yield surface remains convex, it requires that 0.778 ≤ k ≤ 1. In this paper, soil is considered an elastic-plastic continuum that reveals material hardening and the soil is modeled by the linear form of the extended Drucker-Prager model. The yield function [20,25,28] in ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA) is expressed as Equation (2).
where F is the yield function, t is the deviatoric stress given by Equation (3), p is the normal stress given by Equation (4), β and d are the internal friction angle and cohesion of the linear form of extended Drucker-Prager model, and β has a relationship with the internal friction angle ϕ in the Mohr-Coulomb criterion given by Equation (7).
where σ 1 , σ 2 and σ 3 are compressive stress in a triaxial test, r is the third invariant of deviatoric stress, and q is the mises equivalent stress. k is defined as the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression, which controls the dependence of the yield surface on the value of the intermediate principal stress. To ensure that the yield surface remains convex, it requires that 0.778 ≤ k ≤ 1.
The soil used in the simulations is sampled at depths of 0.1-0.3 m from the soil bin in the agricultural machinery laboratory at Jilin University, Jilin Province, China. The Young's modulus of soil E and Poisson's ratio of soil µ are measured by triaxial compression tests under the confining pressure of 100 kPa, which is close to atmospheric pressure [29]. The internal friction angle ϕ of Mohr-Coulomb is calculated on the basis of data obtained in triaxial compression tests under the confining pressure of 50 kPa, 100 kPa and 150 kPa, and then the flow stress ratio k and the internal friction angle of extended Drucker-Prager model β are calculated by Equations (6) and (7), respectively. The dilation angle ψ determines the direction of the plastic strain gradient with respect to the yield surface in the p-t plane and thus introduces different plastic flow rules, such as the non-dilation flow rule that requires ψ = 0 [28], the non-associated flow rule that requires 0 < ψ < β [30] and the associated flow rule that requires ψ = β [20,31,32], which has a significant effect on the prediction of soil strength and element deformation. Based on our previous studies [33], the non-associated flow rule at ψ = 35 • is adopted in this paper. Furthermore, shear damage is used to simulate the fracture of soil in this study [32], and the required fracture strain ε s is set as 0.15, while the damage evolution type is set as displacement. The displacement at failure, which is related to element size and generally adopted in a range smaller than the element size, is set as 9 mm. The material of potato digging shovel is defined as 65 Mn steel [34], a kind of widely used material for tillage tools, and its property parameters are listed in Table 1. In addition, the default values in ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA) are adopted for other parameters.

Setup of the FE Simulations
The 3D model of potato digging shovels and the soil box are created and then imported into ABAQUS software (Dassault Systèmes Simulia Corp., Providence, RI, USA) to set up their FE models. For the convenience of modeling, the potato digging shovel and the connecting frame are immerged into a single part. The potato digging shovel is set as a rigid body with a reference point (RP) located at its tip. Shovels are meshed with C3D4 elements (4-node linear tetrahedron element) and the element size is 10 mm, while C3D8 elements (8-node linear brick element) are employed for soil boxes.
The soil box in the FE model for the interaction between soil and a single shovel is a cuboid with a dimension of 1500 mm × 1000 mm × 600 mm (length × width × height), as shown in Figure 3a,b. The dimension of the shank' cross section is shown in Figure 3c. To balance the accuracy and calculation time in FE simulations, a finer mesh with an element size of 10 mm is used for the domain in the middle of the top layer, while the element size for the remainder in the soil box is 50 mm, as shown in Figure 3a,b. The bottom left-side wall and right-side wall in the soil box are totally fixed, while other surfaces are not constrained. Meanwhile, the potato digging shovels are controlled to move only along the X direction with different speeds. In order to reach the needed speed gradually to avoid the divergence in calculations, the speed is multiplied by an amplitude function. The contact between shovel's rigid surfaces and soil nodes is defined as a surface-to-surface contact using a penalty function method, and the rigid surfaces are selected as master surfaces and the deformable soil nodes in the finer mesh domain as slave surfaces. The normal behavior between contacted surfaces is defined as hard contact, while the tangential behavior is defined as the penalty friction formulation with a friction coefficient of 0.42.
Appl. Sci. 2022, 12, 1761 6 of 21 wall and right-side wall in the soil box are totally fixed, while other surfaces are not constrained. Meanwhile, the potato digging shovels are controlled to move only along the X direction with different speeds. In order to reach the needed speed gradually to avoid the divergence in calculations, the speed is multiplied by an amplitude function. The contact between shovel's rigid surfaces and soil nodes is defined as a surface-to-surface contact using a penalty function method, and the rigid surfaces are selected as master surfaces and the deformable soil nodes in the finer mesh domain as slave surfaces. The normal behavior between contacted surfaces is defined as hard contact, while the tangential behavior is defined as the penalty friction formulation with a friction coefficient of 0.42.  Figure 4 shows the FE models for the interactions between soil and grouped shovels, which are employed to predict the field working performance of shovels. According to the requirements of ridge planting for potatoes in Gansu province, a ridge (red line in Figure 4a,b) is added on the soil surface (green line in Figure 4a,b). The top width, bottom width and height of the ridge are 650 mm, 900 mm and 150 mm, respectively. Thus, the soil box in the FE models for the interactions between the soil and grouped shovels is a cuboid (length × width × height = 1000 mm × 1300 mm × 450 mm) with a ridge on its top, as shown in Figure 4a,b. Similarly, the whole soil box is partitioned into two layers by a partition plane (black line in Figure 4a,b), and a finer mesh with an element size of 10 mm is used for the top layer, while the element size for the bottom layer is 50 mm. Moreover, the interaction parameters used in the FE models for the interactions between soil and grouped shovels are kept identical with those of the FE models for the interactions between soil and a single shovel. The bottom of the soil box is totally fixed, while other surfaces are not constrained. In addition, the potato digging shovels are also set to move only along X direction. Figure 4d,e show the finite element mesh of the soil box being cut with grouped DZ and YS, respectively.  Figure 4 shows the FE models for the interactions between soil and grouped shovels, which are employed to predict the field working performance of shovels. According to the requirements of ridge planting for potatoes in Gansu province, a ridge (red line in Figure 4a,b) is added on the soil surface (green line in Figure 4a,b). The top width, bottom width and height of the ridge are 650 mm, 900 mm and 150 mm, respectively. Thus, the soil box in the FE models for the interactions between the soil and grouped shovels is a cuboid (length × width × height = 1000 mm × 1300 mm × 450 mm) with a ridge on its top, as shown in Figure 4a,b. Similarly, the whole soil box is partitioned into two layers by a partition plane (black line in Figure 4a,b), and a finer mesh with an element size of 10 mm is used for the top layer, while the element size for the bottom layer is 50 mm. Moreover, the interaction parameters used in the FE models for the interactions between soil and grouped shovels are kept identical with those of the FE models for the interactions between soil and a single shovel. The bottom of the soil box is totally fixed, while other surfaces are not constrained. In addition, the potato digging shovels are also set to move only along X direction. Figure 4d,e show the finite element mesh of the soil box being cut with grouped DZ and YS, respectively.
The horizontal force parallel to the cutting direction is the key component of the draught force [35]; thus, the force parallel to the X direction is regarded as the draught force in this paper. As shown in Figure 5, after being compacted and cut by shovels, the displaced soil is moved forward, sideward and upward, and finally a crescent soil failure zone is formed in front of the shovels. The soil rupture distance or soil rupture distance ratio (defined as the ratio of the rupture distance r to the working depth D) is always selected to characterize the feature of this zone [36][37][38]. Thus, the draught force (F D ) under the steady state is selected as the primary indicator, and the soil rupture distance ratio (f ) as the secondary indicator to evaluate the performance of potato digging shovels in this study. Here, the soil rupture distance is measured by the border between the rupture zone and non-disturbance zone in the displacement contour plot (as shown in Figure 5b). Figure 5c shows the soil rupture observed in soil bin experiments.

Simulations of the Interactions between Soil and Potato Digging Shovels
The response surface methodology (RSM) is efficient at establishing the model between response values and selected factors as well as illustrating the optimal solution to the problem, which has been used to solve soil tillage problems [39,40]. Accordingly, in the simulations of the interaction between soil and single shovels, the effects of the cutting speed, digging depth and mounting angle, which affect the working performance of potato digging shovels significantly, on the draught force are studied by the Box-Behnken RSM, with the three factors and three levels as listed in Table 2. The levels of these three factors are selected according to the design parameters of the 4U-110 potato harvesters. Then, the soil cutting process of grouped DZ and grouped YS is simulated and compared at cutting speed of 0.56 m/s, digging depth of 0.25 m and mounting angle of 45 • based on the design parameters of the 4U-110 potato harvester.  The horizontal force parallel to the cutting direction is the key component of the draught force [35]; thus, the force parallel to the X direction is regarded as the draught force in this paper. As shown in Figure 5, after being compacted and cut by shovels, the displaced soil is moved forward, sideward and upward, and finally a crescent soil failure zone is formed in front of the shovels. The soil rupture distance or soil rupture distance ratio (defined as the ratio of the rupture distance r to the working depth D) is always selected to characterize the feature of this zone [36][37][38]. Thus, the draught force (FD) under the steady state is selected as the primary indicator, and the soil rupture distance ratio (f) as the secondary indicator to evaluate the performance of potato digging shovels in this study. Here, the soil rupture distance is measured by the border between the rupture zone and non-disturbance zone in the displacement contour plot (as shown in Figure 5b).    The horizontal force parallel to the cutting direction is the key component of the draught force [35]; thus, the force parallel to the X direction is regarded as the draught force in this paper. As shown in Figure 5, after being compacted and cut by shovels, the displaced soil is moved forward, sideward and upward, and finally a crescent soil failure zone is formed in front of the shovels. The soil rupture distance or soil rupture distance ratio (defined as the ratio of the rupture distance r to the working depth D) is always selected to characterize the feature of this zone [36][37][38]. Thus, the draught force (FD) under the steady state is selected as the primary indicator, and the soil rupture distance ratio (f) as the secondary indicator to evaluate the performance of potato digging shovels in this study. Here, the soil rupture distance is measured by the border between the rupture zone and non-disturbance zone in the displacement contour plot (as shown in Figure 5b).  . r refers to the soil rupture distance, D refers to the working depth, ζ refers to the soil rupture angle, and α refers to the rake angle.

Simulations of the Interactions between Soil and Potato Digging Shovels
The response surface methodology (RSM) is efficient at establishing the model between response values and selected factors as well as illustrating the optimal solution to the problem, which has been used to solve soil tillage problems [39,40]. Accordingly, in the simulations of the interaction between soil and single shovels, the effects of the cutting

Soil Bin Experiments
In order to validate the FE model for the interaction between soil and a single shovel, the simulated draught forces are compared with those obtained by a soil bin tester (TCC-3, Bona Science & Technology Co. Ltd., Beijing, China) in agricultural machinery laboratory at Jilin University, Jilin Province, China. The soil bin is 40 m × 2.5 m × 1.5 m (length × width × depth, the effective travel is 30 m) and it is filled with black soil. The experimental setup corresponds to that of the simulations for the interaction between soil and a single shovel. As shown in Figure 6b, the force sensors consist of an upper pull rod sensor (BK-1-LG, China Academy of Aerospace Aerodynamics, Beijing, China) and two lower suspension pin sensors (BK-5-XG, China Academy of Aerospace Aerodynamics, Beijing, China), whose maximum range are all 15 kN. Figure 6c shows the connecting shank. By adjusting the screwing length of the double-end stud (8) in the joint bearing (7), the lower shank (10) will rotate around the rotating shaft (9) and then the mounting angle of the potato digging shovel will be changed. The cutting depth is controlled by the vertical elevator on the soil bin test vehicle, while the cutting speed is controlled by the operation panel of the vehicle. The resistance signals are recorded and exported by the data acquisition and processing system on the soil bin test vehicle. Figure 6d,e demonstrates the DZ and YS used in the soil bin experiments, respectively.

Soil Bin Experiments
In order to validate the FE model for the interaction between soil and a single shovel, th simulated draught forces are compared with those obtained by a soil bin tester (TCC-3, Bon Science & Technology Co. Ltd., Beijing, China) in agricultural machinery laboratory at Jili University, Jilin Province, China. The soil bin is 40 m × 2.5 m × 1.5 m (length × width × depth the effective travel is 30 m) and it is filled with black soil. The experimental setup correspond to that of the simulations for the interaction between soil and a single shovel. As shown i Figure 6b, the force sensors consist of an upper pull rod sensor (BK-1-LG, China Academy o Aerospace Aerodynamics, Beijing, China) and two lower suspension pin sensors (BK-5-XG China Academy of Aerospace Aerodynamics, Beijing, China), whose maximum range are a 15 kN. Figure 6c shows the connecting shank. By adjusting the screwing length of the double end stud (8) in the joint bearing (7), the lower shank (10) will rotate around the rotating sha (9) and then the mounting angle of the potato digging shovel will be changed. The cuttin depth is controlled by the vertical elevator on the soil bin test vehicle, while the cutting spee is controlled by the operation panel of the vehicle. The resistance signals are recorded an exported by the data acquisition and processing system on the soil bin test vehicle. Figure 6d demonstrates the DZ and YS used in the soil bin experiments, respectively. When an experiment begins, the soil in the soil bin is first treated by rotary tillage, repression and water spray, in succession. Then, the volumetric water content (VWC) of soil is measured by a soil moisture meter (Field Scout TM TDR 300, Spectrum Technologies, Inc., Plainfield, IL, USA), and the cone index of soil is measured by a soil compaction meter with a 1 2 " diameter cone tip (Field Scout TM SC 900, Spectrum Technologies, Inc., Plainfield, IL, USA). These two indicators within the depth of 0.30 m need to be controlled in the range of 30-40% and 400-500 kPa, respectively. In addition, the experimental values of the draught force are acquired when the cutting process reaches a steady state, and each experiment is repeated three times.

Preparations and Soil Conditions
Field experiments for grouped shovels include the field draught force test and the field fuel consumption test, and they are conducted to further validate the FE model of the soil-shovel interaction and test the performance of potato digging shovels in Gansu province, China. The distance between two ridges is 1.2 m, and the dimensions of the cross section of the ridge correspond to those of the ridge in Figure 4. The field experiments are carried out under the same operation condition as that of the simulations for grouped shovels after the potato stalks have been removed. Before the field experiments start, the soil cone indices for eight sampling points are measured by a soil compaction meter with a 1 2 " diameter cone tip (Field Scout TM SC 900, Spectrum Technologies, Inc., Plainfield, IL, USA) as shown in Figure 7. gies,Inc., Plainfield, IL, USA), and the cone index of soil is measured by a soil compaction meter with a ½" diameter cone tip (Field Scout TM SC 900, Spectrum Technologies, Inc., Plainfield, IL, USA). These two indicators within the depth of 0.30 m need to be controlled in the range of 30-40% and 400-500 kPa, respectively. In addition, the experimental values of the draught force are acquired when the cutting process reaches a steady state, and each experiment is repeated three times.

Preparations and Soil Conditions
Field experiments for grouped shovels include the field draught force test and the field fuel consumption test, and they are conducted to further validate the FE model of the soil-shovel interaction and test the performance of potato digging shovels in Gansu province, China. The distance between two ridges is 1.2 m, and the dimensions of the cross section of the ridge correspond to those of the ridge in Figure 4. The field experiments are carried out under the same operation condition as that of the simulations for grouped shovels after the potato stalks have been removed. Before the field experiments start, the soil cone indices for eight sampling points are measured by a soil compaction meter with a ½" diameter cone tip (Field Scout TM SC 900, Spectrum Technologies, Inc., Plainfield, IL, USA) as shown in Figure 7.  Table 3.   Table 3. Note: The soil moisture content is measured by oven drying method and calculated as VWC (%) = (the wet weight of soil sample-the dry weight of soil sample)/(the volume of soil sample × the density of water) × 100%.

Field Draught Force Test
Figure 8a demonstrates the testing system used in the field draught force test, which consists of a tractor (YFX704, Jiangsu Yueda Intelligent Agricultural Equipment Co. Ltd., Yancheng, China, Figure 8a), a 4U-110 potato harvester and a field kinetic parameters telemetry system (Bona Science & Technology Co. Ltd., Harbin, China). The maximum range of the three force sensors used in the test is 30 kN. The values of the draught force under steady state are chosen as the experimental results, and each test is repeated three times. As shown in Figure 8b, the length of the ridges in the field draught force test is 150 m. Considering the instability of the tractor's speed during the beginning and ending phase, each ridge is equally divided into three parts, and every part consists of a 10-meter-long preparation area, S 1 , a 30-meter-long sampling area, S 2 , and a 10-meter-long finishing area, S 3 .
lemetry system (Bona Science & Technology Co. Ltd., Harbin, China). The maximum range of the three force sensors used in the test is 30 kN. The values of the draught force under steady state are chosen as the experimental results, and each test is repeated three times. As shown in Figure 8b, the length of the ridges in the field draught force test is 150 m. Considering the instability of the tractor's speed during the beginning and ending phase, each ridge is equally divided into three parts, and every part consists of a 10-meterlong preparation area, S1, a 30-meter-long sampling area, S2, and a 10-meter-long finishing area, S3. Sampling area partition for field draught force experiments (b): L1, L2 and L3 represent the experimental travelling distance, S1 represents the preparation area, S2 represents the sampling area and S3 represents the finishing area.

Field Fuel Consumption Test
The field fuel consumption test using a 4U-110 potato harvester and a MF404 tractor (YTO Group Corporation, Luoyang, China) is conducted afterwards. The fuel consumption in unit distance is chosen as the indicator, and the travelling distance in each test is 500 m. Each test is also repeated three times.

Simulations vs. Experiments for Single Shovels
Based on the Box-Behnken theories, there are 12 combination trials along with five repeated trials for central points (Nos. [13][14][15][16][17], and the simulated results are listed in Table  4. It is clear that the values of FD vary from 1.65 kN to 4.44 kN while the values of f vary from 0.95 to 2.47. The statistical analyses on the draught force of DZ and YS are performed Sampling area partition for field draught force experiments (b): L 1 , L 2 and L 3 represent the experimental travelling distance, S 1 represents the preparation area, S 2 represents the sampling area and S 3 represents the finishing area.

Field Fuel Consumption Test
The field fuel consumption test using a 4U-110 potato harvester and a MF404 tractor (YTO Group Corporation, Luoyang, China) is conducted afterwards. The fuel consumption in unit distance is chosen as the indicator, and the travelling distance in each test is 500 m. Each test is also repeated three times.

Simulations vs. Experiments for Single Shovels
Based on the Box-Behnken theories, there are 12 combination trials along with five repeated trials for central points (Nos. [13][14][15][16][17], and the simulated results are listed in Table 4.  Tables 5 and 6, respectively. The F value represents the contribution indices of the source, while the p value represents the significance of the corresponding source. Since the simulated results of five central points are identical, the F value and the p value for the lack of fit are invalid.

Simulations of the Interactions between Soil and a Single DZ
The ANOVA (analysis of variance) of the simulated draught force of DZ is listed in Table 5, and the corresponding response surfaces are shown in Figure 9. The regression equation related to the response surfaces is obtained and expressed as Equation (8),  This regression equation is highly significant at p < 0.0001, and the value of its R 2 is 0.9942. The statistical analysis for DZ demonstrates the keen sensitivity of F D to all factors (i.e., cutting speed v, digging depth D and mounting angle γ (see Figure 3a)), which can also be supported by Figure 9a-c in which the increases in these three factors all result in the increase in the draught force. The need to accelerate more soil mass and overcome increased friction is the major reason for increasing the draught force caused by the increase in cutting speed and digging depth [41]. The higher the γ, the larger the projected area facing soil mass, and thus, the draught force increases, which was reported by McKyes and Ali [36], Tamas et al. [42] and Barr et al. [43]. Moreover, Figure 9 shows that the response surfaces evolve sharply along the change both in D and v, but it turns tame along the change in γ, which can also be concluded based on their F values in Table 5, and thus the contributions of the factors to the draught force follow the order D > v > γ. Additionally, the p values of the sources indicate that the effects of v, D, γ, vD, vγ and Dγ are highly significant, while the effects of the other sources (i.e., v 2 , D 2 and γ 2 ) are not. Notes: ** means highly significant (p < 0.01), * means significant (0.01 < p < 0.05).

Simulations of the Interactions between Soil and a Single DZ
The ANOVA (analysis of variance) of the simulated draught force of DZ is listed in Table 5, and the corresponding response surfaces are shown in Figure 9. The regression equation related to the response surfaces is obtained and expressed as Equation (8), This regression equation is highly significant at p < 0.0001, and the value of its R 2 is 0.9942. The statistical analysis for DZ demonstrates the keen sensitivity of FD to all factors (i.e., cutting speed v, digging depth D and mounting angle γ (see Figure 3a)), which can also be supported by Figure 9a-c in which the increases in these three factors all result in the increase in the draught force. The need to accelerate more soil mass and overcome increased friction is the major reason for increasing the draught force caused by the increase in cutting speed and digging depth [41]. The higher the γ, the larger the projected area facing soil mass, and thus, the draught force increases, which was reported by McKyes and Ali [36], Tamas et al. [42] and Barr et al. [43]. Moreover, Figure 9 shows that the response surfaces evolve sharply along the change both in D and v, but it turns tame along the change in γ, which can also be concluded based on their F values in Table 5, and thus the contributions of the factors to the draught force follow the order D > v> γ. Additionally, the p values of the sources indicate that the effects of v, D, γ, vD, vγ and Dγ are highly significant, while the effects of the other sources (i.e., v 2 , D 2 and γ 2 ) are not.

Simulations of the Interactions between Soil and a Single YS
The ANOVA in Table 6 and the response surfaces in Figure 10 indicate that the draught force of YS is also sensitive to the cutting speed v, the digging depth D and the mounting angle γ. The regression equation with the value of R 2 being 0.9610 is expressed as Equation (9), and it is also highly significant at p = 0.0004.
In addition, the contributions of these factors follow the same order with that of DZ: The p values in Table 6 indicate that the effects of v, D, γ and vγ are highly significant and the effect of γ 2 is significant, while the effects of the other sources (vD, Dγ, v 2 and D 2 ) are not.
Compared with DZ, the draught force of YS also increases as v and D increase, but the relation between the draught force of YS and γ is quite different. As it can be seen from Figure 10b, at fixed D and low v, the draught force of YS decreases but then increases as γ increases, whereas the draught force of YS at fixed D and high v decreases and tends to be stable with the increase in γ. This similar trend was also reported by Ucgul et al. [44] and Saunders et al. [45]. In addition, at fixed v (Figure 10c), regardless of whether D is high or low, the draught force of YS decreases and then increases as γ increases, but the variation trend is tamer than that at fixed D as shown in Figure 10b. Indications mentioned here demonstrate that the interaction between v and γ is worthy of note for YS, which is different to DZ and can be confirmed by Tables 5 and 6. significant and the effect of γ 2 is significant, while the effects of the other sources (vD, Dγ, v 2 and D 2 ) are not.
Compared with DZ, the draught force of YS also increases as v and D increase, but the relation between the draught force of YS and γ is quite different. As it can be seen from Figure 10b, at fixed D and low v, the draught force of YS decreases but then increases as γ increases, whereas the draught force of YS at fixed D and high v decreases and tends to be stable with the increase in γ. This similar trend was also reported by Ucgul et al. [44] and Saunders et al. [45]. In addition, at fixed v (Figure 10c), regardless of whether D is high or low, the draught force of YS decreases and then increases as γ increases, but the variation trend is tamer than that at fixed D as shown in Figure 10b. Indications mentioned here demonstrate that the interaction between v and γ is worthy of note for YS, which is different to DZ and can be confirmed by Tables 5 and 6.

Soil Bin Experiments for Single Shovels
It is obvious from Table 7 and Figure 11a that the variation trend of the draught force obtained from the soil bin experiments matches well with those in simulations, and the experimental values with a range of 1.51-3.61 kN are slightly lower than the simulated values, but the relative errors are within 15%; thus, the simulated results are reasonable and acceptable accordingly.

Soil Bin Experiments for Single Shovels
It is obvious from Table 7 and Figure 11a that the variation trend of the draught force obtained from the soil bin experiments matches well with those in simulations, and the experimental values with a range of 1.51-3.61 kN are slightly lower than the simulated values, but the relative errors are within 15%; thus, the simulated results are reasonable and acceptable accordingly. For the soil rupture distance ratio (Figure 11b), the variation trend of the experimental results obtained from the soil bin test is generally consistent with that of the simulated results. However, great errors exceeding ±40% are seen at deep digging depth D or high cutting speed v. On the one hand, it is difficult to keep the soil conditions everywhere in each experiment perfectly identical, and the soil rupture presents high randomness under these extreme operation conditions. On the other hand, those small displacements which can be easily observed in simulations are hardly measurable in practical experiments [20]. In this regard, it is necessary to improve soil preparation technics, and the soil rupture properties under more other operation conditions also deserve further evaluation.
For the soil rupture distance ratio (Figure 11b), the variation trend of the experimental results obtained from the soil bin test is generally consistent with that of the simulated results. However, great errors exceeding ±40% are seen at deep digging depth D or high cutting speed v. On the one hand, it is difficult to keep the soil conditions everywhere in each experiment perfectly identical, and the soil rupture presents high randomness under these extreme operation conditions. On the other hand, those small displacements which can be easily observed in simulations are hardly measurable in practical experiments [20]. In this regard, it is necessary to improve soil preparation technics, and the soil rupture properties under more other operation conditions also deserve further evaluation.  Specifically, for the combinations No. 5, No. 6, No. 9, and No. 10 of YS, there are no valid test data in soil bin experiments, and the corresponding simulated values for the draught force are higher than those of DZ. For these circumstances, they can be explained based on the report of Hettiaratchi et al. [46], as shown in Figure 12. For YS, there is an angle of 23.25 • between the plane of P-a and P-c. The soil mass applies resistance onto the obverse side when γ is larger than 23.25 • (e.g., combination No. 8), which ensures that YS works at a steady depth like DZ. However, when γ is smaller than 23.25 • (e.g., combination No. 6), the rake angle α of YS exceeds 90 • , and then a newly formed soil wedge below the reverse side of YS produces a new interface, which applies large resistance forward and upward to the reverse side of YS. In simulations, potato digging shovels are constrained to move only along forward direction, which guarantees YS a steady working depth. However, in soil bin experiments, the three-point linkage of the soil bin test vehicle is floating in some range, and thus, YS is unable to work in the set depth and even cannot stick into soil, especially at a shallow depth or high speed. Therefore, it is essential to adopt γ more than 23.25 • for YS to avoid the formation of a soil wedge.
Compared with DZ, it also can be concluded from Tables 4 and 7 and Figure 11 that YS increases the draught force significantly at γ = 16 • but decreases the draught force at γ = 30 • (especially for combination No.12 where the draught force is reduced by more than 8% in both simulations and soil bin experiments), while the draught force for γ = 23 • is not significantly different, which means that YS is suitable to be adopted on the 4U-110 potato harvester at a deep digging depth and great mounting angle.
Soil failure can be defined as the permanent deformation of the soil, and thus, plastic strain can be used to evaluate the soil failure [31]. It can be concluded from Figure 13 that the great equivalent plastic strain is obvious both around the tip of shovels and in front of the shank, but the great equivalent plastic strain area of DZ is evidently larger than that of YS. Furthermore, the soil rupture angle of YS (66 • ) is 4 • lower than that of DZ (62 • ). All of the evidence means that the soil influenced by YS is more likely to move upward. Obviously, the resistance in the upward direction is much lower than the reaction force from undisturbed soil in the forward direction.

Simulations vs. Field Experiments for the Grouped Shovels
As is shown in Figure 14, the soil failure evolution under the influence of potato digging shovels (PDS) is a continuous cycle and can be summarized as follows: (1) Firstly, after the P-a of PDS sticks into and cuts soil, elastic deformation occurs and then evolves into plastic deformation [47], which spurs the emergence of horizontal cracks (Figure 14a).
(2) Next, as PDS keeps moving forward, the horizontal crack extends, the soil mass is compressed and lifted upward and backward, and then an oblique crack emerges (Figure 14b). Here, the compaction to soil accounts mainly for the increase in the draught force [48].
(3) Finally, as a critical intensity of the combination of the normal and shear stresses is reached, the oblique crack extends to the soil surface and a failure plane is formed [48,49], which leads to the rupture of soil mass and its separation with undisturbed soil (Figure 14c). Subsequently, the draught force decreases and the above process repeats, and thus the signal fluctuates. Appl. Sci. 2022, 12, 1761 Specifically, for the combinations No. 5, No. 6, No. 9, and No. 10 of YS, there are no valid test data in soil bin experiments, and the corresponding simulated values for the draught force are higher than those of DZ. For these circumstances, they can be explained based on the report of Hettiaratchi et al. [46], as shown in Figure 12. For YS, there is an angle of 23.25° between the plane of P-a and P-c. The soil mass applies resistance onto the obverse side when γ is larger than 23.25° (e.g., combination No. 8), which ensures that YS works at a steady depth like DZ. However, when γ is smaller than 23.25° (e.g., combination No. 6), the rake angle α of YS exceeds 90°, and then a newly formed soil wedge below the reverse side of YS produces a new interface, which applies large resistance forward and upward to the reverse side of YS. In simulations, potato digging shovels are constrained to move only along forward direction, which guarantees YS a steady working depth. However, in soil bin experiments, the three-point linkage of the soil bin test vehicle is floating in some range, and thus, YS is unable to work in the set depth and even cannot stick into soil, especially at a shallow depth or high speed. Therefore, it is essential to adopt γ more than 23.25° for YS to avoid the formation of a soil wedge. Compared with DZ, it also can be concluded from Tables 4 and 7 and Figure 11 that YS increases the draught force significantly at γ = 16° but decreases the draught force at γ = 30° (especially for combination No.12 where the draught force is reduced by more than 8% in both simulations and soil bin experiments), while the draught force for γ = 23° is not significantly different, which means that YS is suitable to be adopted on the 4U-110 potato harvester at a deep digging depth and great mounting angle.  Figure 15 demonstrates the simulated and field measured draught force as well as the soil failure properties of grouped shovels. It is very clear that the draught forces for both DZ and YS begin to stabilize at a displacement of about 140 mm but decrease to 0 gradually at a displacement of about 560 mm. Meanwhile, the signals of the draught force are fluctuating constantly, which corroborate the reports by Kasisira and du Plessis [48], Wang et al. [47] and Kobayakawa et al. [50]. the great equivalent plastic strain is obvious both around the tip of shovels and in front o the shank, but the great equivalent plastic strain area of DZ is evidently larger than tha of YS. Furthermore, the soil rupture angle of YS (66°) is 4° lower than that of DZ (62°). Al of the evidence means that the soil influenced by YS is more likely to move upward. Ob viously, the resistance in the upward direction is much lower than the reaction force from undisturbed soil in the forward direction.

Simulations vs. Field Experiments for the Grouped Shovels
As is shown in Figure 14, the soil failure evolution under the influence of potato dig ging shovels (PDS) is a continuous cycle and can be summarized as follows: (1) Firstly after the P-a of PDS sticks into and cuts soil, elastic deformation occurs and then evolve into plastic deformation [47], which spurs the emergence of horizontal cracks (Figure 14a) (2) Next, as PDS keeps moving forward, the horizontal crack extends, the soil mass i compressed and lifted upward and backward, and then an oblique crack emerges (Figur 14b). Here, the compaction to soil accounts mainly for the increase in the draught forc [48]. (3) Finally, as a critical intensity of the combination of the normal and shear stresse is reached, the oblique crack extends to the soil surface and a failure plane is formed [48,49], which leads to the rupture of soil mass and its separation with undisturbed soi (Figure 14c). Subsequently, the draught force decreases and the above process repeats and thus the signal fluctuates.

Simulations vs. Field Experiments for the Grouped Shovels
As is shown in Figure 14, the soil failure evolution under the influence of potato digging shovels (PDS) is a continuous cycle and can be summarized as follows: (1) Firstly, after the P-a of PDS sticks into and cuts soil, elastic deformation occurs and then evolves into plastic deformation [47], which spurs the emergence of horizontal cracks (Figure 14a).
(2) Next, as PDS keeps moving forward, the horizontal crack extends, the soil mass is compressed and lifted upward and backward, and then an oblique crack emerges ( Figure  14b). Here, the compaction to soil accounts mainly for the increase in the draught force [48]. (3) Finally, as a critical intensity of the combination of the normal and shear stresses is reached, the oblique crack extends to the soil surface and a failure plane is formed [48,49], which leads to the rupture of soil mass and its separation with undisturbed soil (Figure 14c). Subsequently, the draught force decreases and the above process repeats, and thus the signal fluctuates.  In Figure 15, there is a clear strain concentration band near the tip of shovels, and the strain near the shovel surface and the soil surface is higher than that in the inner soil. In addition, the great strain area of DZ is significantly larger than that of YS, and the thickness of the furrow slice produced by DZ (about 143 mm) is also smaller than that of YS (about 169 mm). The results mentioned here indicate that the grouped YS also reduces the compaction to soil and thus reduces the draught force. In addition, the biomimetic spatial surface of YS is able to change the stress state and motion state of the lifted soil mass continually, and this kind of change is conducive to reduce the friction and adhesion between soil-shovel interfaces [13], which improves soil mobility and makes the soil mass easier to be broken; this is also an important resistance reduction reason for biomimetic spatial surfaces. Moreover, compared with single YS, the grouped YS shows much more significant draught force reduction performance without the shank's interaction with the soil.
As shown in Figure 16, the simulated and field measured draught forces of the grouped YS are 6.40 kN and 9.00 kN, respectively, while those of the grouped DZ are 7.99 kN and 10.39 kN, respectively. It is noticeable that the draught force of YS in simulations and field experiments are 19.85% and 13.33% lower than those of DZ, respectively. In addition, the fuel consumption of YS is 9.18% lower than that of DZ, which further verifies the simulated and field measured draught forces for the grouped shovels. Therefore, YS exhibits significant performance in energy consumption reduction. Figure 15 demonstrates the simulated and field measured draught force as well as the soil failure properties of grouped shovels. It is very clear that the draught forces for both DZ and YS begin to stabilize at a displacement of about 140 mm but decrease to 0 gradually at a displacement of about 560 mm. Meanwhile, the signals of the draught force are fluctuating constantly, which corroborate the reports by Kasisira and du Plessis [48], Wang et al. [47] and Kobayakawa et al. [50].  In Figure 15, there is a clear strain concentration band near the tip of shovels, and the strain near the shovel surface and the soil surface is higher than that in the inner soil. In addition, the great strain area of DZ is significantly larger than that of YS, and the thickness of the furrow slice produced by DZ (about 143 mm) is also smaller than that of YS (about 169 mm). The results mentioned here indicate that the grouped YS also reduces the compaction to soil and thus reduces the draught force. In addition, the biomimetic spatial surface of YS is able to change the stress state and motion state of the lifted soil mass continually, and this kind of change is conducive to reduce the friction and adhesion between soil-shovel interfaces [13], which improves soil mobility and makes the soil mass easier to be broken; this is also an important resistance reduction reason for biomimetic spatial surfaces. Moreover, compared with single YS, the grouped YS shows much more significant draught force reduction performance without the shank's interaction with the soil.
As shown in Figure 16, the simulated and field measured draught forces of the grouped YS are 6.40 kN and 9.00 kN, respectively, while those of the grouped DZ are 7.99 kN and 10.39 kN, respectively. It is noticeable that the draught force of YS in simulations and field experiments are 19.85% and 13.33% lower than those of DZ, respectively. In addition, the fuel consumption of YS is 9.18% lower than that of DZ, which further verifies the simulated and field measured draught forces for the grouped shovels. Therefore, YS exhibits significant performance in energy consumption reduction.

Conclusions
In this paper, the mole rat's digits are taken as the biomimetic prototype to design a biomimetic potato digging shovel (YS), and finite element analysis, soil bin experiments as well as field experiments for the planar digging shovel (DZ) and YS are conducted to compare its performance with that of DZ. Results confirm the feasibility and effectiveness

Conclusions
In this paper, the mole rat's digits are taken as the biomimetic prototype to design a biomimetic potato digging shovel (YS), and finite element analysis, soil bin experiments as well as field experiments for the planar digging shovel (DZ) and YS are conducted to compare its performance with that of DZ. Results confirm the feasibility and effectiveness of the biomimetic macroscopic surface modification method to reduce the draught force of potato digging shovels, and the following conclusions can be drawn.
(1) The finite element soil-shovel interaction model predicts the draught force well, but slightly higher errors are seen for the soil rupture distance ratio. Hence, more efforts should be made to improve the soil preparation technics and further calibrate the parameters used in the soil-shovel interaction model. All factors have a great influence on the draught force of DZ and YS, and the order of the contribution is as follows: digging depth > cutting speed > mounting angle. For DZ, the increase in cutting speed, digging depth and mounting angle all result in the increase in the draught force; however, for YS, the draught force increases as the cutting speed and digging depth increase, but decreases as the mounting angle increases.
(2) YS, compared with DZ, illustrates brilliant draught force reduction performance. For single shovels, YS increases the draught force at low mounting angle but decreases the draught force significantly by more than 8.41% when the mounting angle is higher than 30 • . For grouped shovels, YS decreases the draught force and fuel consumption by more than 13.33% and 9.18%, respectively.
(3) The draught force reduction reasons of YS are that the biomimetic spatial surface of YS leads the soil mass to move upward and changes the soil's stress state with the motion state continually, which reduces the friction and adhesion between soil-shovel interfaces, improving soil mobility, and finally reducing soil compaction and making soil easier to be broken.
This paper provides a reference for the application of finite element method in soil tillage and draught force reduction design for tillage tools. Considering the limitation of this study, further work will be conducted to explore the effect of the soil types and conditions (e.g., soil moisture, and soil cone index) as well as other constitutive model parameters on the interaction between soil and single shovels; in addition, different operation conditions on the interaction between soil and grouped shovels will also be further studied.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.