Understanding of penetrating resistance on an inclined plate to cohesive soil using two-dimensional discrete element method

In order to know the effect of varied grousers on the reaction of soil, the reaction of soil by a plate at different angles of inclination using the two-dimensional discrete element method (DEM) has been investigated. Either the simulation or the experiment, a flat iron plate with 5 × 77 × 100 (mm) was used, and the inclination angle between the plate and the soil surface ranged from 0 to 60 degrees with 15-degree intervals. According to the result, the cohesive soil model introduced in the 2D DEM could successfully replicate the experimental results considering the effects of different inclination angles of the plate and the values of the root mean square error between the simulation and the experiment were found as 6.62, 36.60, 58.63, 60.84, and 66.25 N for inclination angles from 0 to 60 degrees in the horizontal direction. Moreover, the RMSE values of 127.91, 216.06, 75.88, 110.40, and 25.49 N were obtained for the vertical direction. Overall, a minimal discrepancy was observed between the results of DEM simulation and those of experiments in the vertical component of reaction force when the inclination angle of the plate was increased.


Introduction
Soil reaction to running gears or tools contacting with terrain has been one of popular research subjects in terramechanics. Tractive performance should be improved for off-road vehicles, and, similarly, the resistance force on a soil-working tools needs to be reduced from a viewpoint of economical utilization of fossil fuels (Dechao & Yusu, 1992;Naderi-Boldaji et al., 2014;Yang et al., 2014;Johnson et al., 2015;Ge et al., 2016;Jiang et al., 2018).
When we focus on an interaction between soil and grouser, the soil reaction on the grouser plate has been studied by applying a semi-empirical, or parametric, approach (Bekker, 1956;Bekker, 1960;Bekker, 1969;Wong, 2001;Tiwari et al., 2010). The effect of link pitch for rigid link tracks was analyzed, and the longer link pitch caused an increase in the pull coefficient (Gao & Wong, 1993;Wong, 2010). The most effective gross traction was reported to be obtained with a ratio of grouser pitch to height in the range of 3-4 (Hata & Hosoi, 1981). Muro et al. (1988) confirmed experimentally that the ratio of grouser pitch to height to be about 3.2 to obtain maximum tractive effort in cohesive soil by using model track. By applying the parametric approach, a bounding limit of grouser height for maximum gross tractive effort was also reported (Wang et al., 2002;Wong & Huang, 2006). In addition, the effect of inclination angle on loads has also been studied in other research areas, such as Ou Scholars investigated the responses of inclined loaded piles in layered foundations, this paper proposes a generalized solution based on the principle of Minimum Potential Energy, the study is different from the object of this paper, and the research idea is the same (Ou et al., 2022).
With recent development of computer technology (Zang and Zhao, 2013), the discrete element method (DEM) proposed by Cundall & Strack (1979) has become popular for studies on contact interactions in terramechanics to fully consider the discontinuous nature of soil (Smith & Peng, 2013). For soilgrouser interactions, Oida et al. (1997) first analyzed the interaction between a single grouser shoe and cohesive soil by DEM. By adjusting the density of the element and the spring constant for contact model, they could obtain a result of gross traction comparable to their experimentally obtained results. Asaf et al. (2006), Asaf et al., 2007) applied a commercially available DEM code, PFC2D, to the analysis of link-track performance, and reported that the maximum shear stress showed linearly increasing behavior with grouser height up to 26 mm. By applying a quasi-2D experiment where a soil model of consisting of aluminum cylinders of two radii, Nakashima et al. (2015) showed that their in-house 2D DEM analysis could reproduce the experimentally obtained results of gross tractive effort. Moreover, the effect of open spaces between grousers on the gross traction of a track shoe was investigated by 2D DEM (Hettiaratchi et al., 1966;Yokoyama et al., 2020).
In general, a typical vertical grouser on the track shoe has often been studied previously as stated above (Keen et al., 2013;Zhang et al., 2018). There is also other type of grouser, called triangle grouser, used especially for weak soil condition. The prediction of gross traction for this type of grouser has not yet been studied sufficiently.
Since the oblique face of triangle grouser has two roles---one for supporting contact load and the other for generating gross traction, the observation of these vertical penetration and horizontal translation using an inclined plate should bring an insight into the interaction between the oblique face of triangle grouser and soil. As the first attempt for the numerical observation on an interaction of the triangle grouser, DEM may be applied to vertical penetration and then to horizontal translation of inclined plate. Vertical penetration of inclined plate becomes similar to plate penetration test in parametric approach when the angle of inclination is set to 0 deg (Yong & Hanna, 1977). Moreover, the horizontal translation of inclined plane resembles to the cutting of soil by flat tool when the rake angle is set smaller or equal to 90 deg (Hermawan et al., 1996;Hermawan et al., 1997;Watyotha et al., 2001). As recently studied by Murino Kobayakawa, the relationships between different rake angles and drag force changes of plates is explored using DEM from the certain point of the variation of the local volume fraction inside the shear band. And then, the reason for the observed change in drag force is explained by using a three-dimensional wedge model (Kobayakawa et al., 2020), Many previous reports could be found for such rake angles for soil-tool interactions, such as buckets (Cleary, 1998;Kuczewski & Piotrowska, 1998;Coetzee et al., 2007), plane tool (Momozu et al., 2003;Tanaka et al., 2007;Ono et al., 2013;Jiang et al., 2017;Kobayakawa et al., 2018), bulldozer blade Tsuji et al., 2012), and subsoiler implement (Chen et al., 2013;Wang et al., 2019a).
The properties of soil are always affected by the consisted elements (Coetzee, 2016;Coetzee, 2017). Moreover, Ono et al. (2013) investigated the effects of elemental shape on the cutting resistance of soil by a narrow flat blade using 3D DEM. Wang et al. (2019b) also investigated that the varying particle radii affected soil subsoiler interactions in the DEM models and the 7-mm-radius particles were recommended as the best choice for the Hertz-Mindlin contact with a bonding model.
The purpose of this study is to investigate the effect of the inclination angle of a flat plate on the soil reaction in the plate penetration and translation processes using 2D DEM with a cohesive soil model. Since the computational load is still high for 3D DEM, the 2D DEM was applied in this study (Nishiyama et al., 2016). For calibration of DEM parameters and comparison of accuracy, previous experimental data on plate penetration has been used.

