Calibration of granular material parameters for DEM modelling and numerical verification by blade–granular material interaction

Available online at www.sciencedirect.com

Journal of Terramechanics Journal of Terramechanics 46 (2009) 15–26 www.elsevier.com/locate/jterra

Calibration of granular material parameters for DEM modelling and numerical verification by blade–granular material interaction C.J. Coetzee *, D.N.J. Els Department of Mechanical Engineering, University of Stellenbosch, Private Bag X1, Matieland 7602, South Africa Received 15 February 2007; received in revised form 25 August 2008; accepted 14 December 2008 Available online 3 February 2009

Abstract The discrete element method (DEM) is a promising approach to model blade–granular material interactions. The accuracy of DEM models depends on the model parameters. In this study, a calibration process was developed to determine the parameter values. The particle size was the same as the real material and the particle shape was modelled using two spherical particles rigidly clumped together to form a single grain. Laboratory shear tests and compressions tests were used to determine the material internal friction angle and stiffness, respectively. These tests were replicated numerically using DEM models with different sets of particle friction coefficients and particle stiffness values. The shear test results are found to be dependent on both the particle friction coefficient and the particle stiffness. The compression test results show that it is only dependent on the particle stiffness. The combination of shear test and compression test results can be used to determine a unique set of particle friction and particle stiffness values. The calibration process was validated experimentally and numerically by modelling a blade moving through granular material. Results show that the forces acting on the blade can be accurately modelled with DEM and the maximum error is found to be 26%. The relative particle-blade displacements were used to predict the position and shape of the shear lines in front of the blade. A good qualitative correlation was achieved between the experiments and the DEM simulations. Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved.

1. Introduction Earthmoving equipment plays an important role in the agricultural, earthmoving and mining industries. The equipment is highly diverse in shape and function, but most of the equipment can be categorised into one of three principal classes, namely blades, rippers and buckets (shovels). This paper focuses on the interaction between granular material and blades. Bulldozers make use of blades for material cutting in a process which involves many material and blade properties [1]. The material strength, the build-up heap in front of the blade and the blade geometry are all important parameters. In order to improve the design of blades, reliable models are needed which take all the above parameters into consideration. *

Corresponding author. Tel.: +27 21 808 4239; fax: +27 21 808 4958. E-mail address: [email protected] (C.J. Coetzee).

The design of effective and efficient blades begins with the analysis forces and energy required by the machinery. Karmakar and Kushwaha [2] summarised the different approaches used to model the cutting process. Five major methods can be identified, namely empirical and semiempirical analytical methods, finite element methods (FEMs), discrete element methods (DEMs) and computational fluid dynamics (CFD). Semi-empirical analytical models ([3–7]) are based on Terzaghi’s passive earth pressure theory and the tool configuration [8]. The methods of Hansen [9] and Perumpral et al. [10] assume that the material in front of the blade fails on a pre-determined failure plane, and then moves as a rigid body. The sum of all the forces acting on this rigid body can be used to solve for the forces acting on the blade. Sokolovski [11] used a different method based on the theory of plasticity. This theory is based on the assumption that a state of failure exists at any point within a certain area (rupture zone) or on a certain curve (rupture line).

0022-4898/$36.00 Ó 2009 ISTVS. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jterra.2008.12.004

16

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

A set of differential equations, derived from the MohrCoulomb yield criterion and the equations of equilibrium, is numerically solved using the method of characteristics and appropriate boundary conditions and material properties. The analytical methods are based on passive earth pressure theory and assumptions of a preliminary material failure pattern [2]. Dynamic effects, complicated tool geometry (such as buckets) and large deformations are difficult to consider in analytical modelling [2]. The finite element method can overcome some of these difficulties. The finite element method supplies more information regarding the progressive material failure zone, material stresses, displacement, velocity and acceleration [2]. Difficulties encountered by FEM include large deformation and complex flow patterns which can lead to severe mesh distortion [12]. Continuum based methods, such as FEM, cannot easily model granular phenomena such as material cracking, mixing of different material layers and segregation. These disadvantages of FEM can be solved using DEM. It has been shown that CFD can be used to model soil– tool interaction [2]. The preliminary results presented by Karmakar and Kushwaha are very encouraging, but further research is needed before CFD can be used to model soil–tool interaction with confidence. The discrete element method is a promising approach to model soil-implement interaction [13] and can be used to overcome some of the difficulties encountered by analytical methods and FEM. In DEM, the failure patterns and shape of the up heave are not needed in advance. The tools are modelled using a number of flat walls and the complexity of the tool geometry does not complicate the DEM model. Large deformation and the development of the material free surface (build-up heap) are automatically handled by the method. The main problem with DEM is to determine the material micro-properties (spring stiffness, friction coefficient, damping values, bonding strength). The micro-properties should be specified such that the flow of thousands of particles behaves in the same way as the real granular material (macro- oedometer properties). Laboratory experiments (e.g. shear tests, biaxial/triaxial tests and tests) or in-situ tests are necessary to determine these properties before any useful modelling and predictions can be made [14]. Very few researchers have tried to determine the microproperties experimentally. Vu-Quoc et al. [15] measured the coefficient of restitution in soybeans by performing simple drop tests from various heights. Data from high-speed video recordings was used to determine the particle properties. Tanaka et al. [16] conducted bar penetration tests and compared the results with those obtained from DEM simulations. They measured the material density, but the contact stiffness was chosen without any experimental validation. By comparing the movement of the particles during the experiment with the movement of the particles during the simulation, Tanaka et al. [16] could determine the friction coefficient that gave the best results.

