Simulation of Track-Soft Soil Interactions Using a Discrete Element Method

: With the development of unmanned tracked vehicles, soil model predictions of soft terrains are becoming more essential. In order to accurately simulate the interaction characteristics between soil particles and the track, soil modeling with a discrete element method (DEM) is proposed. Volume-based scaled-up modeling of DEM soil particles and the calibration of DEM input parameters were investigated as a feasible approach to realizing many particle calculations. Calibration of DEM input parameters can solve the distortions between actual and DEM particle sizes. Cohesion and friction parameters of the scaled-up soil particle model were recalibrated by the shape accumulated through the virtual design of the experiment. Soil DEM particles were scaled up to 1 mm spherical particles, and recalibrated DEM parameter values were used to match the actual accumulated soil shape. Three calibrated scaled-up soil models were used for the shear stress–displacement DEM simulation of a track segment, and the mean absolute percentage error ( MAPE ) was less than 11% compared with the actual shear stress–displacement test. The parameter value of soil traction performance empirical model of a tracked vehicle is modified according to the soil shear stress–displacement DEM simulation. Comparative analysis was performed for travel test results of a tracked vehicle; the relative error of the soil traction prediction results to actuals was less than 16.8%. This showed that the volume-based particle scaling technique is an effective DEM for the mechanical simulation of soil.


Introduction
With the development of special unmanned tracked vehicles, more unmanned tracked vehicles are being used in soft soil environments, i.e., river dredging tracked equipment.In order to estimate the trafficability of soft soil, the research focus of the terramechanics of tracked vehicles has gradually shifted to the detailed modeling of soft soil.
In soft ground conditions, the tractive performance of a tracked vehicle is determined by the interaction between the track and the soil.Previous researchers studied the effect of the geometry of the track segment on tractive performance [1,2], such as grouser height, pitch, thickness and pitch-to-height ratio.Most of the contemporary research has focused on land soil.For soft soil with a high moisture content, the traditional method is still used to test the mechanical properties of the soil.Due to the high sinkage of soft soil, traditional soil shear tests are unsuitable.
Most previous studies on the interaction between soil and track have been performed through a soil bin in the laboratory or on an actual field.With the development of computer technology, the application of DEM to interaction problems in terramechanics has become popular.DEM was first proposed and introduced into the study of vehicle terramechanics [3].Nakashima used a two-dimensional DEM to simulate soil, and verified the effect of track traction in combination with a soil bin test [4].Du assessed the mechanical properties of soft ground with low moisture, low viscosity and high friction [5].Peters established a vehicle-terrain simulation model based on the DEM and proved the accuracy of traditional drawbar-pull tests [6].Zhou proposed a systematic and accurate discrete element modeling method for sandy roads [7].The verification results showed that the modeled interactions between track and sandy ground based on DEM simulations were consistent with the results of the actual soil bin test.For fine soil particles, a major disadvantage of DEM simulation in practical applications is the challenge of modeling very small particles; it is computationally impossible and impractical to calculate every particle of true size.
In order to reduce the amount of calculation, a feasible method is to apply scaled-up particle technology to DEM simulation.At present, the scaled-up particle method is more commonly used in DEM simulation [8][9][10].If all particle sizes are scaled up by a constant factor, but the material properties remain the same, the dynamics of the granular system will be altered [11].Thus, DEM parameter calibrations are the key challenge for the scaled-up particle method.
Many contemporary scholars have focused on the parameter calibration method of DEM simulation.DEM calibration consists of the careful choice of the material input parameters needed to reflect the real behavior of the particles in the numerical simulation.Careful calibration is the basis for achieving predictive results.The methods of parameter calibration include the angle of repose (AOR) test, penetration test, compression test and shear test [12][13][14][15][16].Due to the large number of material attributes to be set, most of the calibration is focused on selected parameters.More calibration procedures for the DEM parameters have been developed to reduce computational effort [17][18][19].
In the studies reported above, the DEM was applied to the terramechanics of tracked vehicles; the main research objects are limited to land soil and sand.However, the applicability of DEM simulation to soft soil is uncertain.The objectives of this study were: (1) to calibrate scaled-up particle contact parameters consistent with the mechanical properties of soft soil; (2) to validate the shear properties of the scaled-up particle model and (3) to corroborate the soil traction performance of the tracked vehicle according to the modified soil shear DEM simulation.