Plate penetration experiments
A flat iron plate measuring 5 mm × 77 mm × 100 mm (thickness, length, and width, respectively) was used for the experiment. The plate was set at different angles of inclination

Frontiers in Built Environment
frontiersin.org and the vertical and horizontal soil reaction were detected using the extended octagonal ring transducer. A schematic diagram of the device is shown in Figure 1. As shown in the figure, the plate was driven by a hydraulic cylinder, the force and displacement of which were detected by an extended octagonal force sensor (stress concentration type) and a displacement transducer, respectively. The device was also used for validation experiments. The soil was sandy loam, which was prepared from the experimental farm field of Mie University. A direct shear test was carried out to measure the internal and external coefficient of friction, adhesive and cohesive strength. Table 1 lists the parameters of the test soil. The plate experiment was only carried out once for every inclined plate condition.

Cohesive soil model
The soil model was first established and calibrated to give a solid foundation to the experiment (Du et al., 2017;Chi & Kushwaha, 1990;Coetzee, 2019). To calibrate the soil model, we carried out a series of tentative penetration tests. The banding and algorithm model of the compound particles of soil clumps were decided at the beginning.
In general, in DEM, the contact forces between non-cohesive solid particles can be expressed as the Voigt model, as shown in Figure 2, where C n and K n are the coefficients of damping and elasticity in the normal direction, respectively, and C s and K s are the coefficients of damping and elasticity in the tangential direction, respectively. The local coordinate system is defined as in Figure 3. A global Cartesian coordinate system x-y and a relatively moving local coordinate system u s -u n have been established. Based on the Voigt model and the coordinate system, the contact forces exerted on the particle i could be expressed as Eqs 1, 2: where F n represents the normal force, u n represents the displacement and v n represents the velocity in the normal direction. Similarly, the tangential contact force can be expressed as follows: where F s represents the tangential force, μ represents the friction coefficient, f n represents the normal force, v s represents the tangential velocity component, Δx 1 represents the previous displacement of particle i and Δx represents the variation quantity in the tangential direction. There is also a rolling resistance moment M, which could be represented by Eq. 3: where α is the coefficient of resistance moment, and b is half the contact length (Sakaguchi, 1995). However, the target soil in the present study is a cohesive one, and cohesive soil model should be introduced, where soil particles may be clumped together by liquid bridges, as in Figure 4A. In this