Recently, Franco et al. [1] proposed an inverse calibration method to determine the micro-properties for modelling soil–bulldozer blade interaction. In their approach, the material was assumed cohesionless. Based on DEM results, the particle friction coefficient and stiffness were determined from energy principles and direct shear tests. The draft force and vertical force acting on the blade, as predicted by DEM, were then compared with McKeys’s calculation model [17]. The calibration process was based on direct shear test DEM simulations and blade–soil DEM simulations but experimental validation was not performed. Asaf et al. proposed an in-situ method for determining the micro-properties [14,18]. Their method was based on wedge penetration tests and a non-linear optimisation scheme. The methodology was validated by using DEM simulation and in-situ tests. Mishra and Murty [35] determined the contact stiffness and damping of a steel ball by performing drop tests and measuring the contact force using an ultra fast load cell. Landry et al. [36] determined the properties of manure products using direct shear tests. A sensitivity study was used to determine the effect of the various DEM parameters on the shear test results. The first objective of this study was to develop a calibration process in which the DEM material parameters could be determined. For this purpose, experimental and DEM results from shear tests and compression tests were used. The second objective was to validate the calibration method and demonstrate to what degree DEM could accurately predict blade–granular material interaction. For this purpose, results from a DEM simulation of a blade cutting through corn grains were compared to experimental results. The study is limited to two-dimensional models and cohesionless granular material. 2. The DEM model Discrete element methods are based on the simulation of the motion of granular material as separate particles. DEM was first applied to soils by Cundall and Strack [19]. Calculations performed during a DEM simulation cycle alternate between the application of Newton’s second law (applied to the particles) and a force–displacement law (applied at the contacts between particles). In this study, all the simulations were two-dimensional and performed using the commercial DEM software PFC2D [20]. Using the soft particle approach, each contact is modelled with a linear spring in the contact normal (stiffness kn) and a linear spring in the contact tangential (stiffness ks) direction as depicted in Fig. 1. Frictional slip is allowed in the tangential direction with a friction coefficient l. The particles are allowed to overlap and the amount of overlap is used in combination with the spring stiffness to calculate the contact force components. The contact force in the normal direction is given by X k n DU n ð1Þ Fn ¼ 

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

17

In this research, corn grains were used. Since only cohesionless material was under consideration, the material parameters investigated were the particles size and shape, particle density, particle friction coefficient and particle contact stiffness. The default damping C = 0.7 proposed by PFC2D [20] was used in all simulations and the width (dimension into the plane) was taken as 0.2 m (the width of the blade test rig). 3.1. Particle size and shape

kn ks

Fig. 1. DEM contact model.

where DU n is the increment in overlap in the contact normal direction. The contact shear force is given by P k s DU s for jF s j < jlF n j P ð2Þ Fs ¼ lF n signð DU s Þ for jF s j P jlF n j where DU s is the displacement increment in the contact tangential direction. The damping force acts on a particle in the opposite direction to the global particle velocity and is proportional to the resultant force applied on the particle with a damping coefficient C. For a detailed description of DEM, the reader is referred to [19,21–23]. 3. Calibration of material parameters

R 2.5 - 4.5 5-6 8 - 12

Rowlands [24] observed that seed grains are suitable for experimental testing and closely resemble natural soil flow into dragline buckets. The grains have a relatively low friction coefficient with glass, which makes it a good material to be used in the experimental setup (see Section 4). The seed grains were also found suitable for DEM simulations. The size of the particles was large enough to be modelled accurately. If, for example, sand were used, it would not be possible to model each individual grain, and the grain size to be used in the DEM model would become another unknown parameter.

5-9

3.0 - 5.0

Friction μ

The particles were chosen so that it roughly represented the physical shape and size of the corn grains. Particles were constructed by forming clumps. Clumps can be formed by adding two or more particles (discs in 2D and spheres in 3D) together to form one rigid particle, i.e. particles included in the clump remain at a fixed distance from each other [20]. Particles within a clump can overlap to any extent and contact forces are not generated between these particles. Clumps cannot break up during simulations regardless of the forces acting upon them. Asaf et al. [13] investigated the effect of particle structure on internal friction angle. They used two equally sized discs, clumped together (Fig. 2b), and biaxial tests to investigated the effect of particle centre distance on the internal friction angle. Results showed that internal friction increases with an increase in centre distance with an asymptotic behaviour. It was noticeable that beyond a certain centre distance, the value of internal friction remains relatively constant although the centre distance is increased. Ni et al. [34] observed higher peak and internal friction angles for clumped particles when compared to single spheres. This is a result of higher ‘‘interlocking” between particles. Fig. 2 shows the ‘‘tear-drop” shape of the grains and the equivalent DEM clump used in this study. The diameter of the disc shaped particles used in the simulations was chosen from a uniform distribution within the range shown in the figure. Although the grain shapes could be modelled more accurately by using more than two particles in a clump, the representation chosen was thought to be accurate enough and computationally efficient.

