Modeling soil-bulldozer blade interaction using the discrete element method (DEM)

Journal of Terramechanics 88 (2020) 41–52

Contents lists available at ScienceDirect

Journal of Terramechanics journal homepage: www.elsevier.com/locate/jterra

Modeling soil-bulldozer blade interaction using the discrete element method (DEM) Mehari Z. Tekeste a,⇑, Thomas R. Way b, Zamir Syed c, Robert L. Schafer b a

Department of Agricultural and Biosystems Engineering, Iowa State University, 2331 Elings Hall, Ames, IA, USA USDA-ARS National Soil Dynamics Laboratory, Auburn, AL, USA c Department of Civil, Construction and Environmental Engineering and Department of Agricultural and Biosystems Engineering, Iowa State University, USA b

a r t i c l e

i n f o

Article history: Received 23 August 2019 Revised 3 November 2019 Accepted 25 December 2019

Keywords: Discrete element method Similitude Soil horizontal force Soil vertical force

a b s t r a c t Limited studies have been conducted to establish scaling relationships of soil reaction forces and length scales of bulldozer blades using the Discrete Element Method (DEM) technique. With a DEM-based similitude scaling law, performance of industry-scale blades can be predicted at reduced simulation efforts provided a calibrated and validated DEM soil model is developed. DEM material properties were developed to match soil cone penetration testing. The objectives of the study were to develop a DEM soil model for Norfolk sandy loam soil, establish a scaled relationship of soil reaction forces to bulldozer blade length scales (n = 0.24, n = 0.14, n = 0.10, and n = 0.05), and validate the DEM-predicted soil reaction forces on the scaled bulldozer blades to the Norfolk sandy loam soil bin data. Using 3D-scanned and reconstructed DEM soil aggregate shapes, Design of Experiment (DOE) of soil cone penetration testing was used to develop a soil model and a soil-bulldozer blade simulation. A power fit best approximated the relationship between the DEM-predicted soil horizontal forces and the bulldozer blade length scale (n) (R2 = 0.9976). DEM prediction of soil horizontal forces on the bulldozer blades explained the Norfolk sandy loam soil data with a linear regression fit (R2 = 0.9965 and slope = 0.9634). Ó 2019 ISTVS. Published by Elsevier Ltd. All rights reserved.

1. Introduction Utilizing simulation-based engineering tools to design and evaluate performance of earthmoving equipment requires modeling soil-tool interaction. As tools exert loads on soils, soil deformation could occur in a combination of volumetric soil compression, soil distortion in shear failure, and cutting (Gill and Vanden Berg, 1968). Soil reaction varies depending on the soil type and soil conditions, which strongly influence the performance of machine systems, for example the energy required for cutting and moving soils, and productivity of excavator buckets and bulldozer blades during soil filling and flow (unloading). Simulation of soil-tool interaction is thus an essential component of virtual earthmoving equipment design and analysis. The discrete element method (DEM) assumes soil as a discontinuous medium (Cundall and Strack, 1979) and can predict various dynamic soil behaviors (volumetric compression, shear, and flow) as tools interact with soil. Other mesh-based and continuum numerical techniques such as the Finite Element Method (FEM) and Computational Fluid Dynamics (CFD) have ⇑ Corresponding author at: Department of Agricultural and Biosystems Engineering (ABE) at Iowa State University, Ames, IA 20010, USA. E-mail address: [email protected] (M.Z. Tekeste). https://doi.org/10.1016/j.jterra.2019.12.003 0022-4898/Ó 2019 ISTVS. Published by Elsevier Ltd. All rights reserved.

limitations in capturing soil deformation involving cutting, mixing and flow (post-failure soil behavior) (Shmulevich et al., 2007), but these soil behaviors commonly occur in soil-blade interactions. DEM was successfully applied to simulate soil-tool interaction problems from wide cutting blades (Hofstetter, 2002, Asaf et al., 2007, Shmulevich et al., 2007; Obermayr et al., 2014) and from tillage tools (chisel plow, cultivator sweeps) (Chen et al., 2013, Ucgul et al., 2015, Tekeste et al., 2019). Most of the studies reported fairly good agreement between simulation and experiments in predicting soil forces at relatively low tool speed (Hofstetter, 2002, Shmulevich et al., 2007 and Obermayr et al., 2014). In those soil-tool studies, the DEM simulation setup did not always obey soil, geometry, and kinematic (tool velocity) similarities to their respective physical experiments. DEM soil-tool simulation procedures to reproduce the physical experimental setup (geometry, soil and kinematics) varied widely by the analyst’s experience or available computing resources. In an attempt to reproduce the physical experimental setup within the DEM simulation, analysts appear to be driven by ways to reduce computational effort, experience in selecting primitive DEM particle shapes to approximate scaled soil particle physical dimensions (shape and size), number of DEM particles relative to the soil-engaging geometry dimensions (tool width and tool cutting depth), and relative scaled tool velocity to

42

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