JKR Model in the DEM
The Johnson-Kendall-Roberts (JKR) contact model was developed from the Hertz-Mindlin contact model.The Hertz-Mindlin contact model is the basis of DEMs, and is accurate and efficient.Here, the Hertz contact model is applied for the calculation of normal forces [20], and the Mindlin contact model is applied for the calculation of tangential forces [21].As a supplement, the damping component should be considered based on the damping coefficient and recovery coefficient [22].Additionally, tangential friction and rolling friction should complied with the Mohr-Coulomb law and the contact-independent directional constant torque model [23].
Based on the theory of the Hertz-Mindlin model, the load analysis for interactions among particles is derived as five categories: normal elastic force, normal damping force, tangential elastic force, tangential damping force and rolling friction.For soft soil with a high moisture content, the liquid bridge force has a great influence on the interaction between particles, which is ignored by the Hertz-Mindlin model.In the JKR model, the normal elastic force is developed and modified, with a comprehensive analysis of cohesive forces caused by the action of water.The JKR model is a method of contact simulation for cohesive particles, which can be used to simulate the cohesive effects among wet particles, as shown in Figure 1  In order to characterize the cohesion between particles, a new parameter of surface energy is introduced to calculate the cohesive forces.The equivalent JKR normal elastic force can be expressed as [24]: where * E is the equivalent elastic modulus, * R is the equivalent radius, γ is the sur- face energy, δ is the overlap and α is the interaction parameter.
The equivalent elastic modulus and equivalent radius could be calculated with Hertz [20] model.

(
) ( ) where i E and j E are the elastic moduli of two contacted particles; i υ and j υ are the Poisson ratios of two contact particles and i R and j R are the radii of two contacted par- ticles.
In the JKR contact model, even if the particles are not in contact, cohesion forces still occur between particles.The maximum clearance that produces non-zero cohesion force between particles can be expressed as [24]:

Calibration of the JKR Model Parameters
Riverbed surface soil is a complex sediment with a high moisture content.With the increase in moisture content, the soil presents a plastic state, fluid plastic state and flow state.The different content of sand in the soil will also change the soil characteristics.In this paper, the clay and sand were collected from the Huai River, a representative river in northern China, used as the in situ data for tests and simulations.As a supplement, three kinds of clay and sand mixed soil, 1# (75:25 wt.%), 2# (50:50 wt.%) and 3# (25:75 wt.%), were applied to soil shear tests and DEM simulations.The wet densities of mixed soil 1#, 2# and 3# were 1.54 g cm −3 , 1.65 g cm −3 and 1.74 g cm −3 , respectively.
The discretization of DEM particle sizes and their distribution which comprises particle sizes of natural soils is very complex and computationally prohibitive.The main particle size composition of soil can be screened by standard assays.Based on the main particle size of the soil, soil particles are assumed to be uniformly distributed.Under the assumed condition, the scaled-up particle is set to 1 mm in the DEM simulation, which can simplify and accelerate the calculation [25].
In this paper, standard screening and drying experiments were applied to test the physical properties of soil.First, the standard screening of 300 mesh and 200 mesh was used to identify clay particles of sizes 0.05-0.075mm.Meanwhile, the standard screening of 35 mesh and 18 mesh was used to identify soil particles of sizes 0.5-1 mm.Second, the collected soil was placed in a graphite crucible and warmed in the incubator at 100 °C.The weight of the crucible was measured every 24 h until the error of two measurements was less than 1%.The wet density, dry density and moisture content (water mass fraction) of the soil could be obtained through calculations, with the results shown in Table 1.As shown in Table 1, particle sizes of clay and sand were less than 0.075 mm and 1 mm, respectively, and modeling soil based on actual particle sizes requires a significant amount of computing resources.The optimal method is to increase the size of soil particles proportionally and derive the same mechanical properties as the actual soil particles by setting appropriate parameters.In this study, for comprehensive consideration of the actual particle size and the calculation ability of computer, the clay and sand particles were set as spherical particles 1 mm in size in the DEM simulation [26].The characteristic coefficients of particle interactions were determined by comparing the real experiment with the simulation [15,16].
Three types of parameters were indispensable in the JKR contact model [14]: (1) The intrinsic parameters of particles, including wet density, Poisson ratio, shear modulus and particle size, are basic parameters and can be measured by the test.
The parameters are shown in Table 2; (2) The contact parameters of particles, including the static friction coefficient and rolling friction coefficient, are friction characteristic parameters and can be measured by tests or calibrated by simulations.The friction between particles is calculated according to the friction characteristic parameters, and the specific calculation method refers to the Hertz-Mindlin model [20,21]; (3) The JKR models specific particle parameters, i.e., surface energy, which is a cohesive characteristic parameter and is difficult to measure.However, it can be calibrated by DEM simulation.The cohesive force between particles is calculated according to the cohesive characteristic parameters, as shown in Equation ( 1).For soil with a high moisture content, the repose angle of soil particles is affected by the friction and cohesion of particles.Thus, it is necessary to recalibrate the characteristic parameters of scaled-up particles.
In this paper, the lifting cylinder test was used to calibrate the soil repose angle.First, the cylinder filling soil particles was placed on the surface of the test platform and raised vertically at a speed of 5 mm s −1 by a screw device driven by a stepper motor (57BL55S06).The soil particles fell freely through and formed a conical pile on the base plate.Second, the same experimental parameters were set in the DEM simulation.The basic contact parameters of particles are shown in Table 2. Finally, the contact parameters that needed to be calibrated were finally confirmed by comparing the accumulation test and simulation.The dimensions of the cylinder and the cylinder model established in the EDEM software with the same size are shown in Figure 2. As shown in Figure 3e, the repose angle of sand particles is mainly determined by friction and can be measured directly.Unlike the formation principle of the repose angle of sand particles, clay is a kind of soil with a small internal friction angle.The repose angle of clay particles is mainly determined by the cohesion, which causes the irregular shape accumulation and is difficult to measure accurately.
To accurately measure the characteristic parameters of mixed soil and analyze the accumulation of mixed soil with different proportions of sand, a method of image edge extraction and pixel superposition through a saliency detection network and grey-binary processing was proposed with MATLAB software.A schematic diagram of the soil contact parameter calibration method is shown in Figure 4.

Shear Stress-Displacement Model
The locomotion of a tracked vehicle on the surface of soil is mainly based on the shear forces that develop between the track links and surface soil in the longitudinal direction.Soil conditions have significant effects on tractive performance.
According to prior surveys and tests, the characteristics of clay are similar, in that shear stress initially increases sharply and reaches a "hump" of maximum shear stress max τ at a particular shear displacement, and then decreases and approaches a relative constant residual value res τ with further increases in shear displacement, which rela- tionship may be characterized by the following empirical equation proposed by Wong [32,33]: ) where max τ is the maximum shear stress, w K is the shear displacement where the maximum shear stress occurs and r K is the ratio of residual shear stress to maximum shear stress.However, sand has shear stress characteristics which increase rapidly and then approach a relatively constant value, for which the relationship may be characterized by the following empirical equation proposed by Janosi [34]: where K is the shear displacement where the maximum shear stress occurs.Variation trends in shear stress of different types of soil are shown in Figure 5.

Shear displacement (mm)
Shear characteristic curve of dense clay (Wong model) Shear characteristic curve of loose sand (Janosi model) Shear stress (kPa) Figure 5. Schematic of the shear characteristic curve of the track segment.The red curve is derived from Wong's [32,33] model.The blue curve is derived from Janosi's model [34].