R 1.5 - 3.0

4-5 3-6

Fig. 2. (a) Physical grain dimensions and (b) DEM grain model. Dimensions in [mm].

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

3.2. Particle density The material bulk density was measured by filling a cylindrical container with material. In the filling process, the material was allowed to flow from a hopper into the cylindrical container. Excess material above the upper edge of the container was carefully removed. The container with material was weighed and the internal dimensions of the container measured to calculate the volume. The bulk density was defined as the total material mass divided by the volume filled by the material. Three experiments were conducted to ensure repeatability. In the model, the experimental filling process was replicated using a hopper with the same dimensions and height above the container. All the particles were given an initial density (solid density) slightly higher (15%) than the measured bulk density. The assembly of particles was then allowed to reach static equilibrium and the bulk density calculated using the total mass of all the particles and the total volume filled by the assembly of particles. Using an iterative procedure, the density of all the particles was changed until the bulk density was within 0.1% of the measured value. This accuracy could be obtained within three to four iterations. After each iteration, the assembly of particles was allowed to reach a new static equilibrium state. 3.3. Direct shear tests A shearbox was used to experimentally and numerically (DEM) determine the corn internal friction angle and the friction angle between the corn grains and steel (blade). The box had dimensions of 375 by 375 mm with a total fill height of approximately 175 mm and a shear rate of 10 mm/min. Normal stresses were applied in the range 0.5–3 kPa. Three experiments were performed to ensure repeatability. Although the ratio of particle normal stiffness to particle tangential stiffness, kn/ks, has an influence on the materials Young’s modulus and Poisson’s ratio [20], this ratio was taken equal to one, i.e. the particle normal and tangential stiffness were assumed to be the same. Poisson’s ratio given by a 2D biaxial simulation can strictly not be compared with laboratory-measured values since conditions in a 2D biaxial test are neither plane strain nor plane stress according to PFC2D [20]. Consequently, Poisson’s ratio was not used in the calibration process. Asaf et al. [14] have shown that the internal friction angle is dependent on the normal stress applied to the shear test. The internal friction angle decreases with an increase in normal stress and it is important to apply an expected range of normal stress during testing. In this study, the normal stresses imposed during the shear tests were chosen to be within the range experienced during blade cutting experiments. A numerically equivalent DEM model with the same dimensions as the physical apparatus was constructed and tests were performed under the same conditions. The

shearbox was filled with approximately 3000 clumped particles. A constant normal stress was applied to the loading wall using a servo control algorithm available in PFC2D. The bottom half of the box was displaced at the shear rate and the horizontal force needed to move the bottom half of the box recorded. From the shear force, the shear stress could easily be calculated. The average shear stress between 2.1% and 4.8% strain (8 and 18 mm of displacement) was used to calculate the friction angle, which was defined as the average shear stress divided by the average normal stress, Fig. 3. A series of shear tests were performed using a range of particle stiffness and friction coefficient values. The shear simulations showed that the internal friction angle depended on both the particle stiffness and friction coefficient, especially in the low stiffness range. The shear test results are shown in Fig. 4. For stiffness values above 100 kN/m, all particle friction values l resulted in a higher internal friction value /. For example, with the particle friction l = 0.3, all the calculated internal friction angles were greater than tan1(0.3) = 16.7°. This is a direct result of the interlocking phenomenon and also observed by Asaf et al. The interlocking friction angle can be found by calculating the internal friction angle based on a shear test simulation where the particle friction coefficient l is set equal to zero. The interlocking angle is dependent on the particle size, size distribution, porosity and particle stiffness according to Asaf et al. [14]. From Fig. 4, it is clear that both the particle stiffness and friction coefficient had an influence on the material internal friction angle, /. For stiffness values below 100 kN/m the internal friction angle is strongly dependent on the particle stiffness, but for higher stiffness values the dependence decreases, and more so for lower friction coefficients. The reason for this is that with lower particle stiffness values, the overlap between contacting particles is greater and this decreases the interlocking effect. Particles can more easily slip over the other and this lowers the internal friction angle. Thus, the direct shear test cannot be used to deter-

1.0 0.9

DEM

0.8 Shear stress [kPa]

18

Experiment

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

2

4

6

8 12 10 Displacement [mm]

14

16

18

Fig. 3. Typical shear test result for a normal stress of 2 kPa.

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26 40 μ = 0.30

Internal friction angle φ [deg]

35

slight increase in / with an increase in l. They also showed that as the particle size distribution was made wider, the magnitude of / increased.

30

μ = 0.20

3.4. Confined compression test

25

μ = 0.15