reduce simulation time. Hofstetter (2002) simulated soil cutting forces from a straight bulldozer blade (400 mm wide (W), 140 mm height (H) and 25 mm blade cutting depth (z)) using a three-cluster sphere DEM soil particle (10.5 mm particle diameter (d)). The experiment was conducted in a soil bin filled with ‘‘artificial soil” with a soil-internal friction angle of 38° and a cohesion value of 13 kPa. Hofstetter (2002) obeyed blade geometry (blade width and blade cutting depth) similarities between the DEM and the soil bin experiment. They, however, had a distorted soil-blade system using a 10.5 mm DEM soil particle diameter and approximately two DEM particles per cutting depth (z/d = 2.38). They used this to approximate the ‘‘artificial soil” particle (with clay particles measuring approximately less than 2 lm) and an experimental tool cutting depth of 25 mm in the soil bin. DEM simulation of soil-blade cutting of soil particles was conducted in a virtual soil bin filled with 20,000 DEM particles. Hofstetter (2002) reported good agreement between DEM and laboratory measured soil horizontal forces after the tool travelled a distance of 650 mm (L/W = 1.625). Shmulevich et al. (2007) simulated cutting dry sea sand (#40 ASTM standard sieve screened Cesarea sand) by a blade (199 mm wide and a cutting knife length of 30 mm) at two blade cutting tip depths (z) of 25 mm and 50 mm. The DEM simulation by Shmulevich et al. (2007) was conducted using two-dimensional (2D) modeling in PFC2D which is a DEM commercial code (Itasca Consulting Group, Inc., Minneapolis, MN). A total of 15,000 DEM particles were used and each consisted of a two-disk clumped primitive particle (particle disk radius of 0.7–1.1 mm (3.6 mm maximum length) and 0.64 mm disk radius distance between the two disks). They were filled in a virtual soil bin (450 mm long and 130 mm depth). The simulation was run cutting the DEM soil by the blade traveling at 20 mm/s (10 times the laboratory test speed, to reduce simulation time). The DEM predicted soil horizontal force from the blade (W/z = 7.96 for the 25 mm cutting depth) after traveling a distance, L, of 200 mm (L/W = 1.005) showed good agreement with the laboratory-measured values. Obermayr et al. (2014) also used an assembly of DEM particles to simulate soil cutting forces from a blade using three blade widths (W) (50 mm, 100 mm, and 200 mm) and three blade cutting depths (z) (20 mm, 40 mm, and 60 mm). The experiment was conducted in a 300 mm by 300 mm square box filled with a rounded gravel material. Each DEM particle was 1.5 times as large as the equivalent diameter (d ranges from 1.5 to 3 mm) of the rounded gravel. After a 50 mm length (L) travel distance (L/W = 2) using the 100 mm wide (W) blade and 40 mm cutting depth (z), the DEM predicted soil horizontal force from the cutting blade matched fairly well to the laboratory measured data (within one root mean square deviation). Similar to Hofstetter’s (2002) DEM formulation steps, Obermayr et al. (2014) and Shmulevich et al. (2007) did not attempt to create DEM particle counts the same as their experimental granular materials. Hofstetter (2002), Obermayr et al. (2014), and Shmulevich et al. (2007) each used different scaling approaches. Calibration of DEM parameters to match bulk soil properties from direct shear tests in Hofstetter (2002) or simple wedged-tool penetration tests in Asaf et al. (2007) might have helped to obtain good agreement with the experiment data of soil horizontal forces on blades. Even though their numerical analyses had limitations matching soil vertical forces and soil rupture, with the DEM formulations employed in these studies, the relationships between soil reaction forces and tool design parameters (blade width and cutting depth) were successfully established and used for earthmoving equipment design (e.g. Hofstetter, 2002). DEM particle definition, material properties, and the simplification techniques used for the relatively small blade sizes (e.g. 199 mm, 300 mm, and 400 mm wide (W) in the above referenced works) may not necessarily be scaled easily to evaluate perfor-

mance of industry-scale blades. DEM simulations of soil-blade interactions to establish relationships of DEM-prdicted soil reaction forces and length scales of blades have not been investigated. Similitude techniques were applied in Garrity et al. (1968) and Reaves et al. (1969) to evaluate performance of bulldozer blades from soil bin experiments and they established scaling laws between soil horizontal forces and blade length scales from n = 4 to n = 20. We hypothesize those results as a basis to investigate DEM-based scaling law relationships between soil reaction forces on a blade, and a length scale (n). The overall goal of this study was to investigate a DEM-based scaling relationship between soil reaction forces (horizontal and vertical) and length scale (n) of bulldozer blades. The approach used in this study may be adapted to enable proper scaling of performance of soil reaction forces from small-size bulldozer blades to full industry scale bulldozer blades. Specifically the objectives were (1) to develop a DEM soil model in a soil-blade interaction for cohesionless soil, (2) to investigate the relationship of soil forces to length scale of bulldozer blades, and (3) to validate the DEM-predicted scaling predictions of soil forces to soil bin experimental data. 1.1. Similitude technique review for soil-bulldozer blade analysis Reaves et al. (1969) applied similitude to predict the horizontal component of soil force on bulldozer blades from soil and length scale variables which significantly affect bulldozer blade performance. Using the dimensional analysis of a soil-blade system (Fig. 1), a dimensionless dependent variable of soil horizontal force (D/qWd) is related to the dimensionless independent variables of blade geometry, soil, and the system according to Eq. (1).

D ¼f qWd




l d H R A B C E CI xW V 2 ; ; ; ; ; ; ; ; a; b; l; ; ; W W W W W W W W q q gW

 ð1Þ

where, D = horizontal component soil force on the blade (soil horizontal) (F); q = soil strength parameter (F L2) from surcharge soil; W = blade width (L); d = blade cutting depth (L); l = distance of blade loading (L); H = blade height (L); R = radius of curvature of moldboard (L);

Fig. 1. Left side view showing bulldozer blade geometry, after Reaves et al. (1969). Drawing at right side shows detail of cutting edge.

43

472

0.05 0.10 0.14

1,166 4,165 734

Soil cone index (% d.b.)

1.63 72:17:11

As explained in Reaves et al. (1969), all values were averages from the top 150 mm soil layer. Model coefficients and predicted soil horizontal forces using power fit equation in Garrity et al. (1968). b

a

7.6

a

Soil moisture content (Mg/m3) a

Norfolk sandy loam

The model coefficients (a and b) depend on soil type and condition. In order to predict soil horizontal force on soil-bulldozer blade using Eq. (3) (D = a nb), Garrity et al. (1968) conducted a series of soil bin experiments using four geometrically similar and scaled bulldozer blades on Norfolk sandy loam soil (72% sand, 17% silt and 11% clay). The initial soil conditions for the top 150 mm of depth were a dry bulk density of 1.63 Mg/m3, 7.6% d.b. soil moisture content, and a soil cone index of 734 kPa (Reaves et al., 1969). The four geometrically similar bulldozer blades had a length scale (n) of n = 0.24, n = 0.14, n = 0.10, and n = 0.05. A full industry length scale bulldozer blade (n = 1) had a width (W) of 3.2 m and a height (H) of 1.07 m. The tests were operated at a tool speed of 0.11–0.22 m/s and a blade cutting depth of 20% of each of the scaled bulldozer blade heights (H). The power law-fitted relationship resulted in model coefficients of a = 148 kN (33,300 lb) and b = 2.5 (R2 = 0.9827) from the soil bin experiment on Norfolk sandy loam soil. The summary of results from Garrity et al. (1968) and Reaves et al. (1969) is shown in Table 1. The soil horizontal force was approximated from a bulldozer blade at full soil load.

Soil dry bulk density

D = soil horizontal force; n = blade length scale; a = model coefficient that estimates soil horizontal force from a full scale bulldozer blade (n = 1), and b = exponential model coefficient that estimates the slope of a straight line relationship between soil horizontal force and the blade length scale on a logarithmic scale.

% sand:% silt:%clay

where,

Soil properties

ð3Þ

Soil type

Reaves et al. (1969) and Garrity et al. (1968) simplified the piterm relationships between the dependent and the independent variables in Eq. (1) to predict the bulldozer blade soil horizontal force (D) from a length scale (n) according to Eq. (3) using the assumptions on distorted soil-tool systems in Schafer et al. (1969).

Table 1 Predicted soil horizontal force, D (N) using dimensional analysis on Norfolk sandy loam soil, from Garrity et al. (1968) and Reaves et al. (1969).

df = the force prediction factor; Dp = prototype force; Dm = model force; n = blade length scale (n = lp/lm). lp is the characteristic length of prototype and lm is the characteristic length of model, and s = an exponent for a particular soil-blade system.

a

(kPa)

where,

D ¼ anb

0.24 b

ð2Þ

a (kN)

Model coefficientsb

df ¼ Dp =Dm ¼ ns

2.5

Soil horizontal force, D (N)b Length scale (n)

In the soil-blade system, ensuring geometric similarity of the model and prototype, and assuming the distorted model system and negligible effect from tool velocity on bulldozer blade force, Schafer et al. (1969) showed experimental evidence from several soil-machine systems (one of them being a bulldozer blade) that the force prediction factor ratio (df) of the soil horizontal force of the prototype (Dp) to the soil horizontal force of the model (Dm) was assumed to be related to the length scale ratio (n).

148

A, B, C, E are lower tangent length, length of cutting blade, offset of cutting blade, and thickness of cutting edge, respectively (Fig. 1) (L); a = angle of curvature of moldboard; b = lift (rake) angle of cutting blade; l = soil-metal friction coefficient; CI = soil cone penetration resistance (F L2); x = soil bulk density (F L3); V = velocity of blade (L T1), and g = gravitational acceleration (L T2).

85

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

44

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

1.2. DEM approach to approximate similitude relationship of soilbulldozer blade experiment In order to investigate the soil horizontal force and length scale relationship for the soil-blade system in DEM, the geometric and kinematic similarities in Garrity et al. (1968) were used. To maintain the conditionality for the distorted model system in Schafer et al. (1969) with respect to soil properties, we proposed a calibration methodology to determine DEM material properties of Norfolk sandy loam soil, that matches mean soil cone penetration resistance. With a calibrated DEM soil model reproducing the soil cone penetration resistance in Norfolk sandy loam soil, the distorted soil model assumption in Schafer et al. (1969) can be satisfied. Similarly we assumed the power fit relationship between soil horizontal force and the length scale of bulldozer blades in Garrity et al. (1968) and Reaves et al. (1969) can be formulated from DEM simulation of soil-bulldozer blade interaction. 2. Materials and methods The simulations of soil-cone penetrometer and soil-bulldozer blades were implemented in EDEM 2018 (DEM Solutions Ltd., Edinburgh, UK), a commercial DEM code (EDEM, 2011). The simulation work flow consisted of (1) generating CAD geometry surface meshes of the ASABE cone penetrometer (cone base and shaft) (ASABE, 2018a) and scaled bulldozer blades, (2) defining a DEM soil particle shape which approximates the soil aggregates for the soil-cone penetrometer, and the soil-bulldozer blade interactions, (3) assigning material model properties for soil and tool (steel), (4) specifying DEM solver settings including a fixed time step and grid sizes, and (5) post-processing to extract soil reaction forces on the conical tip of the ASABE cone penetrometer and scaled bulldozer blades. The total number of particles to approximate full soil load was also extracted from the soil-bulldozer blade simulation. 2.1. DEM soil particle: Shape and particle size approximation Discretization of DEM particle sizes and their distribution that comprises particle sizes (sand (2 mm) to clay (less than 2 lm) range) of natural soils is very complex and computationally prohibitive. Contrary to, for example, a 4-node quadrilateral rigid or shell element which has proven to be useful in FEM analysis of soil-tool interaction problems (Tekeste et al., 2009), creating DEM particle shapes for soil modeling has not been developed. In

many previous studies (Hofstetter, 2002, Asaf et al., 2007, Shmulevich et al., 2007; Obermayr et al., 2014; Chen et al., 2013, Ucgul et al., 2015, Tekeste et al., 2019), scaled particle shapes and sizes (single-sphere or clumped spheres) are arbitrarily chosen to perform DEM calculations at a reasonable computational effort. Agricultural soils, however, are often found as soil aggregates with sub-angular, angular, granular and platy structure shapes (Hillel, 1998). A systematic shape approximation approach was proposed in our study to create approximate DEM soil particles. In this study, we hypothesized creating scaled DEM particles comprising a single-sphere or clumped spheres that approximate naturally occurring large size soil aggregates. This may provide a better alternative approach to previous studies. A 3D laser scanner was used to generate approximated CAD geometry soil aggregates. DEM clumped (overlap) spheres were generated in EDEM 2018 pre-processing to match the 3D reconstructed shapes of soil aggregates. Agricultural soil from a cultivated field site that exhibits natural soil aggregates was sampled to attempt DEM particle size and shape approximation. We used this approach instead of using arbitrarily selected scaled sizes. Soil samples from a Clarion loam soil (fine loamy, mixed, mesic Typic Haplaquolls) (33% sand, 45% silt and 22% clay) were collected after field cultivation in Boone, Iowa. The soil aggregates after being air-dried at room temperature, were mechanically screened through ASTM standard sieves (12.5 mm, 4.75 mm, and 0.85 mm screen sizes). Samples retained on the 12.5 mm, 4.75 mm, 0.85 mm, and passing through 0.85 mm, were selected for 3D scanning and reconstructing shapes for DEM particle shape approximation (Figs. 2 and 3). Note that we were not attempting to imply that the natural soil particle sizes for sand, silt and clay fractions should match between the 3D scanned soil aggregates and the soil used for validation of the cone penetration testing and soil-blade simulation studies. Five samples from the soil aggregates retained on the 12.5 mm, 4.75 mm, and 0.85 mm screens were 3D-scanned and reconstructed to characterize the aggregated shape for DEM approximations. The 3D scanner (Artec Space Spider model, Artec Studio, Luxembourg) has a resolution of 2.75e-5 mm. The soil aggregate shape representations were reconstructed in ANSYS SpaceClaim 3D modeling software (SpaceClaim Corp., Concord, Mass.). Similarly using a digital caliper with a resolution of 0.01 mm, the major and minor axis dimensions of the sampled aggregates were measured. The accuracy of the 3D scanning image capturing and reconstruction was first verified by comparing the relative error with the laboratory measurement. A relative error was found to be 6 –11%. Fig. 3 shows the 3D reconstructed soil aggregate CAD geometry obtained from the equivalent screen sizes of passing through 0.85 mm, and retained on 0.85 mm, 4.75 mm, and 12.5 mm screens, and DEM spheres created in EDEM 2018 matching the major axis length and inscribing the soil geometry using clumped spheres. 2.2. DEM soil material constitutive model

Fig. 2. ASTM standard sieves screened sample of soil aggregates (A) retained on 12.5 mm, 4.75 mm, 0.85 mm, and passing through 0.85 mm screen sizes from left to right, and (B) dimensional measurement of soil aggregates measured along the minor and major axis using a caliper with a resolution of 0.01 mm.

Hertz-Mindlin (HM) with parallel bond, available in EDEM 2018, was the contact physics law used for DEM calculation of forces and displacements from the interacting soil-soil and soil-tool (steel) elements. Tekeste et al. (2019) reported good prediction of soil draft forces on a cultivator sweep on loam soil using this HM with parallel bond model. This DEM contact model was assumed to reproduce the soil behaviors of Norfolk sandy loam soil during the insertion of the conical tip of the ASABE cone penetrometer and the cutting process by the scaled bulldozer blades. The HM with parallel bond includes Hertz-Mindlin contact theory (Walton and Braun, 1986) for calculating the Hertzian contact normal forces and Mindlin contact tangential forces, and the rolling resistance model explained in Sakaguchi et al. (1993)

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

45

Fig. 3. Clumped-sphere DEM shapes approximating soil aggregate samples passing through a 0.85 mm screen (A), retained on 0.85 mm screen (B), retained on 4.75 mm screen (C), and retained on 12.5 mm screen (D).

46

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

to calculate torque on a particle, and the parallel bond model from beam theory after Potyondy and Cundall (2004). Details of the HM with parallel bond constitutive equations and their application for soil-tool interaction simulation are available in EDEM (2011) and Tekeste et al. (2019). The HM with parallel bond contact model in EDEM 2018 requires (1) material parameter values of Poisson’s ratio, shear modulus (Pa) or Young’s modulus (Pa), particle density (kg/m3) for each material type (soil and tool (steel)); (2) interactional parameter values of coefficient of restitution, coefficient of static friction, and coefficient of rolling friction for soil-soil and soil-steel materials. Values of the parallel bond components, normal stiffness per unit area (Kn) (N/m3), shear stiffness per unit area (Kt) (N/m3), bond disk radius (mm), critical normal stress (Pa), and critical shear stress (Pa), are required.

Table 3 DEM parameters for initial soil cone penetration simulation. DEM input parameter

Value

Reference

Particle (single-sphere) Particle radius (r) (mm)

1.5

Table 1

0.32 2650 (Top layer) 2675 (Middle layer) 2675 (Bottom layer) 2.26e + 07 (Top layer) 2.26e + 07 (Middle layer) 2.26e + 07 (Bottom layer)

Literaturea Calibratedb

Soil material properties Poisson’s ratio Particle density (kg/m3)

Shear modulus (Pa)

2.3. DEM soil model calibration 2.3.1. Simulation setup for DEM parametric sensitivity analysis Ideally all DEM micro-model parameters that affect soil reaction behaviors (forces and flow) in a general soil-machine system could be calibrated. However, the computational efforts to perform calibration of all DEM model parameters are prohibitive. Measurement of the micro-model parameters at the individual soil particle size are also not possible. Researchers first consider a sensitivity analysis of limited DEM model parameters that strongly influence soil-tool responses, and then apply a design-of-experiment (DOE) technique to calibrate the DEM parameters to match bulk experimental data (Sadek and Chen, 2015, Mousaviraad et al., 2017; Tekeste et al., 2019). As part of the DEM soil model calibration approach, DEM parameters that strongly influence the soil cone penetration are identified based on soil characterization and prior experimental or numerical analyses of soil-tool interaction problems on similar sandy soils (Tekeste et al., 2007, Sadek and Chen, 2015) and loam soil (Tekeste et al., 2019). Soil reaction forces on a vertically moving cone penetrometer (Sadek and Chen, 2015) and soil reaction forces on a blade (Tekeste et al., 2019) were found to be sensitive to the bond stiffness parameter. Sadek et al. (2011) successfully calibrated the Particle Flow Code in Three Dimensions (PFC3D) parallel bond model bond stiffness parameter (kn) reproducing soil yield forces from virtual direct shear tests for varying soil conditions (dry to wet soil moisture levels and loose to compacted bulk density levels) of sandy soil (86% sand, 4% silt and 10% clay). In this study, we adopted the Sadek et al. (2011) approach to perform a sensitivity study of the bond stiffness (kn) DEM parameter and to calibrate a DEM soil model to match soil cone penetration data for a Norfolk sandy loam experiment, for dry and wet soil moisture conditions. For all the other parameters, the DEM determination was conducted one parameter at a time using prior experimental and FEM soil constitutive model parameter effects on soil cone penetration, in the Norfolk sandy loam soil. Particle density values were adjusted during particle assembly generation to match laboratory-measured soil bulk density of the three layers (top, middle and bottom) of Norfolk sandy loam. Values for soil

Table 2 Soil bulk density and soil moisture content (mean (standard deviation)) at the three soil layers (top, middle (hardpan) and bottom) of Norfolk sandy loam soil at two soil moisture levels (wet and dry state) (Tekeste et al., 2007). Soil layer

Soil bulk density (Mg m3)

Dry state soil moisture content (%, d.b.)

Wet state soil moisture content (%, d.b.)

Top Middle Bottom

1.27 (0.09) 1.64 (0.11) 1.48 (0.06)

3.02 (0.32) 4.97 (0.21) 5.43 (0.26)

6.08 (1.95) 8.08 (2.22) 9.53 (3.25)

Literaturea

Tool (steel) properties Shear modulus (Pa) Poisson’s ratio Particle density (kg/m3)

8e + 10 0.3 7800

Soil-soil interaction Coefficient of restitution Coefficient of static friction Coefficient of rolling friction

0.01 0.49 0.5

Literatured Literaturea Literaturea

Soil-tool (steel) interaction Coefficient of restitution Coefficient of static friction Coefficient of rolling friction

0.01 0.37 0.05

Literatured Literaturea Literaturec

HM with bonding (initial parameters) Normal stiffness per unit area 1e + 08 (Kn) (N/m3)e Shear stiffness per unit area (Kt) 4e + 07 (N/m3) Bonded disk radius (rb) (mm) 1.5

Assumed Kt = 0.4 Kn Same as particle radius (r)

a

Tekeste et al. (2009). Calibrated implies the particle density was adjusted during particle generation for the top, middle and bottom layers to match laboratory measured soil bulk density according (Mousaviraad et al., 2017). c The material properties for tool (steel) were default values from EDEM 2018 material data base. d Tekeste et al. (2019). e Slope (modulus) estimate from soil cone penetration versus soil depth profile on Norfolk sandy loam soil. The normal stiffness per unit area was assumed 2.5 times the shear stiffness per unit area. b

Poisson’s ratio and soil shear modulus were obtained from finite element analysis of soil cone penetration testing on layered Norfolk sandy loam soil (Tekeste et al., 2009). Soil-soil and soil-tool (steel) interaction parameters of coefficient of restitution and coefficient of rolling friction, were assumed to have less influence on quasi-static cone penetration resistance and were fixed as the values reported in (Tekeste et al., 2019), from the simulation of soil reaction forces on cultivator sweeps in loam soil. For the parallel bond stiffness value in HM with parallel bond, initially the normal stiffness value of 1e + 08 N/m3 was used (Tekeste et al., 2019). Model parameters used in the initial soil cone penetration model are listed in Table 2. A series of cone penetration testing simulations was run at seven values of Kn (N/m3), the bond normal stiffness per unit area, (1.0e + 08, 2.05e + 08, 2.05e + 09, 3.75e + 09, 8.75e + 09, 1.34e + 09, 1.66e + 10, and 2.05e + 10). Maximum and minimum values of the bond stiffness were estimated from Sadek and Chen (2015), as they reported values of equivalent bond normal stiffness per unit area (Kn) of 1.1e + 10 N/m3 (dry-compacted) and 1.4e + 08 N/m3 (wet-loose) for sandy soil, and from the estimated slope of soil cone penetration resistance versus soil depth for Norfolk sandy loam soil (Tekeste et al., 2007).

47

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

(A)

(B) Top layer 36 mm

Middle layer 47 mm

Table 4 DEM parallel bond stiffness per unit area levels for sensitivity study. The parallel bond model was used only for the particle-to-particle interaction. For all the EDEM runs, the HM was used for particle-to-geometry contact physics. All the other DEM parameter values from Table 2 were used for the simulations. DEM input parameter

Value

Reference

Normal stiffness per unit area (Kn) (N/ m3) Shear stiffness per unit area (Kt) (N/m3)

2.05e + 08, 2.05e + 09, 3.75e + 09, 8.75e + 09, 1.34e + 09, 1.66e + 10, 2.05e + 10 1e + 05, 1e + 06, 1e + 07, 1e + 08

To be calibrated

Bottom layer 65 mm Table 5 DEM parameters calibrated to predict soil cone penetration testinga with values matching soil bin data (Reaves et al., 1969). The DEM parameters were used for simulation of the soil-blade interaction. Fig. 4. ASABE standard 30° cone penetrometer (A), DEM particle assembly packed in the top layer, middle layer, and bottom layer (B) similar to experimental test setup in Tekeste et al. (2007). The DEM particle assembly consists of a single sphere DEM shape that approximated the soil aggregate sample passing through a 0.85 mm screen. The cone with the shaft was inserted at 30 mm/s according to the ASABE Standard (ASABE 2018b).

The corresponding values of bond shear stiffness per unit area, Kt (N/m3), were generated assuming bond normal stiffness as 2.5 times that of the bond shear stiffness, to allow a compressive shear failure zone during conical tip indentation into the soil according to the cavity-expansion theory (Bishop et al., 1945). The parallel bond parameters of critical normal and shear stresses were assigned values of 1e + 12 Pa after Tekeste et al. (2019). The bond disk radius multiplier was assumed to be 1.5 mm, the same as the radius of the soil particle (see Table 3). Using the virtual soil cone penetration experiments, with the bond stiffness model parameter as an independent variable and predicted mean soil cone penetration resistance (0–120 mm soil depth) as the response variable, a best regression fitted metamodel was developed. An inverse optimization algorithm after Mousaviraad et al. (2017) with the objective of matching experiment data, was used to determine the optimal bond stiffness value. 2.4. DEM simulation formulation 2.4.1. Soil cone penetration testing for DEM calibration Cone penetrometer testing is a simple and easy method for characterizing in-situ soil strength. Modeling soil penetration resistance from a conical tip penetrometer for DEM material calibration can be helpful to predict dynamic soil behaviors including shear, compression, and cutting (Gill and Vanden Berg, 1968). The ASABE 30° conical tip and a cone base of 12.53 mm diameter (ASABE Standards, 2018a) penetrating a soil at 30 mm/sec was simulated to calibrate a DEM soil model for Norfolk sandy loam soil. Soil cone penetration data on Norfolk sandy loam soil reported in Tekeste et al. (2007) was used for the DEM material properties estimation and calibration to match the soil cone penetration resistance. The DEM soil model calibrated to the cone penetration testing on Norfolk sandy loam soil was assumed to simulate the soil reaction forces from the soil-bulldozer blade experiment data on the same soil as reported in Garrity et al. (1968) and Reaves et al. (1969). A virtual container (115 mm diameter and 148 mm long) with three stratified layers to reproduce the initial soil layers from the soil cone penetration experiment on Norfolk sandy loam soil (Tekeste et al., 2007) (Table 2) was created filling each layer using the random particle factory generator in EDEM 2018. The virtual

DEM input parameter

Value

Particle (single-sphere) Particle radius (r) (mm)

1.5

Soil material properties Poisson’s ratio Particle density (kg/m3)

Shear modulus (Pa)

0.32 2650 (Top layer) 2675 (Middle layer) 2675 (Bottom layer) 2.26e + 07 (Top layer) 2.26e + 07 (Middle layer) 2.26e + 07 (Bottom layer)

Tool (steel) properties Shear modulus (Pa) Poisson’s ratio Particle density (kg/m3)

8e + 10 0.3 7800

Soil-soil interaction Coefficient of restitution Coefficient of static friction Coefficient of rolling friction

0.01 0.565 0.5

Soil-tool (steel) interaction Coefficient of restitution Coefficient of static friction Coefficient of rolling friction

0.01 0.44 0.05

HM with bonding (initial parameters) Normal stiffness per unit area (Kn) (N/m3) Shear stiffness per unit area (Kt) (N/m3) Bonded disk radius (rb) (mm)

5.5e + 09 2.2e + 09 1.5

a The calibration approach deployed includes matching the soil-to-soil and soilto-tool DEM coefficients of restitution and friction (static and rolling) to soil-to-tool interaction modeling explained in Tekeste et al. (2009) and Tekeste et al. (2019), matching the initial bulk density of soil layers (Tekeste et al., 2007), and applying DOE experiment and optimization to obtain normal stiffness per unit area and shear stiffness per unit area values. The accuracy of the soil cone penetration resistance refers to comparing the DEM-predicted soil cone penetration profile and RMSE to laboratory data on Norfolk sandy loam soil.

container was filled with 60,739 particles of the single-sphere DEM particle diameter (Pd = 3 mm). The 3 mm diameter singlesphere DEM particle was approximated from the smallest 3D scanned soil aggregate (Fig. 3A). Considering all the 3D-scanned and DEM approximate soil shapes (Fig. 3), the 3 mm particle (Fig. 3A) produced the 4.2 ratio of the characteristic lengths of the ASABE cone (CPbase = cone base diameter = 12.5 mm) to the DEM particle diameter (Pd = 3 mm) and is assumed to ensure three particles randomly in contact with the cone at once. Other studies (Jiang et al., 2006; Falagush et al., 2015, Janda and Ooi, 2016; and Syed et al., 2017) showed a cone base diameter to the DEM particle size ratio (CPbase/Pd) within 1.8–3.0 gave desired soil penetration resistance. Iterative DEM particle generation and packing according to Tekeste et al. (2019) was performed to create an initially stable

48

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

Table 6 Tool dimensions and test parameters. The particle geometry mean (d

geo-mean)

was 19.9 mm estimated from the soil aggregates retained on the screen size of 12.5 mm.

Blade numbera

Length scale (n)b

Blade width (mm)

Blade height (mm)

Cutting depth (mm)

‘‘4” ‘‘5” ‘‘6” ‘‘7”

0.24 0.14 0.10 0.05

785 470 328 164

262 156 109 54.0

52.4 31.2 21.8 10.8

(3dgeo-mean) (1.6dgeo-mean) (1.1dgeo-mean) (0.5dgeo-mean)

a

Blade number in this study corresponds blade number 4, 5, 6, and 7 in Reaves et al. (1969). Blade scale length (n) is the ratio of width and height of the scaled blades to the market bulldozer blade dimensions of 3277 mm width and 1090 mm height (Reaves et al., 1969). b

DEM soil particle assembly and values of predicted soil bulk density within each layer which were close to the experimental results. A ratio of kinetic energy to the potential energy as low as 1e-08 was assumed to be stable, similar to results from Janda and Ooi (2016). The initial soil bulk density values in the DEM soil column chamber were 1.35 Mg/m3, 1.65 Mg/m3, and 1.51 Mg/m3 in the top (36 mm height), middle (hardpan) (47 mm height), and bottom (65 mm height) layers, respectively (Fig. 4). The initial soil bulk density of each layer in the experiment (Tekeste et al., 2007) was reproduced with a relative error of 6.3%, 0.3%, and 2.2% for the top, middle (hardpan), and bottom layers, respectively. Stainless steel (ASABE, 2018a) material properties were used for the cone penetrometer, the blade, and the soil box walls. 2.4.2. DEM soil-bulldozer blade simulation procedure For the soil-bulldozer blade interaction simulation, a DEM particle bed assembly was created in EDEM Academic (EDEM, 2011) using a DEM soil shape that approximated the 3D scanned and reconstructed soil aggregate sample retained on the 12.5 mm ASTM sieve size (Fig. 3D). The virtual soil bin dimensions were 1000 mm long (50dgeo-mean), 1000 mm (50dgeo-mean) wide and 254 mm (13dgeo-mean) depth. The total number of the particles was 43,451 clumped 3-sphere DEM particles. Four scaled blades were generated in CAD for the analysis with the tool parameters in Table 6 according to the specifications in Reaves et al. (1969). During the simulation, the tool cutting depth was kept constant at 20% of the blade height, a 45° blade lift angle was used, and the travel speed was 0.22 m/s, based on Garrity et al. (1968). The material properties calibrated using soil cone penetration testing were used for the soil-bulldozer blade simulation. 2.4.3. DEM fixed time step estimation For reasonable DEM calculation accuracy, EDEM (2011) recommends a time step 20% to 30% of the Rayleigh time step (TR) based on a wave speed across the smallest sphere (Rsmallest in Fig. 3D was 5 mm) during contact collision as proposed by Thornton and Randall (EDEM 2011). For Euler-based central time difference explicit DEM calculation in EDEM 2018, the DEM particle assembly exploded at the 20% TR (2.07e-05 s) making the simulation unstable. Fixing the time step (Dt) was instead approximated for wave speed propagation across the parallel bond beam which showed better stable simulation from bonded particle simulation (Guo et al., 2013). Eq. (4) was used to estimate a fixed time step using the bond shear modulus (G) instead of particle-to-particle contact shear modulus. Shear modulus (G) was estimated from G = E/(2* (1 + t)). Young’s modulus (E) of the bond was estimated as E = Kn*lbond/Abond, where Kn is bond stiffness, lbond (bond length) = 2rb and Abond (bond cross sectional area) = p r2b, and rb is disk radius (EDEM, 2011). Poisson’s ratio (t) was assumed 0.3 and particle density (q) was assumed as 2650 kg/m3 for time step estimation. We fixed the time step (Dt) at 1e-06 s (approximately 1.78% of TR) for soil-cone penetrometer simulation (soil Rsmallest was 1.5 mm) and 3.19e-06 s (3.09% of TR) for the soil-bulldozer blade simulation (soil Rsmallest

was 5 mm). Both soil-cone penetrometer and soil-blade simulation with their respective fixed time steps gave a stable DEM simulation, and the ratio of kinetic energy to potential energy was less than 1e-08 according to Janda and Ooi (2016).

qffiffiffi

pR qG Dt <¼ T R ¼ 0:1631m þ 0:8766

ð4Þ

where, TR is Rayleigh time step; R = Radius of the smallest particle (Rsmallest); q = Particle density; G = Shear modulus, and m = Poisson’s ratio. All the DEM simulations were run in EDEM 2018 using 12-CPU cores on a Dell Precision T7910 (2.30 GHz processor speed and 64 GB RAM). Running the simulations at their fixed time steps, it took 23 h/s and 22 h/s normalized elapsed computation time (CPU) to complete a 3 s blade simulation time and 5 s cone penetration testing, respectively. For the soil cone penetration simulation, the vertical reaction force on the cone tip was sampled every 0.01 s. Soil cone penetration resistance (kPa) was calculated by dividing the vertical reaction force on the cone by the 123 mm2 cone base area (ASABE Standards, 2018a). Similarly from the soil-blade simulation, the soil horizontal and vertical reaction forces on the blade were sampled every 0.01 s. From the blade force versus time series data, steady-state forces that estimate the full soil load forces were extracted from every blade simulation. The total number of soil DEM particles on top of the original soil surface was obtained to estimate the full soil load by multiplying the total number of soil particles by the single DEM soil particle mass of 6.4311e-3 kg. 2.5. Data analysis In the JMP Pro 14.1 software (JMPÒPro14.1. SAS Institute Inc., Cary, NC, 1989–2007), a regression model using the GLM procedure was created to determine the meta-model relationship between the DEM bond normal stiffness and mean soil cone penetration resistance (0–120 mm depth). A power fit equation of a form similar to Eq. (3) was fit to the data of the DEM-predicted soil horizontal forces on the bulldozer blade. We used the length scale (n) and a nonlinear curve fitting technique, the LevenbergMarquardt optimization algorithm, which was a default backward error propagation algorithm in MATLABÒ software (Math Works, Inc., Natick, Mass.). The nonlinear least squares fitting was the preferred option because it gave a minimum Root Mean Square Error (RMSE) and best coefficient of determination (R2). Comparison between the DEM-predicted soil forces and the soil bin experimental data was performed using a linear regression curve fitting in MATLABÒ software.

49

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

3. Results and discussion 3.1. DEM model parameters sensitivity analysis The relationship between soil cone index (0–120 mm depth) from the DEM soil cone penetration resistance and the DEM parallel bond normal stiffness was best explained using CIDEM = 232 + 9. 14e-08Kn,DEM (R2 = 0.975859). The best regression meta-model with the initial parameter values as in Table 4 was selected for the optimization procedure. Using the best least square regression model, the ‘‘Desirability” optimization procedure in the JMP Pro 14.1 statistical analysis package was used to determine the optimized DEM bond stiffness. During the optimization step, the DEM-predicted response of soil cone penetration was targeted to the laboratory soil cone index (734 kPa), a mean value from the Norfolk sandy loam soil bin test (Reaves et al., 1969). A bond stiffness (Kn) of 5.5e-09 N/m3 was found to be the optimal value with desirability of 0.98. All simulations of the scaled bulldozer blades were run with DEM parameter values (Table 5), that gave the best estimate of the soil responses from soil cone penetration testing and the mean soil cone index value for the soil bin data reported in Reaves et al. (1969) for the Norfolk sandy loam. The initial estimates of the coefficients of static friction for the soil-soil and soil-

Fig. 6B. Norfolk sandy loam soil with middle layer (bulk density of 1.64 Mg/m3) in a relatively dry soil condition (soil moisture content of 4.17%, d.b.). The RMSE between the test and simulation for 0–12 cm depth was 850 kPa.

Blade- “4” (Length scale, n = 0.24) (L/W =0.57)

Blade- “5” (Length scale, n = 0.14) (L/W = 0.83) Fig. 5. Soil failure patterns as the cone tip (with shaft) advances into the soil layers from DEM simulation of soil cone penetration testing with the initial soil layering as in Fig. 4. Soil failure modes (A) Cone within the top soil layer with spiral-soil failure mode, (B) Cone within the middle layer with cavity-expansion soil failure mode, and (C) Cone within the bottom layer with cavity-expansion soil failure mode.

Blade-“6” (Length scale, n = 0.10) (L/W = 0.87)

Blade-“7” (Length scale, n = 0.05) (L/W = 1.8) Fig. 7. Soil cutting in front of scaled blades 4, 5, 6, and 7 at the steady-state conditions.

Fig. 6A. Norfolk sandy loam soil with middle layer (bulk density of 1.64 Mg/m3) in a relatively wet soil condition (soil moisture content of 8.78%, d.b.). The RMSE between the test and simulation for 0–12 cm depth was 48 kPa.

tool (steel) from Tekeste et al. (2007) were also adjusted to improve the prediction of the soil cone penetration data. The soil deformation in front of the advancing cone in the top layer (Fig. 5) appeared to exhibit behavior similar to the bearing capacity failure theory (Yu and Mitchell, 1998). As the cone fully penetrated into the denser middle and bottom layers, the soil failure resembles the cavity-expansion theory from indention (Bishop et al., 1945). Soil cone penetration for the relatively wet soil moisture condition (Fig. 6A) and dry soil moisture condition (Fig. 6B) with the soil middle (hardpan) bulk density of 1.64 Mg/m3 showed good agree-

50

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

Fig. 8. Soil horizontal forces from DEM simulation of blade-4, blade-5, blade-6, and blade-7 running at 0.22 m/s and 20% percent of the blade height.

Fig. 9. Comparison of DEM-simulated and similitude model prediction from Norfolk sandy loam soil bin data (Reaves et al., 1969).

Fig. 10. Soil horizontal forces from the DEM simulation of soil-bulldozer blades and best power fit Eq. (5) predicting the soil horizontal forces as a function of bulldozer blade length scale (n).

ment between DEM and experimental data for Norfolk sandy loam soil. For the dry soil condition, DEM was not able to capture the soil cone penetration resistance profile for depths greater than 8 cm. The reason why the DEM underpredicted the cone penetration resistance for depths greater than 8 cm was not obvious. FEM analysis of cone penetration (Tekeste et al., 2009) using the DruckerPrager soil constitutive model also was not able to predict the soil cone penetration resistance for relatively dry soil at similar deep depths in the soil profile. Numerical representation of the natural soil drying process from the experiment may require further investigation.

scale and compared with the Norfolk sandy loam soil bin data (Garrity et al., 1968). The soil bin data were estimated from a curve fit using Eq. (3). As shown in Fig. 7, soil builds in front of the bulldozer blades during the cutting process. The load curve showed the full blade width of cut at different L/W ratios (Fig. 7). The load-growth curves from the soil horizontal forces appeared to show a similar trend as the data observed by Reaves et al. (1969). The oscillation of the soil horizontal force data showed the soil filling and spillage from the blades during the cutting process (Fig. 8). The mean soil horizontal reaction forces from the DEM simulation and the dimensionless prediction according to Eq. (3) after Garrity et al. (1968) were transformed to natural logarithmic soil horizontal data and plotted as a function of the length scale (n). The DEM-predicted soil horizontal forces as functions of blade length scaling (from 0.05 to 0.24) were in good agreement with the dimensional analysis-predicted model from Norfolk sandy

3.2. DEM soil-bulldozer blade simulation and verification with experimental tests At the steady-state cutting of soil (Fig. 7), the soil horizontal forces from the DEM simulation are plotted as a function of length

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

Fig. 11. The soil weight pushed by fully loaded blades, from DEM simulation and power fitted models to the data from the DEM simulation of soil-bulldozer blades.

loam soil bin data (linear regression with slope = 0.9634 and R2 = 0.9965, Fig. 9). At the length scale (n) of 0.05 (blade7), the DEM overpredicted the steady-state soil horizontal forces which might be attributed to having larger DEM particles relative to the blade parameters (164 mm blade width and 54 mm blade height). In the DEM simulation, the blade reached full load at a shorter blade travel length and large DEM particles were shown to overtop the blade. Reaves et al. (1969) also noted soil loading behavior from blade-7 consistently different from the other blades. The soil conditions from this blade test might have been more heterogeneous than for the other blade tests. The accuracy of soil horizontal forces on the 164 mm blade (n = 0.05) predicted by DEM needs further improvement; for example, running the DEM simulation using a smaller particle size which may require greater CPU time per simulation time. Blade-7 is a small blade, so a 3D-printed blade can also be run to obtain soil bin data. The relationship between the DEM-predicted soil horizontal force (DDEM) and length scale (n) was estimated using a power fit equation Eq. (5). The nonlinear best fit to Eq. (5) resulted in a RMSE of 94.14 N (R2 = 0.9976) with model coefficient estimates of ADEM = 93,670 N and BDEM = 2.284.

DDEM ¼ ADEM ðnÞBDEM

ð5Þ

Using Fig. 10, DEM simulation of a model bulldozer can be used to predict soil horizontal force performance of prototype bulldozer blades on Norfolk sandy loam soil. After the blades were fully loaded, the soil full loads were also estimated by multiplying the individual DEM soil particle mass by the total number of DEM soil particles above the original soil surface in the DEM soil box. Reaves et al. (1969) assumed fully loaded blades when the soil spillage around the width of the blade was equal to the new soil supplied to the blade pile from the constant depth of cut. Even though no report was available about soil weight on the fully loaded blades from the soil bin experiment in Garrity et al. (1968) and Reaves et al. (1969), a power fit similar to the soil horizontal force (Eq. (5)) estimated the DEM soil weight pushed by the fully loaded blades and the length scale (Fig. 11) with the model coefficients of ADEM soil = 17,460 N and BDEM soil = 2.822, with RMSE = 17.15 N and R2 = 0.9903. 4. Conclusion Simulation of soil-bulldozer blade interaction was successfully performed using the DEM technique. Approximating soil DEM shapes for this soil-tool interaction problem was done using 3D-

51

scanned and reconstructed soil aggregates. Using a soil model represented as 3 mm diameter single spheres (Cone base diameter (CPbase)/Particle diameter (Pd) = 4.2) and DEM material properties calibrated to the mean soil cone penetration resistance from experimental data, DEM predicted the soil cone penetration profile in good agreement with the experimental data. Using the calibrated soil DEM material properties and a larger particle size (d geo-mean of 19.9 mm), a power fit scaling law was established between soil reaction forces and the length (n) for four scaled bulldozer blades. Comparison of the DEM power fit showed good agreement with dimensional analysis-predicted performance of the blades on Norfolk sandy loam soil bin data from the literature (Garrity et al., 1968, and Reaves et al., 1969). Performance of machine-scale bulldozer blades interacting with cohesionless soil (e.g. Norfolk sandy loam) can be virtually evaluated using the power fit scaling from DEM simulation. Future research to validate the simulationbased dimensional analysis of soil-bulldozer blade on cohesive soils can be done using a similar DEM technique. Declaration of Competing Interest We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. References ASABE Standards, 2018a. S313.3 FEB1999 (R2018). Soil cone penetrometer. St. Joseph, Mich.: ASABE. ASABE Standards, 2018b. EP542 FEB1999 (R2018). Procedures for using and reporting data obtained with the soil cone penetrometer. St. Joseph, Mich.: ASABE. Asaf, Z., Rubinstein, D., Shmulevich, I., 2007. Determination of discrete element model parameters required for soil tillage. Soil Tillage Res. 92 (1–2), 227–242. Bishop, R.F., Hill, R., Mott, N.F., 1945. The theory of indentation and hardness. Proc. Phys Soc. 57 (3), 147–159. Chen, Y., Munkholm, L.J., Nyord, T., 2013. A discrete element model for soil-sweep interaction in three different soils. Soil Tillage Res. 126, 34–41. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. EDEM, 2011. EDEM theory reference guide. Edinburgh, UK: DEM Solutions. Falagush, O., McDowell, G.R., Yu, H.S., 2015. Discrete element modeling of cone penetration tests incorporating particle shape and crushing. Int. J. Geomech. 15 (6), 04015003. Garrity, R.J., Kolthoff, C.P., Reaves, C.A., Schafer, R.L., 1968. A test comparison of model and full-size bulldozer blades. SAE Trans. 77, 2495–2511. Gill, W.R., Vanden Berg, G.E., 1968. Soil dynamics in tillage and traction. Handbook 316, Agr. Res. Service, U.S. Dept. Agriculture, Washington, D.C. Guo, T., Wassgren, C., Hancock, B., Ketterhagen, W., Curtis, J., 2013. Validation and time step determination of discrete element modeling of flexible fibers. Powder Technol. 249, 386–395. Hillel, D., 1998. Environmental Soil Physics. Academic Press, San Diego, CA. Hofstetter, K.W., 2002. Analytic method to predict the dynamic interaction of a dozer blade with earthen material. In : Proc. 14th International Conference of the International Society for Terrain-Vehicle Systems, Vicksburg, MS. Janda, A., Ooi, J.Y., 2016. DEM modeling of cone penetration and unconfined compression in cohesive solids. Powder Technol. 293, 60–68. Jiang, M.J., Yu, H.-S., Harris, D., 2006. Discrete element modelling of deep penetration in granular soils. Int. J. Numer. Anal. Meth. Geomech. 30 (4), 335– 361. Mousaviraad, M., Tekeste, M., Rosentrater, K.A., 2017. Calibration and validation of a discrete element model of corn using grain flow simulation in a commercial screw grain auger. Trans. ASABE 60 (4), 1403–1415. Obermayr, M., Vrettos, C., Eberhard, P., Däuwel, T., 2014. A discrete element model and its experimental validation for the prediction of draft forces in cohesive soil. J. Terramech. 53, 93–104. Potyondy, D.O., Cundall, P.A., 2004. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 41 (8), 1329–1364. Reaves, C.A., Schafer, R.L., Garrity, R.J., Kolthoff, C.P., 1969. Similitude of bulldozer blades. Trans. ASAE 12 (5), 577–579. 583. Sadek, M.A., Chen, Y., Jude, L., 2011. Simulating shear behavior of a sandy soil under different soil conditions. J. Terramech. 48, 451–458. Sadek, M.A., Chen, Y., 2015. Feasibility of using PFC3D to simulate soil flow resulting from simple soil-engaging tool. Trans. ASABE 58 (4), 987–996. Schafer, R.L., Reaves, C.A., Young, D.F., 1969. An interpretation of distortion in the similitude of certain soil-machine systems. Trans. ASABE 12 (1), 145–149.

52

M.Z. Tekeste et al. / Journal of Terramechanics 88 (2020) 41–52

Sakaguchi, E., Ozaki, E., Igarashi, T., 1993. Plugging of the flow of granular materials during the discharge from a silo. Int. J. Mod. Phys. B 7, 1949–1963. Shmulevich, I., Asaf, Z., Rubinstein, D., 2007. Interaction between soil and a wide cutting blade using the discrete element method. Soil Tillage Res. 97, 37–50. Syed, Z., Tekeste, M., Way, T.R., 2017. Discrete element modeling (DEM) of cone penetration testing on soil with varying relative soil density. ASABE Paper 1701608. St. Joseph, Mich.: ASABE. Tekeste, M.Z., Raper, R.L., Tollner, E.W., Way, T.R., 2007. Finite element analysis of cone penetration in soil for prediction of hardpan location. Trans. ASABE 50 (1), 23–31. Tekeste, M.Z., Tollner, E.W., Raper, R.L., Way, T.R., Johnson, C.E., 2009. Non-linear finite element analysis of cone penetration in layered sandy loam soilconsidering precompression stress state. J. Terramech. 46 (5), 229–239.

Tekeste, M.Z., Balvanz, L.R., Hatfield, J.L., Ghorbani, S., 2019. Discrete element modeling of cultivator sweep-to-soil interaction: Worn and hardened edges effects on soil-tool forces and soil flow. J. Terramech. 82, 1–11. Ucgul, M., Fielke, J.M., Saunders, C., 2015. Defining the effect of sweep tillage tool cutting edge geometry on tillage forces using 3D discrete element modelling. Inf. Process. Agric. 2 (2), 130–141. Walton, O.R., Braun, R.L., 1986. Stress calculations for assemblies of inelastic spheres in uniform shear. Acta Mech. 63, 73–86. Yu, H.S., Mitchell, J.K., 1998. Analysis of cone resistance: Review of methods. J. Geotech. Geoenviron. Eng. 124 (2), 140–149.