FIGURE 3
Definition of the local coordinate system in DEM.

Frontiers in Built Environment
frontiersin.org study, a simple model using a cohesive force f c has been added (Bui, Kobayashi, Fukagawa, Wells, 2009) (see Figure 4B).
The cohesive force f c can be acquired through the coefficient c ′ n . and the contact area of two particles (Bui, Kobayashi, Fukagawa, Wells, 2009). Thus, the normal and maximum tangential forces acting on particle i in Figure 4B are expressed by the following equations, respectively: in which φ is the internal friction angle between the particles. The motions of an arbitrary element should be governed by Newton's second law for the DEM algorithm. If there is an element d, the specific calculation of the motions can be expressed by Eq. 6:

Analysis of the movement and forces of the plate
The plate used for penetration tests had a dimension of 5 mm × 77 mm × 100 mm (thickness, length and width, respectively). A schematic diagram of the plate movement with its inclination angle is shown in Figure 5.
As shown in Figure 5, the depth of penetration in the simulated tests corresponds to previous experimental tests that used sandy clay. The forces and moment of a single particle can be calculated with Eqs 3-5. If there are m particles in contact with the plate, the resultant force enacted on the plate can be expressed as follows: where F x p and F y p represent the force components acted on the plate in the x and y directions, respectively; F x i and F y i represent the force components of a single contacting particle with the plate in the x and y directions, respectively. DEM simulation was repeated three times for a given inclination condition.

Calibration of DEM parameters
2.3.1 Behaviour of the soil model Figure 6 shows the result of the collapse simulation test for the non-cohesive soil model and the cohesive soil model. The initial configurations of the soil models were both at a scale of 270 mm × 500 mm (height × length), and cells had the dimensions of 100 mm × 100 mm. The cohesive force coefficient c ′ n of the soil was 3.58 kPa. As shown in the figure, the soils were marked by different colours in order to readily observe the failure surface of the soil. The soil particles stuck to each other, as shown in Figure 6B, thus a fragment forming on the collapsed surface was observed. The collapsed distance of the cohesive soil model became shorter. This is a result of the effect of the interlocking force between the soil particles, which was caused by the introduction of the cohesive soil model.

Calibration of the soil parameters
To calibrate the DEM parameters, series of penetration tests were deployed through a plate with a 0°inclination angle. The device and instruments were as same as those which were introduced in Section 2.1.

Data processing in DEM
In this study, the dimensions of the numerical soil bin were set at a length of 0.5 m, a width of 0.3 m and a height of 0.9 m.

FIGURE 8
A series of trial runs for determining the parameters of the inhouse program.

FIGURE 9
Initial configuration of the soil and the plate model in DEM calculation.

Frontiers in Built Environment frontiersin.org
Three different radii of particles, 3.0, 2.0, and 1.2 mm, were put into the soil bin at a ratio of 1:2:3 according to the number of particles. The total number of particles for each DEM was approximately 9400. Initially, the particles were distributed randomly in the soil bin with no overlaps. The time step could be estimated using Eq. 8 and was 1 × 10 −5 s in this study. Because of the existence of gravity, particles went down freely and piled up, as shown in Figure 7.
Δt < 0.385 m min K n where Δt represents the timestep, m min represents the minimum mass of element. Frontiers in Built Environment frontiersin.org 06 As Figure 7 shows, the soil particles were uniformly and randomly distributed in the soil bin at the initial stage. At 0.3 s, the packing procedure was completed when the soil surface was at a height of approximately 270 mm, and the bulk density became approximately 2,011.0 kg·m −3 . According to the global coordinate established in Figure 5, the reaction force of soil on the plate can be decomposed into two directions: the horizontal one, f x , and the vertical one, f y .