In the confined compression test (also called the oedometer test), stress is applied to the specimen along the vertical axis, while strain in the horizontal directions is prevented. Shear stresses and shear strains as well as compressive stresses and volume changes occur in this test, but since the material is prevented from failing in shear, compression is the dominant source of strain. The experiments were performed with a 250 mm diameter cylinder and a specimen height of 80 mm. The cylinder was carefully filled and a load applied to the loading plate. The load and displacement was recorded. All experiments were performed at low compression rates (±2 mm/min) to avoid any dynamic effects, i.e. quasi-static conditions were assumed. The experiment was replicated numerically (Fig. 6) and since side friction will always play a role in this test it had to be included in the model. The wall friction was set to the specific material-steel friction coefficient. Although the experiment was performed in a cylindrical container and a 2D model was used, it is assumed that the strain is one-dimensional and a comparison can thus be made. Fig. 5 shows typical experimental results obtained from unloading and reloading (hysteresis). Roughly 23 of the strain that occurred during initial loading was recovered during subsequent unloading. The strains that resulted from sliding between particles were largely irreversible. The rebound upon unloading was caused by the elastic energy stored within the particles, as the assembly was loaded. However, there could be some reverse sliding between particles during unloading. Only a small amount of permanent strain resulted from the next two cycles. The numerical results (Fig. 6) showed little to no hysteresis. This might be due to the initial porosity of

μ = 0.10 20

15 Measured internal friction angle 10 0

100

200 300 400 500 Particle stiffness kn = ks [kN/m]

600

Fig. 4. The effect of particle stiffness and friction coefficient on the material internal friction angle.

mine a unique set of parameters. It is clear from Fig. 4 that more than one combination of the particle friction coefficient and stiffness values can result in an internal friction angle similar to the measured value (indicated by the dotted line). An additional test, the confined compression test, is used in the next section to determine the unique set of parameters. The friction angle with steel was measured by replacing the lower part of the shearbox with a steel sheet. The shear results correspond to the results by Asaf et al. [14] who have shown that the internal friction angle increases with an increase particle stiffness. Corriveau et al. [26] also studied the effect of particle friction on the internal friction angle numerically by using a linearspring-dashpot model and the Walton and Braun contact model [15] which allows for plastic deformation. They performed biaxial tests with an assembly of 1250 and 300 circular discs. Their results showed that the internal friction angle depended strongly on the particle friction coefficient in the lower range of l. In the upper range there was only a -90

-90

-80

-80

-70

-70

-60

-60

Axial stress [kPa]

Axial stress [kPa]

19

-50 -40 -30

Load

-50 -40 -30

-20

-20 -10

-10 0

-0.5

-1

-1.5

-2

-2.5

-3

-3.5

-4

Axial strain [%]

Fig. 5. Experimental load–displacement curve from confined compression test.

0

-0.5

-1

-1.5 -2 Axial strain [%]

-2.5

-3

Fig. 6. DEM load–displacement curve from confined compression test.

20

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

the material. During the experiments and particularly the first loading cycle, relative particle movement might occur, decreasing the porosity and compressing the material. During subsequent cycles, the particles have settled and relative particle movement is limited. The porosity was not measured and should be added to future calibration methods. Experiments and simulations were stopped after three cycles. The second and third cycles followed a stable hysteresis loop. These stable cycles were used to determine the stiffness of the system of particles. During a uni-axial tensile test, stress is applied in the axial z-direction and Young’s modulus of elasticity is defined as rz ð3Þ E¼ ez During confined compression, the direct strain in the xand y-directions is constrained to be zero and it can be shown that the relation between the confined Young’s modulus E0 and the usual Young’s modulus E, is given by   ð1 þ mÞð1  2mÞ ð4Þ E ¼ E0 1m where v is Poisson’s ratio and E is given by Eq. (3). The material’s confined Young’s modulus was calculated from the last loading cycle using r80  r10 E0 ¼ ð5Þ e80  e10 where r80 is an axial compression stress of 80 kPa and e80 the corresponding strain and r10 an axial compression stress of 10 kPa and e10 the corresponding strain. Although this approach ignores the non-linear load–displacement effect, the same procedure is used for both the experimental and numerical results. A series of simulations, each with different particle stiffness and friction values, were performed. Stiffness values of 10, 50, 100, 300, 450 and 500 kN/m were used with friction coefficients of 0.10, 0.12, 0.20 and 0.30 for each stiffness value, respectively. The results are summarised in Fig. 7.