Shear Stress-Displacement Test
The shear test apparatus consisted of an electrical step motor which produced a constant horizontal displacement in the soil bin (600 mm × 300 mm × 300 mm) with a soil height of 30 mm and a resistance strain force sensor, as shown in Figure 6.The force sensor (DYMH-106) is a strain sensor with a comprehensive error of less than 0.3% F.S. To ensure the track segment shear stability process, the motor speed was set to 0.05 m/s and the sampling frequency was set to 50 Hz.The test platform was connected to a computer for data acquisition and processing.It should be noted that after each stressdisplacement test, the soil and water were thoroughly mixed and pre-compressed according to the grounding contact pressure of the tracked vehicle.Schematic diagrams of the experimental platform and track segment.Parameters of the track segment structural dimensions are shown in Table 3.

Support frame Soil
Step motor module Screw

Shear Stress-Displacement DEM Simulation and Modification Process
According to the soil shear stress-displacement test, the virtual soil bin in the DEM simulation was set to 130 mm × 120 mm × 22 mm containing 400,000 soil particles, considering the computer's capability.The computer's CPU was an Intel 7700 processer with four cores and eight threads, and the simulation time was 42 h.
The virtual soil bin and the shear simulation process between the soil and track segment with DEM simulation are shown in Figures 8 and 9, respectively.The number of soil particles that produced the bulldozing effect above the track segment was about 165,000, and the soil volume is 85.3 cm 3 .The movement speed of soil particles on the track segment was about 0.036 m s −1 , which was less than the movement speed of the track segment 0.05 m s −1 .The difference in the velocity among soil particles led to the bulldozing resistance of the track segment.In order to remove the bulldozing resistance caused by soil accumulation, the calculation domain of the particles was adjusted.The soil shear test was performed based on adjusting the particle calculation domain.The virtual soil bin and simulation of the shear process between the clay and track segment with the EDEM software are shown in Figure 10.In the virtual shear simulation, there was no soil accumulation in front of the track segment.The image of the motion state of clay particles was derived, as shown in Figure 10.The motion of soil particles near the inner wall of the soil bin was almost static, which indicated that soil particles near the inner wall did not participate in the shear process, and the size of the soil bin was suitable.The motion state of clay particles participating in the shear simulation was about 0.045 m s −1 , which is consistent with the moving speed of the track segment-0.05m s −1 .According to the motion of the particles in front of the track segment, there were about 23,500 particles, and the effective soil volume involved in the shear process was 24 cm 3 .Compared with the soil shear DEM simulation without modification, the volume of accumulated soil causing the error was 60 cm 3 .This shows that the error caused by soil accumulation could not be ignored.
The mean absolute percentage error (MAPE) was applied to evaluate the accuracy of the simulation model.When the MAPE value was close to zero, the DEM simulation model was more consistent with reality.
where s n is the number of samples, i y is the test value of samples and ^i y is the sim- ulation value of samples.

Traction Calculation and Travel Test of Tracked Vehicle
The shear stress at any point of contact can be determined based on the hypothesis that contact pressure is distributed uniformly.With Equation ( 7), the traction force t F can be deduced by the integral [32,33]: 2 where A is the contact area between track and soil b is the track width, L is the contact length between track and soil, h is the grouser height, W is the gravity of tracked vehicle and i is the slip.In order to verify the accuracy of the DEM and evaluate the soil tractive performance, external load estimated by the motor current was carried out.Firstly, the corresponding relationship between motor current and load torque was analyzed based on the motor characteristics test data.Secondly, the soil tractive performance could be calculated according to the load torque of the motor and the structure of the tracked vehicle.A schematic diagram of the tracked vehicle is shown in Figure 11.

Drive Motor Characteristic Tests
A motor characteristics test platform was designed, as shown in Figure 13.The test platform was equipped with a high-precision DC brushless servo motor, a planetary gear reducer, a torque sensor, and a magnetic powder brake.The specifics of the motor test platform are shown in Table 5.

Servo brushless DC motor
Planetary Reducer  For a comprehensive consideration of traveling speed and capacity of the tracked vehicle, the sampling frequency of the motors was set at 10 Hz, and the motor power was set at 60 W. The relationship between the motor current and external load torque was obtained at a series of certain rotational speeds.The relationship between the motor current and external load torque is shown in Figure 14.According to the test data, the relationships between the torque, speed and current of the motor are considered comprehensively, and the load torque T can be expressed as:

Magnetic powder brake
where T is the load torque of the motor, I is the current of the motor and n is the rotary speed of the motor.

Mechanical Property Verification of Scaled-Up Particles
The accumulated shape, obtained by pixel superposition, is shown in Figure 15.With the increase in sand particle content, the accumulation height of soil decreases gradually, and the flow state becomes more obvious.When the content of sand particles is low, the interaction between soil particles is mainly a cohesive force.However, with the increase in the content of sand particles, pore water film between the clay particles is broken.This leads to the molecular force between clay particles becoming smaller, and soil cohesion also decreases.Moreover, when the soil is mainly composed of sand particles, the pores between sand particles are filled with clay particles.These high-moisture-content clay particles affect lubrication, which causes the soil to become a flow state.As analyzed above, the AOR of clay particles and sand particles is different.The simulation of soil accumulation can reproduce the accumulation test results with the properly adjusted friction and cohesion parameters of soil.Friction and cohesion characteristics of different soils are shown in Figure 16.The above results indicate that both cohesive and frictional properties exist in clay particles.However, friction is the main characteristic of sand.Meanwhile, the friction characteristic curves of soils with different compositions are nonlinear.Clay and sand are modeled as spherical particles, which means that the friction parameters between particles are set as larger than that of the actual situation.The significance of accumulation is to explore appropriate contact parameters according to different soil compositions.
It should be noted that the friction and cohesion characteristics in this study were based on a soil particle size of 1 mm.Additionally, the calibration of friction and cohesion are greatly affected by particle size.The coordinate ordinate axes represent the friction coefficient and surface energy of soil particles based on spherical scaled-up particles.

Comparative Analysis of the Results of the Soil Shear Test and DEM Simulation
As the shear displacement of the track segment increases, more soil is accumulated in front of the track segment (marked by the red circle in Figure 7).Additionally, the test data collected by the force sensor include the shear force between the track segment and the soil and the resistance caused by the soil accumulation in front of the track segment.The shear stress-displacement test result has great errors without completely removing the influence of "bulldozing".However, it is impossible to completely remove the accumulated soil in front of the track segment in the test.As shown in Figure 17, the shear stress values of clay test data increase continuously.The accumulated soil in front of the track segment produced a great impact in the measurement and was difficult to remove completely.However, the accumulated soil in front of the track segment could be effectively removed by changing the calculation domain in DEM, which could modify the shear relationship model of the soil and track segment accurately.
The MAPE of clay, sand and mixed soil were 7.5%, 11% and 9%, which indicated that the prediction model is accurate and effective.Shear stress-displacement fitting data from the test and simulation are shown in Figure 18 (a) (b) (c)

Modified DEM Simulation of Track Segment Shear Process
The stress along the shear direction of the track segment was converted into a line graph, as shown in Figure 19.Through the comparative analysis of Figures 18 and 19, it can be seen that the bulldozing resistance caused by soil accumulation in front of the soil plate was effectively eliminated.The force variation trend was more consistent with the shear regularity of the clay and track segment.The stress trend in the track segment following the shear regularity could also be obtained from the sand and mixed soil simulations.The key parameters in the classical model are summarized in Figure 19 and detailed in Table 6.20; the stress variation in the DEM simulation is similar to that of the classical model.Meanwhile, the maximum shear stress and residual shear stress in the DEM simulation are consistent with the classical model.According to the shear stress-displacement data and fitting curve, the necessary parameters of the shear stress-displacement model could be obtained, as shown in Table 6.However, as the maximum shear stress changed to the residual shear stress, errors were obvious between the simulation and classical model.The pores between the scaled-up particles in the DEM simulation were larger than the actual pores of soil particles, which mean that the particles in the DEM simulation are easier to compress and rebound.In the terramechanics of the tracked vehicle, maximum shear stress, residual shear stress and location of maximum stress determined the accuracy of soil traction calculation.Shown in Figure 20, the values of the maximum shear stress, residual shear stress and location of the maximum stress fitted well in the DEM simulation and classical model.Errors between the simulation and classical model will not affect the derivation of soil traction of the track.