Result of calibration
As shown in Figure 8, gravity acceleration of 1 G (−9.81 m·s −2 ) was employed in DEM runs 2, 6, 7 and 8, and an increased gravity of 2 G was exerted on the particles in DEM runs 1, 3, 4 and 5. Meanwhile, the normal stiffness was changed from 2.5 × 10 4 to 5.0 × 10 4 kPa. According to Figure 8, a series of trial tests were carried out using the established 2D DEM model. Supplementary  Table S2 summarises the root-mean-square error (RMSE) of the results between the laboratory experiment and the DEM simulation. From this table, the moduli of the soil in DEM run 2 was determined for DEM, where the normal stiffness was 1.0 × 10 5 N·m −1 and the tangential stiffness was 2.5 × 10 4 N·m −1 , whereas the density of the model soil was set to 2,500 kg·m −3 and the coefficient of friction was set to 0.8.

Soil reaction on inclined plate
As in Figure 9, the soil particles were settled in the soil bin. Then, the plate penetrated the soil at a speed of 0.01 m·s −1 . The behaviour of the soil particles and the plate at the 70-mm penetration depth at different inclined angles are depicted in Figures 10A−E. The figure shows that the colours of the soil particles are related to its angular velocity: the redder the colour, the higher the angular velocity. Figures 12-15, show comparisons of f x and f y between the simulation and experimental results.  In these figures, the simulation is named as Sim.01, Sim.02, and Sim.03. To compare the results of the experiment and the simulation, the average value of the three simulations has been calculated at each equivalent depth of the experiment.
As shown in Figure 11, the simulation well replicated the experiment in all plate reaction forces before the penetration depth of 40 mm when the inclination angle of the plate was 0°.
For the horizontal force fx, the results of the experiment and the simulation were all close to nil N when the plate was sinking into the soil. In the experiment, the vertical force fy proportionally increased with the increasing depth of penetration, and the maximum value reached 680.9 N at a depth of 78 mm.
In comparison, the vertical force of the simulation increased more slowly with the increase in the depth of penetration, and the average maximum force was 393.7 N at a depth of 72.4 mm. For the plate with 0°angle of inclination, the RMSE of fx was 4.57 N and that of fy was 85.0 N. Figure 12 shows the results for the plate with a 15°i nclination angle.
In the case of the plate with 15°inclination, the maximum depth of penetration in the experiment and simulation was 120 mm. As can be seen in Figure 12. At a depth of 120 mm, the minimum fx was −13.9 N for the experiment and −48.9 N at a 94.7 mm depth for the average DEM result. The RMSE of the fx component was 25.88 N.
In the vertical direction, all force components f y lie in a positive quadrant, The f y components in the experiment and the simulation increased along with an increased depth of penetration. At a depth of 120 mm, the maximum value was 885.2 N for the experiment, and at a depth of 119.5 mm, 483.1 N for the average simulation result. The RMSE of the fy component was 216.06 N.
The result for the plate with 30°angle is shown in Figure 13. As the figure shows, the maximum penetration depth was 145 mm for either the experiment or the simulation.
For the x-axis force component, the results of both the experiment and the simulation decreased with an increased plate depth. At a depth of 145 mm, the minimum f x in the experiment was −18.9 N, while the minimum in the simulation reached −111.2 N at 140.9 mm. Almost all the f x values in the simulation were less than those of the

FIGURE 13
Reaction against penetration depth to a 30°inclined plate in cohesive soil. (A) Horizontal reaction against penetration depth to a 30°inclined plate in cohesive soil. (B) Vertical reaction against penetration depth to a 30°inclined plate in cohesive soil.