It was found that the friction coefficient had little to no effect on the system stiffness and a linear relation existed between particle stiffness and the confined Young’s modulus. Deformation is predominantly due to distortion of the individual particles (particle stiffness) and not due to relative particle movement (friction). With both the shear test and compression test results available, a unique set of micro-parameters could be determined. The compression test showed to be independent of the particle friction within the range tested. A particle stiffness of 450 kN/m resulted in a confined Young’s modulus close to the measured value (Fig. 7). With the stiffness known, the shear results (Fig. 4) could be used to determine the particle friction coefficient. A value of l = 0.12, resulted in an internal friction angle / = 24° which is close to the measured value of / = 23°. 3.5. Angle of repose With the particle size, shape, stiffness and friction values known, the angle of repose /r was determined experimentally and numerically as a final test. One thousand grains were allowed to flow from a hopper at a fixed height of 100 mm and allowed to reach a static equilibrium state under the action of gravity. The test was repeated using DEM and the angle of repose measured as shown in Fig. 8. It is also known that for frictional cohesionless granular material, the angle of repose is a good indication of the internal friction angle [27]. 3.6. Summary of corn properties Table 1 summarises the corn material properties. The measured values were determined experimentally. The DEM macro-properties were obtained numerically using the calibrated micro-properties. The values of the angle of repose, average density and corn-steel friction coefficient correspond well to the values measured by Reimbert and Reimbert [28].

1.8

4. Blade experimental setup

Experiment

1.6

According to Maciejewski et al. [29], in practical cases when the three-dimensional motion of a bucket or bull-

1.4 E’ [MPa] = 0.0034 kn/s [kN/m]

E’ [MPa]

1.2 1.0 0.8 0.6

μ= 0.10 μ= 0.12

0.4

μ= 0.20

0.2

0

μ= 0.30 50

100

150 200 250 300 350 Particle stiffness kn = ks [kN/m]

400

450

500

Fig. 7. Effect of particle stiffness and particle friction on the confined stiffness.

Fig. 8. Angle of repose.

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26 Table 1 Measured corn properties and calibrated DEM parameters.

2nd Arm

Macro-property

Measured

DEM

Reimbert and Reimbert [28]

Internal friction angle, / Angle of repose, /r Average density, q Confined bulk modulus, E0 Material-steel friction, /s Material-glass friction, /g

23°

24°

26–29°

25 ± 2° 778 kg m3 1.60 MPa

24 ± 1° 778 kg m3 1.52 MPa

– 780–820 kg m3 –

14°

14°

17–23°

12°

12°

–

21

1st Arm

Direction of motion

Force transducer α

Calibrated micro-properties Particle stiffness, kn = ks Particle density, qp Particle friction coefficient, l

450 kN/m 855 kg/m3 0.12 Fig. 9. Blade experimental setup.

Other properties Damping, C Model width

0.7 0.2 m

dozer blade is discussed, plane strain conditions apply only in some deformation regions. The plane strain solution for such tools can be assumed only with limited accuracy. Maciejewski et al. [29] also investigated the assumption of plane strain conditions in soil bins where the soil and tool motion is constrained between two transparent walls. For measurements in such a bin, the forces acting on the tool due to the friction between the soil and the sidewalls has to be estimated or neglected. In this paper the assumption of two-dimensional deformation is made and the frictional forces against the side panels ignored in all DEM simulations. Two parallel glass panels were fixed 200 mm apart to form the soil bin. The blade profile was fixed to a trolley. The trolley moved on two linear bearing and was driven by a ball screw and stepper motor. The arm mechanism (first arm and second arm in Fig. 9) was used to set the blade depth and orientation. The initial blade depth (immersed depth into the material) could be set by changing and fixing the position of the second arm. The blade rake angle a could be set independently as indicated. Spring loaded Teflon wipers were used to seal the small opening between the blade and the glass panels. A force transducer was designed and built to measure the drag force, vertical force and moment acting on the blade. Three sets of strain gauges were bonded to a steel beam of which the position is shown in Fig. 9. The one set of two strain gauges were used to measure the moment acting on the beam, another set of four strain gauges were used to measure the horizontal force independent of the application point and the last set of strain gauges was used to measure the vertical force. The three sets of strain gauges were calibrated and the calibration checked regularly to avoid drift in the measurements.

The effect of tool velocity on draft force was investigated within the velocity range 1–100 mm s1. The velocity could be accurately set and controlled using the stepper motor. The results indicated that within this range, the velocity had little to no effect on the measured draft force [25]. This is in agreement with experimental results by Albert et al. [30]. Bohatier and Nouguier [31] have shown that within the velocity range used, the flow can be expected to be quasi-static in nature and hence the draft force is independent of the velocity. Using the equations derived by Siemens et al. [32], it can be shown that using the material, blade size and velocities presented in this paper, the dynamic force (acceleration) would account for less than 1% of the total draft force. A velocity of 10 mm s1 was used in all the experiments presented in this paper and quasi-static conditions assumed. In the DEM model the same blade velocity, size of the blade and the width of the test rig were used as in the experiments. 5. Results and discussion The calibrated micro-properties were used in blade–corn DEM simulations and the results compared to experimental measurements and observations. Fig. 10 shows the blade at an initial depth of h = 200 mm and the rake angle a = 90°. Snapshots were taken at displacement increments of 100 mm. The predicted shape of the heap that forms in front of the blade corresponds well to the experiment. The predicted heap is however, somewhat smaller. At a displacement of 400 mm, the maximum measured vertical displacement is 175 mm while the model predicts 145 mm (17% error). This can be attributed to a difference in the initial material porosity. The porosity was not included in the calibration method and needs further investigation. Fig. 11 shows the measured and simulated draft force on the blade with the initial depth h = 200 mm and the rake angle a = 90°. The simulated force follows the measured