Soft Soil Traction Verification
In order to reduce the influence of noise, the test data were applied with a moving average filter using the MATLAB software, as shown in Figure 21.The rotary encoder obtained the rotation speed signals of two motors, as shown in Figure 22.In the tracked vehicle travel tests, clay was difficult to mix uniformly, which caused the rotation speed signal to fluctuate at some positions.In order to eliminate influences of motion state change to current measured values, the acceleration and deceleration processes were excluded when analyzing the test data.The valid interval of test data was selected as from 5 to 60 s.The external load was calculated according to the characteristic relationship formula of the valid interval, as shown in Table 7.According to the characteristic of the motor, the motor current was converted to external load torque, as shown in Figure 23.The external resistance forces estimated by the motor current data and Wong model (Equation ( 9)) were 130.2 N and 108.3 N, respectively.The relative error between the soil traction performance calculation based on the scaled-up particle DEM and actual value was 16.8%, which is suitable for engineering applications.

Figure 2 .Figure 3 .
Figure 2. Bottomless cylinder in accumulation test and simulation.Trial-and-error methods are feasible to change the parameters and comparison results constantly until the simulation results conform to the actual experimental results.Different accumulations were measured through the above experimental steps.Accumulation tests and DEM simulations were carried out, as shown in Figure3

Figure 4 .
Figure 4. Flow chart of parameters calibration in test and simulation.

Figure 8 .
Figure 8.The virtual model of the soil bin and track segment.(a) DEM soil bin model with scaled-up particles.(b) Virtual track segment model (80 mm length, 60 mm width and 12 mm grouser height).

Figure 9 .
Figure 9. DEM simulation of soil shear test.

Figure 10 .
Figure 10.Shear process of track segment and clay in the DEM simulation.

Table 4 .
Specifics of the tracked vehicle.tests were performed on a clay terrain, as shown in Figure12.In the tracked vehicle travel test, the rotation speed of both motors was set to 50 rpm, and the tracked vehicle speed was 0.028 m s −1 .Limited by the soil bin, the tracked vehicle's total distance and test time were about 1.7 m and 60 s, respectively.

Figure 12 .
Figure 12.Tracked vehicle travel test on soft soil.

Figure 15 .
Figure 15.Accumulations of tests and simulations with five types of soil.The edge curves of soil accumulation and DEM particle accumulation shapes were extracted.Two edge curves are compared by pixel superposition.(a) Pixel superposition of sand.(b) Pixel superposition of clay.(c) Pixel superposition of mixed soil 1#.(d) Pixel superposition of mixed soil 2#.(e) Pixel superposition of mixed soil 3#.

Figure 16 .
Figure 16.Friction and cohesion characteristics of different soils.The cohesion and friction characteristics of scaled-up particles were calibrated by DEM simulation.The two ordinates are the cohesion and friction of scaled-up soil particles.The abscissa is the soil composition.

Figure 17 .
Figure 17.Shear test data and curves of three types of soil.

Figure 18 .
Figure 18.Curve fitting of shear stress-displacement in the test and simulation.(a) Trends in the test and simulation for clay.(b) Trends in the test and simulation for sand.(c) Trends in the test and simulation for mixed soil 2#.

Figure 20 .
Figure 20.Curve of simulation data and classic model.(a) Curve fitting of clay.(b) Curve fitting of sand.(c) Curve fitting of mixed soil 2#.

Velocity of particle i and j m s − 1 i ω , j ω Angular velocity of particle i and j rad s − 1 nFn
Elastic modulus of particle i and jPa i υ , j υPoisson ratio of particle i and j i R , j RRadius of particle i and j mRotary speed of the motor rpm

Table 1 .
The characteristics of clay and sand.

Table 2 .
DEM parameters used in the simulations.

Table 5 .
Specifics of the motor test platform.

Table 6 .
Parameters of the shear stress-displacement model.

Table 7 .
Parameters of motor characteristics.