FIGURE 14
Soil reaction varying with penetration depth by 45°inclined plate.
(A) Relationship of horizontal reaction to penetration depth when a 45°i nclined plate was utilized. (B) Relationship of vertical reaction to penetration depth when a 45°inclined plate was utilized.
Frontiers in Built Environment frontiersin.org experiment when the plate penetrated the soil. Although the plate's inclination angle was increased in 15°increments, the f x components between the plates of 15°and 30°were close to each other in the experiment. However, the simulation results reveal an obvious situation of descent on the force-depth relation when increased the penetration depth. The RMSE of the f x component was 59.95 N, which means that there is a large margin of disagreement between the results of the experiment and the simulation.
In the y-axis direction, the force component has a maximum value of 700.4 N at a 145-mm depth in the experiment and 479.0 N at a 143.4-mm depth in the average simulation result. The RMSE of the f y component was 100.25 N, implying that the results of these two methods were closer to each other than the result at the inclination angle of 15°. Nevertheless, they matched each other until 110-mm depth; then, there was an increasing deviation. Figure 14 shows the results of the reaction to a 45°angle of inclination. As the figure demonstrates, although the horizontal force component for the experiment was less than that of the 15°and 30°angles at the same depth, there were still no significant changes.
The lowest f x value of this 45°angled plate at a depth of 138 mm was −42.2 N for the experiment and −146.0 N at 136.4 mm depth for the average DEM result. Essentially, the f x component decreased with the increased penetration depth for both the experiment and the simulation.
The force component along y-axis increased proportionally with the increase of the penetration depth for either the experiment or the simulation. When the penetration depth was 138 mm, the vertical force was 661.5 N, the highest value in the experiment. Similarly, the highest value for the f y component in the average simulation result was 337.5 N at a depth of 137.6 mm. The RMSE was 163.1 N for the vertical force components. The vertical force component was replicated well in advance of the depth at 60 mm ( Figure 14B). Thus, the experimental result shows a more rapid increase than that in the simulation.
The result of the reaction to a 60°angle of inclination is shown in Figure 15. The horizontal force f x on the plate in the experiment slightly decreased when the plate was penetrating the soil; by contrast, f x in the simulation decreased rapidly with the increased depth of penetration. Specifically, the minimum value for the experiment was −32.9 N at a depth of 130 mm, whereas the minimum for the simulation was −115.8 N at a 129.1-mm depth.
In Figure 15B, the simulation tests have similar results and close values to those of the experiment when they were at the same depth. Similar to the case in the previous tests of other angles of inclination, the greater the penetration depth, the larger the value of the vertical force components for the plate inclined at 60°. The maximum value was 272.3 N at 130 mm for the experiment and 232.9 N at 130 mm for the average simulation. As the results of the simulation and experiment matched well, the RMSE was 35.7 N for the vertical force.
In sum, if the plate penetrated at an equal depth, the greater the inclination angle of the plate, the lower the amount of vertical force; nevertheless, the horizontal force was increased in the negative direction according to the simulation results. In terms of the RMSE values, the vertical force has been simulated more effectively than the horizontal force. The reason for this phenomenon is perhaps the different shapes of the soil clumps that formed beneath the plate that moved together with the plate in the simulation and experiments. In addition, the experiment was under the 3D condition, and DEM was under the 2D condition. Because of the difference in degrees of freedom between 2D and 3D, the particles' movement may cause different reactions in the vertical and horizontal directions. These differences lead to different replicated performances in the simulations of the vertical and horizontal force components. The results suggest that the 2D DEM disagreed with the experiment for f x ; meanwhile, the 2D DEM agreed well in advance of the 40-mm depth between them for f y .

Concluding remarks
In this study, a DEM soil model was calibrated using the results of an experiment on the penetration of a plate with no inclination. After this, a series of simulation tests were carried out on the

FIGURE 15
Soil reaction varying with penetration depth by 60°inclined plate. (A) Horizontal reaction against penetration depth to a 60°inclined plate in cohesive soil. (B) Vertical reaction against penetration depth to a 60°inclined plate in cohesive soil.

Frontiers in Built Environment
frontiersin.org penetration of the plate inclined at angles between 0°and 60°. Based on the global coordinate system, the reaction force on the plate has been broken down into horizontal and vertical components, and the results were compared with that of the experiment. From the comparison, the following conclusions can be drawn: (1) The 2D DEM approach could be sufficient for the demands of researching the soil reaction forces acting on a plate, which were affected by different angles of inclination. The whole process was performed successfully in this study.
(2) The cohesive soil model introduced in the 2D DEM could successfully replicate the experimental results considering the effects of different inclination angles of the plate, and from the RMSE values, the results of the soil reaction in the DEM simulation in the vertical force component were more accurate than those of the horizontal force component.
(3) Base on the limitation in this study, the greater the inclination angle of the plate, the lower the amount of vertical reaction force; whereas, horizontal force insensitively affected by the plate's movement.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions
JG, the corresponding author, is in charge of parts like data processing, and program compiling. HN, supervise the study. XW, provided the experimental data.

Funding
This research was partly supported by China Scholarship Council (CSC) and Entrepreneurship and Innovation Program of Hefei.