22

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

force closely up to a displacement of 120 mm. During further displacement the model, however, predicts a draft force lower than the measured. At a displacement of 400 mm, the error is roughly 26%. Fig. 12 shows the vertical force on the blade. The measured force shows large fluctuations which makes a direct comparison difficult. The predicted force, however, shows the same general trend. Fig. 13 shows the moment acting on the blade. The moment is taken around the point where the force transducer is mounted (Fig. 9). Again, the simu-

lation predicts a moment lower than the measured moment, but the general trend is captured. Taking the total blade force and the moment, the point (moment arm) where the resultant force acts on the blade can be calculated. Fig. 14 shows the force application point measured from the bottom edge of the blade and expressed as a percentage of the total blade height. The application point shifts upward during the experiment/simulation as material builds up in front of the blade. The application point was found to be always between 30% and 50% of

Fig. 10. The blade moving through the material.

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26 350

23

150

300

Experiment

Experiment 250

Moment [Nm]

Draft force [N]

100 200

150

Simulation

Simulation 50

100 50

0 50

100

150 200 250 Displacement [mm]

300

350

0

400

50

0

-10

Vertical force [N]

Simulation -20

-30 Experiment -40

-50

100

150 200 250 Displacement [mm]

Fig. 12. Blade vertical force.

300

350

400

300

350

400

40 35 Force application height [%]

the immersed or covered blade depth. According to plasticity theory the pressure against the blade is a linear function of the depth which means that the force application point is equal to one-third of the blade depth, measured from the bottom edge [11]. Analytical methods do not take the up heave of material into account in predicting the blade forces. The effect of the heap can only be incorporated if the exact shape is known in advance. In order to make comparisons, the experiments were repeated with different blade depths (100, 150, 200, 250, 300, 350 mm) and rake angles (60°, 70°, 80°, 90°). The blade was moved forward at 0.5 mm s1 and the excess material (up heave) carefully removed. This procedure is similar to that used by Osman [6]. The same procedure was used during simulations where all particles above a certain height were simply deleted from the model. The draft forces were recorded and compared to that predicted by plasticity theory. The force on the blade can be calculated using plasticity theory and the soil stresses in front of the blade. Using the method of characteristics (finite difference method) nodes

50

150 200 250 Displacement [mm]

Fig. 13. Blade moment.

Fig. 11. Blade draft force.

0

100

Experiment

30 25

Simulation

20 15 Resultant force

10

Application height

5 0

50

100

150 250 200 Displacement [mm]

300

350

400

Fig. 14. Blade force application point.

are placed along the length of the blade. Once a solution is obtained, the stresses along the blade can be integrated over the blade area to get the total force acting on the blade. The results are summarised in Fig. 15 with the linear correlation equation fitted to each data set shown. Most of the forces, calculated using plasticity theory, are lower than the equivalent measured and DEM simulated forces. Perumpral et al. [10], Swick and Perumpral [7] and Osman [6] all reported the theoretical predictions to be lower than the measured values of draft force. The side panel friction is expected to increase the measured forces. The abovementioned authors, however, did not make use of a 2D rig with side panels and still found that the theoretical models underestimated the draft force. A further investigation showed that the accuracy of the analytical method decreased with a decrease in the rake angle a. This is also reported by Osman [6]. Fig. 16 shows the comparison between the DEM simulated and measured draft forces. There is a good agreement for measured draft forces lower than 150 N. With measured draft forces higher than

24

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26 400 y = 1.30 ⋅ x

y = 1.15 ⋅ x

y = 1.13 ⋅ x

350

350 Experimental draft force [N]

Experimental and simulated draft force [N]

400

300 250 200 150 100

300 250 200 150 100

Experiment 50

DEM Simulation

50 0 50

100

150

200

250

300

350

Plasticity theory draft force [N]

Fig. 15. Blade draft force comparison: experiment, DEM and plasticity theory.

150 N, DEM predicts forces lower than the measured forces. The reason for this might be that when the blade depth is higher, the amount of material and hence the frictional force on the side panels is greater. Since the frictional force was not included in the simulation, it was expected to predict lower forces compared to the measured forces. Soil failure is usually characterised by specific failure mechanisms. Definite shear lines or rupture lines can be identified. The position and shape of the shear lines depend on the soil properties and the tool geometry. The shear lines also change position and shape as the tool is displaced. The ability of DEM to predict the shear lines was investigated. The particle displacement ratio (PDR) is defined as the ratio of the magnitude of the particle absolute displacement vector to the magnitude of the blade displacement vector. The particles were then coloured according to their individual PDR values. The shear lines can also be determined from plasticity theory. Solving the hyperbolic partial differential equation using the method of characteristics, the so-called character-

0

50

100

150

200

250

300

350

DEM draft force [N]

Fig. 16. Blade draft force comparison: experiment and DEM.

istic curves are obtained. The most outer of these curves is the shear line. Fig. 17 shows the shear line for a rake angle a = 90°. A PDR P 0.15 was found to result in a shear line with its shape and position close to the observed shear line and the shear line predicted by plasticity theory. In Fig. 17 all the dark particles have a PDR P 0.15. During the experiment, the flow of particles could be observed and with the blade velocity low, the shear line could be traced onto the glass pane. This method has some degree of error, but the general position and shape of the shear line is indicated in Fig. 17. Using plasticity theory, the shear line is a direct result of the calculations. In Fig. 17 the total shear zone is slightly larger than the experimental one, the shape of the shear lines is very similar. The shear line predicted by plasticity theory seems to be less curved compared to the observed and simulated lines. The analytical approach does not include the size of the individual particles. In reality the particles have a finite size and this influences the position and width of the shear line, especially close to the lower edge of the blade where the theory assumes the shear line

Fig. 17. Shear line using a particle displacement ratio PDR = 0.15 and rake angle a = 90°.

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

25

Fig. 18. Shear line using relative particle displacements of 2% and rake angle a = 60°.

to start at the edge. In reality, the shear line does not start at the edge due to the finite size of the particles. A different approach to identify the shear lines is to make use of relative particle displacements. The blade was given a displacement of 150 mm to make sure failure has occurred. Using a ‘‘particle-child” data structure, the position of all the particles and their neighbouring particles was recorded and stored. The blade was then given a further displacement of 10 mm and the new position of each particle recorded. The new distance between the centres of two original parent- and child-particles was then determined and compared to the original centre distance. If the centre distance increased by 2% or more, the colour of both particles was changed. The result is shown in Fig. 18 for a rake angle a = 60°. Again good correlation between the observed, simulated and analytical results was achieved. The shear line became a straight line with the decrease in rake angle. The theoretical angle of this line based on analytical models is equal to (45°  //2) = 34° while the modelled angle is roughly 32°. From Figs. 17 and 18 it can be seen that the first method predicts a shear line which is much smoother and continuous when compared to the shear line predicted by the second method. This was observed for all rake angles tested. The second method, however, is more realistic when keeping in mind that a shear line is formed by a velocity discontinuity (strong or weak) between material layers. 6. Conclusion The main objective of this paper was to develop and validate a calibration process to determine the micro-material parameters needed in DEM simulations. The study was limited to cohesionless granular material and two-dimensional models. Relatively large particles were used which allowed the real size of the particles to be modelled. The calibration process was based on experimental and DEM shear tests and compression tests.

The conclusions of the paper are: 1. The internal friction angle resulting from DEM shear tests depends on both the particle stiffness value and the particle friction coefficient. From these results alone a unique set of stiffness and friction values cannot be determined. The results from the compression test show that the simulated macro-stiffness is a linear function of the particle stiffness and not influenced by the particle friction coefficient. Thus, from the compression test, the particle stiffness can be determined and then from the shear test results the particle friction coefficient can be determined. This results in a unique set of parameters which yields macro-values close to the measured values. 2. Using the calibrated material parameters, a flat blade moving through corn grains were modelled and the results compared to experimental results. Results show that during the initial stages of blade displacement, DEM can accurately predict the blade forces and moments. During further displacements and subsequent up heave of material in front of the blade, DEM however predicts forces/moments lower than the measured forces/moments. At the end of the blade displacement the error is 26%. 3. The up heave of material was prevented to make comparisons with plasticity theory. In general, plasticity theory predicts draft forces lower than the measured forces with a maximum error of 33%. With no up heave, there is a good correlation between the DEM simulated and experimental draft force with a maximum error of 22%. 4. Two methods were used to predict the shear lines. Results from both methods are in good agreement with experimental observations. The first method gives a much smoother and continuous shear line compared to the second method. 5. Further study is needed to include soils such as sand and cohesive materials such as clay. In this study, corn grains were used for ease of experiments and DEM simulations. The next step would be to make use of cohesionless soil

26

C.J. Coetzee, D.N.J. Els / Journal of Terramechanics 46 (2009) 15–26

such as sand. Thereafter cohesive soils such as clay should be modelled. When these types of materials are modelled, it is no longer practical to model the size of the physical particles accurately [33]. In DEM, larger particle sizes should be used and the influence of the particle size in combination with the other micro-parameters on the macro-behaviour needs to be studied. 6. In the proposed calibration method, the material porosity was not measured. Future methods should include the porosity in the calibration method so that there is a close agreement between the measured porosity and the modelled porosity. The particle size distribution and the accuracy in shape representation would influence the modelled porosity. Also, in 3D modelling the porosity can be modelled more accurately than in the 2D modelling used here.

[15]

[16]

[17] [18]

[19] [20] [21]

[22]

References [1] Franco Y, Rubinstein D, Shmulevich I. Determination of discrete element model parameters for soil–bulldozer blade interaction. In: Proceedings of the 15th international conference of the ISTVS, Hayama, Japan, September 25–29, 2005. [2] Karmakar S, Kushwaha R. Dynamic modelling of soil–tool interaction: an overview from a fluid flow perspective. J Terramech 2006;43:411–25. [3] Dechao Z, Yusuf Y. A dynamic model for soil cutting by blade and tine. J Terramech 1992;29(3):317–27. [4] Ibarra SY, McKeys E, Broughton RS. A model of stress distribution and cracking in cohesive soils produced by simple tillage implements. J Terramech 2005;42:115–39. [5] Onwualu AP, Watts KC. Draught and vertical forces obtained from dynamic soil cutting by plane tillage tools. Soil Till Res 1998;48:239–53. [6] Osman MS. The mechanics of soil cutting blades. J Agric Eng Res 1964;9(4):313–28. [7] Swick WC, Perumpral JV. A model for predicting soil–tool interaction. J Terramech 1998;25(1):43–56. [8] Terzaghi K. Theoretical soil mechanics. New York: Wiley; 1943. [9] Hansen JB. Earth pressure calculation. Copenhagen: The Danish Technical Press; 1961. [10] Perumpral JV, Grisso RD, Desai CS. A soil–tool model based on limit equilibrium analysis. Trans ASAE 1983:991. [11] Sokolovski VV. Statics of soil media. London: Butterworth Scientific Publishers; 1954 [Translated 1960]. [12] Mu¨hlhaus H-B. Discrete and continuous models for cohesive frictional materials. In: International symposium on continuous and discontinuous modelling of cohesive frictional materials, Stuttgart, 27–28 April, 2000. [13] Asaf Z, Rubinstein D, Shmulevich I. Evaluation of link-track performances using DEM. J Terramech 2006;43:141–61. [14] Asaf Z, Rubinstein D, Shmulevich I. Determination of discrete element model parameters using in-situ tests and inverse solution

[23]

[24] [25] [26]

[27] [28] [29]

[30] [31]

[32] [33] [34]

[35]

[36]

techniques. In: Proceedings of the 15th international conference of the ISTVS, Hayama, Japan, September 25–29, 2005. Vu-Quoc L, Zhang X, Walton OR. A 3-D discrete-element method for dry granular flows of ellipsoidal particles. Comput Methods Appl Mech Eng 2000;187:483–528. Tanaka H, Momozu M, Oida A, Yamazaki M. Simulation of soil deformation and resistance at bar penetration by the distinct element method. J Terramech 2000;37:41–56. McKeys E. Soil cutting and tillage. Amsterdam: Elsevier; 1985. Asaf Z, Rubinstein D, Shmulevich I. Determination of discrete element model parameters required for soil tillage. Soil Till 2007;92:227–42. Cundall PA, Strack ODL. A discrete numerical method for granular assemblies. Geotechnique 1979;29:47–65. Itasca. PFC2D theory and background manual 1999, Version 2.0. Available from: http://www.itascacg.com. Cleary PW, Sawley ML. Three-dimensional modelling of industrial granular flows. In: Second international conference on CFD in the minerals and process industries, 1999. p. 95–100. Hogue C. Shape representation and contact detection for discrete element simulations of arbitrary geometries. Eng Comput 1998;15(3):374–90. Zhang D, Whiten WJ. The calculation of contact forces between particles using spring and damping models. Powder Technol 1996;88:59–64. Rowlands JC. Dragline bucket filling. Ph.D. Thesis, University of Queensland, Queensland, Australia; 1991. Coetzee CJ. Forced granular flow. M.Sc. Thesis, Mechanical Engineering, University of Stellenbosch, Stellenbosch, South Africa; 2000. Corriveau D, Savage SB, Oger L. Internal friction angles: characterization using biaxial test simulations. In: IUTAM symposium on mechanics of granular and porous materials. Netherlands: Kluwer Academic Publishers; 1996. Lambe TW, Whitman RV. Soil mechanics. New York: John Wiley & Sons Inc.; 1969. Reimbert ML, Reimbert AM. Silos theory and practice. Clausthal, Germany: Trans Tech Publications; 1976. Maciejewski J, Jarzebowski A, Trampczynski W. Study on the efficiency of the digging process using the model of excavator bucket. J Terramech 2004;40:221–33. Albert R, Pfeifer MA, Shiffer P, Barabasi A. Drag force in a granular medium. Unpublished 1998. Bohatier C, Nouguier C. Dynamic soil–tool interaction forces and flow state. In: International conference on applied mechanics, SACAM 2000, Durban, South Africa, January 2000. Siemens JC, Weber JA, Thornburn TH. Mechanics of soil as influenced by model tillage tools. Trans ASAE 1965;8:1–7. Zhang R, Jianqiao L. Simulation on mechanical behaviour of cohesive soil by distinct element method. J Terramech 2006;43:303–16. Ni Q, Powrie W, Zhang X, Harkness R. Effect of particles properties on soil behaviour: 3-D numerical modeling of shearbox tests. Geotechn Spec Pub 2000;96:58–70. Mishra BK, Murty CVR. On the determination of contact parameters for realistic DEM simulations of ball mills. Powder Technol 2001;115:290–7. Landry H, Lague C, Roberge M. Discrete element representation of manure products. Comput Electron Agric 2006;51:17–34.