- Email: [email protected]
Simulation of soil deformation and resistance at bar penetration by the Distinct Element Method
Journal of Terramechanics 37 (2000) 41±56
www.elsevier.com/locate/jterra
Simulation of soil deformation and resistance at bar penetration by the Distinct Element Method$ Hiroaki Tanaka a,*, Masatoshi Momozu b, Akira Oida b, Minoru Yamazaki b a
Laboratory of Agricultural Machinery, Department of Hilly Land Agriculture, Shikoku National Agricultural Experiment Station, 2575, Ikano, Zentsuji, Kagawa 765, Japan b Agricultural Systems Engineering Section, Graduate School of Agriculture, Kyoto University, Kyoto 606-01, Japan
Abstract A mechanical model of soil is constructed using the Distinct Element Method (DEM) which makes it possible to analyze the discontinuous property of soil. To discuss the applicability of the soil model by the DEM, a bar penetration test was conducted and the result was compared with the simulation results. From the results of the behavior of elements, it could be said that the mechanical model by the DEM could well simulate the discontinuous behavior of soil and the parameters used in the simulation play important roles to make the soil model useful. As for the penetrating resistance, some problems which lie in the present DEM model are discussed and the key to solving these problems is indicated. Moreover, the method to determine the time interval used in the DEM simulations is mentioned in terms of the stability of the solution in the calculation. # 2000 ISTVS. All rights reserved.
1. Introduction There are many interactive phenomena between soil and farm machines. To investigate mechanical interactions between soil and machines is very important for designing farm machines or evaluating tillage methods. In terms of computational
$
This report was awarded ``best paper'' at the 12th International Conference of the ISTVS (Beijing, 1996). * Corresponding author. Tel.: +81-877-62-0800; fax: +81-877-62-1130. E-mail address: [email protected]rc.go.jp (H. Tanaka). 0022-4898/00/$20.00 # 2000 ISTVS. All rights reserved. PII: S0022-4898(99)00013-0
42
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Nomenclature dn t ; ds t en t ; en tÿt es t ; es tÿt Et ; Etÿt fn t ; fs t g Ii knp ; ksp knw ; ksw Mi t mi ri ; rj rij t ui ; uj un ; us ui t ; ui tÿt u_ i t ; u_ i tÿt u i t vi ; vj vi t ; vi tÿt v_i t ; v_i tÿt vi t Xi t xi t ; xi tÿt Yi t yi t ; yi tÿt ij np ; sp nw ; sw 'i ; 'j 'i t ; 'i tÿt 'it ;'itÿt '_ i t ; '_ i tÿt ' i t p ; w
normal and tangential damping forces on an element at time t normal spring forces on an element at time t and t ÿ t tangential spring forces on an element at time t and t ÿ t energies of a system at time t and t ÿ t normal and tangential contact forces on an element at time t gravitative acceleration moment of inertia of the element i normal and tangential spring constants between elements normal and tangential spring constants between element and wall resultant moment on the element i at time t mass of the element i radii of the element i and j distance between the centers of the elements i and j time interval horizontal displacements of the elements i and j during t normal and tangential relative displacements of the element i to the element j during t horizontal displacements of the element i at time t and t ÿ t horizontal velocities of the element i at time t and t ÿ t horizontal acceleration of the element i at time t vertical displacements of the elements i and j during t vertical displacements of the element i at time t and t ÿ t vertical velocities of the element i at time t and t ÿ t vertical acceleration of the element i at time t horizontal resultant force on the element i at time t horizontal positions of the element i at time t and t ÿ t vertical resultant force on the element i at time t vertical positions of the element i at time t and t ÿ t angle that the line connecting the centers of the elements i and j makes with the X axis normal and tangential coecients of viscous damping between elements normal and tangential coecients of viscous damping between element and wall rotational displacements of the elements i and j during t rotational positions of the element i at time t and t ÿ t rotational displacements of the element i at time t and t ÿ t rotational velocities of the element i at time t and t ÿ t rotational acceleration of the element i at time t coecients of friction between elements and between element and wall
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
43
mechanics, there have been many reports which have dealt with soil±machine mechanical interactions but almost all of them treated soil analyzed as a continuum. Soil originally consists of granular assemblies which can be cut and separated by the action of machine parts. In such a case, the continuum model is hard to apply. In 1979, Cundall and Strack [1] developed a discrete model to analyze behaviors of granular assemblies. It is called the Distinct Element Method (DEM). The DEM is applied in several ®elds such as powder technology, civil engineering [2] and soil mechanics. The mechanical model by the DEM has the capability to analyze the discontinuous property of the materials which is usually dicult to be analyzed by the continuum model. However, at the same time, the DEM also has some problems that the concrete method to determine the parameter has not been established yet, much time should be taken for calculation and so on. This report will begin with the explanation of the theory for construction of the soil model by the DEM and lead to the discussion about the constructed DEM model of soil by comparing with the result of the experiment. If the soil model by the DEM is constructed with high accuracy, it could be applied to many mechanical and dynamic problems between soil and machines in the ®eld of Terramechanics. 2. Mechanical model of soil using the DEM 2.1. Expression of soil by circular elements The most characteristic feature in the DEM is to assume that the object analyzed is an assembly of elements. In this study, the soil is represented as the assembly of two dimensional circular elements as shown in Fig. 1. It becomes possible to analyze the behavior of the soil by tracing the motion of each element which is derived from the mechanical interaction between elements or between the element and wall of the container or penetrating bar.
Fig. 1. The soil model by circular elements.
44
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
2.2. Mechanical interaction between elements Each element receives contact forces from contacting elements and/or sides of the walls. The magnitude of the contact force is determined by the relative displacement and relative velocity of the element. In the DEM simulation, a spring with the spring constant k, a dashpot with the coecient of viscous damping and a slider with the coecient of friction are supposed to be between elements. The same assumption is put between the element and the wall. The mechanical relationships are shown in Fig. 2. 2.3. Judgement of contact The judgement of contact of elements is based on the geometrical relationship between the element and the element/wall. In other words, if the distance between the centers of two elements i and j is smaller than the sum of the radii of these two elements, it follows that the element i contacts with the element j. Also in the case of the contact between element i and the wall of the container or the bar, if the distance between the center of the element i and the side of the wall or bar is smaller than the radius of the element i, it follows that the element i contacts with the wall/bar. In this study, the side of the metal bar is treated as the ``wall'' in the calculation. 2.4. Principle of motion of elements In the DEM, the positions of each element are determined step-by-step at intervals of t s. The position and contact force on each element determine the motion of the element. Simulation is conducted by repeating the calculations of these two values. The concept of this algorithm is shown in Fig. 3. The details of these calculations are mentioned in later sections. 2.5. Calculation of the contact force on each element Contact forces which act on a contacting element are proportional to the relative displacement and velocity of the element to the other contacting element. In this section, the calculation of contact force is formulated. As shown in Fig. 4, two elements
Fig. 2. Mechanical relationships between the element and the element/wall.
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
45
Fig. 3. The concept of the simulation algorithm.
Fig. 4. The displacements of the contact elements (ui ; vi are positive, uj ; vj are negative in this case).
i and j which contact with each other are taken as an example. This ®gure shows the positions of the elements i and j at time t, respectively. The tracks of motion of each element during t are represented by arrows. The horizontal (the X direction) and vertical (the Y direction) displacements of the elements i and j during t are expressed by the notations ui ; vi ; uj ; vj , respectively. Here the angle and directions are de®ned, which concern the formulation in this section. As shown in Fig. 5(a), considering a line which connects the centers of the elements i, j, and taking the center of the element i as the vertex, the angle that the line makes with the X axis is de®ned as ij . The angle ij is taken as positive in the counterclockwise direction. The sign rule for relative displacement of the element i to the element j is shown in Fig. 5(b). In the normal direction, un is taken as positive for compression. In the tangential direction, us is taken as positive in the
46
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Fig. 5. De®nition of the angle ij [(a)] and the sign rule [(b), (c) and (d)].
clockwise direction around the element j. In Fig. 5(c), the sign rule for contacting force on the element i given by the element j is indicated. fn t and fs t are taken as positive in the directions opposite to un and us . As for the rotational displacement i and the moment Mi t , each value is taken as positive in the counterclockwise direction as shown in Fig. 5(d). In accordance with these de®nitions, the relative displacement of the element i to the element j in the normal direction is expressed by Eq. (1), and in the tangential direction by Eq. (2). ÿ ÿ un ui ÿ uj cos ij vi ÿ vj sin ij
1
ÿ ÿ us ÿ ui ÿ uj sin ij vi ÿ vj cos ij ri 'i rj 'j
2
If these two elements contact for the ®rst time at time t, relative displacements in each direction are expressed by Eqs. (3) and (4), respectively, which are equal to the overlap quantities of the two elements. ÿ u0n ri rj ÿ rij u0s us
u0n un
3 4
From the relative displacements during t, the contact forces on the element i at time t are obtained based on the mechanical relationships in Fig. 2. In the normal direction for example, the spring force at time t, which is proportional to the overlap quantities, is obtained by adding the increment of the spring force knp un to the spring force at the previous step en tÿt as follows: en t en tÿt knp un
5
Therefore, the normal contact force on the element i given by the element j is obtained by adding the force by the dashpot dn t , which is proportional to the
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
47
relative velocity of the element i to the element j during t, to the force by the spring en t as follows: ÿ np un fn t en t dn t en tÿt knp un t
6
In the same way, the contact force in the tangential direction is calculated as follows: ÿ sp us fs t es t ds t es tÿt ksp us t
7
As shown by the no-tension joint in Fig. 2, the tensile force between contacting elements is not considered in this model, therefore: If en t < 0;
dn t e s t ds t 0
8
The magnitude of contact force by the spring in the tangential direction can not exceed the maximum force caused by friction, therefore: If
9 jes t j > p en t ; fs t es t p en t sig es t and ds t 0 where, sig es t is the sign of es t . Transforming the contact forces in the normal and tangential directions into the forces in the horizontal, vertical and rotational directions and summing up the contact forces caused by all of the elements which contact with the element i, resultant forces which act on the element i in the horizontal, vertical and rotational directions are obtained as follows: Xi t Yi t
Xÿ j
Xÿ j
Mi t ÿri
ÿ fn t cos ij fs t sin ij ÿ fn t sin ij ÿ fs t cos ij
Xÿ j
fs t
10 11 12
The same procedure is operated in the calculation between the element and wall. 2.6. Calculation of the position of each element The accelerations of the element i in the horizontal, vertical and rotational directions are obtained by dividing the contact forces in each direction [Eqs. (10)±(12)] by the mass or moment of inertia of the element i as follows:
48
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
u i t
Xi t mi
13
vi t
Yi t mi
14
' i t
Mi t Ii
15
Integrating the accelerations over the time interval t and adding the values to the previous velocities, the velocities of the element i in each direction at time t are obtained as follows: u_ i t u_ i tÿt u i t t
16
ÿ v_i t v_i tÿt vi t ÿg t
17
'_ i t '_ i tÿt ' i t t
18
The displacements of the element i in each direction at time t are given by integrating the velocities over the time interval t. However, the average of the displacements at time t and t ÿ t was adopted as the displacement at time t in actual calculation in order to get the stability of solution [2] as follows: ui t
1ÿ ui tÿt u_ i t t 2
19
vi t
1ÿ vi tÿt v_i t t 2
20
'i t
1ÿ 'i tÿt '_ i t t 2
21
Finally, adding the displacements of the element i at time t to the positions of the element i at time t ÿ t, the renewed positions of the element i at time t is obtained as follows: xi t xi tÿt ui t
22
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
49
yi t yi tÿt vi t
23
'i t 'i tÿt 'i t
24
Adopting these procedures to all of the elements in turn, the contacting forces, the displacements and the positions of all elements at time t are calculated. The values of the displacements ui t , vi t , 'i t calculated by Eqs. (19)±(21) are used as ui ; vi ; 'i in Eqs. (1) and (2) in the next step of the calculation to progress the simulation. 2.7. Stability of solution in the DEM simulation It is generally said that one of the DEM parameters eects so complicatedly the simulation result that the method to determine the DEM parameters is not clear. It is also said that the stability of solution is aected especially by the normal spring constant and the mass of element. In this study, the spring constant was set prior to the time interval and the value of the time interval was decided subsequently in order p to retain the stability of solution. The coecient p of viscous damping was set as 1= 2 times of the critical damping coecient 2 mk [3]. In this section, the interrelationship between normal spring constant k, mass of element m and time interval t is discussed and the limitation of the time interval is proposed. To simplify the problem, a certain system is assumed as shown in Fig. 6. This ®gure means that an element on the boundary of the wall at time t ÿ t moves perpendicularly to the wall during the time interval t, and overlaps with the wall at time t. The limitation of the time interval is to be considered in terms of the energy of the system at these two points in time t ÿ t and t. Assuming the displacement of element during t as un , the velocity of the element at t ÿ t is: u_ tÿt
un t
25
Fig. 6. Collision of an element to the wall.
50
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
At the point in time t ÿ t, as the element does not contact with the wall, the energy of the system at t ÿ t equals the kinetic energy of the element as follows: 1 Etÿt mu_ 2tÿt 2
26
After the time interval t, the element contacts with the wall and contact force generates. From the Eqs. (6),(10),(13) and (16), the velocity of the element at time t is:
un t u_ t u_ tÿt ÿ kun m t
27
At this point in time, the energy of the system equals the sum of the kinetic energy of the element and the potential energy of the spring as follows: 1 1 un un t 2 2 1 2 1 2 ÿ kun Et k un mu_ i k un m 2 2 2 2 m t t
28
Here, it is plain that the energy of the system is consumed by the dashpot during the time interval t. Accordingly, a certain condition is given to the energy of these two points in time as follows: Et < Etÿt
29
namely, 1 1 un un t 2 1 un 2 ÿ kun k un 2 m < m t t t 2 2 m 2
30
p In Eq. (30), the coecient of viscous damping can be replaced by 2mk as mentioned above. And Eq. (30) leads to a cubic inequality of t as follows: r pm32 p m 2 m <0 t t ÿ 2 2 t 2 2 k k k 3
31
The solution of this inequality within the range of 0 < t is: t < 0:75
r m k
32
In this way, the limitation of the time interval becomes clear as shown in Eq. (32). In actual calculation, some kinds of values are used as k and m because two kinds of normal spring constant (knp and knw ) are assumed in Fig. 2 and two kinds of
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
51
diameters of the elements are assumed in the simulation. Therefore, the least value of the limitations calculated by Eq. (32) should be taken as the limitation of the time interval in the calculation. 3. Experiment and simulation 3.1. Experimental procedure A laboratory experiment was conducted to con®rm the soil model by the DEM. As shown in Fig. 7, the soil was put into a container with the volume of 15151 cm3. The metal bar was penetrated vertically into the soil at the speed of 1 cm/s, and the resistant force which acts on the bar and the penetration depth of the bar were measured by the sensors. The behavior of the soil was photographed continuously. The properties of the soil are shown in Table 1. 3.2. Simulation procedure In the simulation of this experiment, the sizes of the container and the bar were set up the same as those in the above experiment. The diameters of the elements were determined to be 0.4 and 0.5 cm. The density of the element was determined to be 1.53 g/cm3 which is equal to that of soil in the experiment. These elements were located at random as the initial condition of the simulation. The number of elements
Fig. 7. The apparatus of the experiment.
52
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
was counted as 944. According to the parameters used in this simulation as shown in Table 2, the time interval t was set to 110ÿ4 s considering the stability of solution mentioned above. The outputs of calculation (the positions of elements and resistant force) were obtained every 1000 cycles of the calculation which equals 0.1 s of the simulation time. 4. Results and discussion 4.1. Result of the experiment (behavior of soil) The soil behavior at several points in time is shown in Fig. 8 as an example. As the soil was marked by powdered chalk every 3 cm from the bottom of the Table 1 Physical properties of soil used in the experiment Composition
Moisture content
Liquid limit
Plastic limit
Sand 57% Silt 31% Clay 12%
8.82%
37%
19.2%
Table 2 Parameters used in the simulation (cf. Fig. 2) Element±element knp ksp np sp p
1000 250 0.774 0.387 0.0, 0.25, 0.5, 0.75
Element±wall [N/m] [N/m] [N.s/m] [N.s/m] [ÿ]
knw ksw nw sw w
1000 250 0.774 0.387 0.0, 0.25, 0.50, 0.75
Fig. 8. The soil behavior in the experiment.
[N/m] [N/m] [N.s/m] [N.s/m] [ÿ]
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
53
container, it can be seen that the soil near the bar moved downward following the movement of the bar. 4.2. Result of the simulation (behavior of elements) Some results of the simulation are shown in Fig. 9. The elements near the bar moved downward following the movement of the bar. It can be said that the tendency of the behavior of elements was the same as that of the soil in the experiment. 4.3. The eect of the coecient of friction between the element and wall In this study, as the ®rst step to build up the concrete method to determine the parameters, the eect of one of the parameters (i.e. the coecient of friction) was examined. Fig. 10 shows some results of the experiment. These results show the positions of the elements at the same elapsed time of 9.0 s when the coecient of friction between element±element was assumed to be 0.75. But the coecient of friction between the element and wall was varied as shown in the ®gure. The eect of the coecient of friction clearly appeared in the behavior of elements. It can be seen that as the coecient of friction becomes higher, the movement
Fig. 9. Some results of the simulation.
Fig. 10. The eect of the coecient of friction between the element and wall.
54
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
of the elements near the bar comes larger. It also can be seen that as the coecient of friction becomes higher, a surface bump appeared clearer. Comparing these results with those of the experiment, the simulation with the coecient of friction between the element and wall of 0.25 can best simulate the behavior of the soil in the experiment. Therefore, it could be con®rmed that the value of parameters have essential eects on the behavior of the elements and the accuracy of the simulation becomes higher if the value of each parameter is selected appropriately. 4.4. The result of penetration resistance (soil and simulation) The penetration resistance in the experiment and the simulation are shown in Fig. 11. Each graph shows the relationship between the penetration depth and the resistance. The tendency that the penetration resistance becomes larger according to the penetration depth is shown in both results of the experiment and simulation, but the dierence between experiment and simulation appeared. The penetration resistance in the simulation considerably ¯uctuated. One reason could be that the diameter of the elements is large, and the other reason could be that the model is lacking the eect of cohesion because of the no-tension joint (cf. Fig. 2) which is put in the mechanical model, so that the frequent change of the resistance given by one element would cause a big change of the penetration resistance. It also can be seen that the value of penetration resistance in the simulation is apparently smaller than that of the experiment. One reason could be that the value of normal spring constant was not appropriate, and the other reason is possibly that the interlocking between elements is dicult to be achieved by using the circular elements which are unstable and are apt to be stirred by the penetrating bar through the calculation. These characteristic behaviors are the problem with establishing the accurate soil model by the DEM. These problems are to be discussed more in the next section. 4.5. Additional experiment using rigid balls and discussion As mentioned in the previous section, there is a large dierence between the result of penetration resistance of soil and that of simulation. Therefore, it is necessary to discuss the reason for such a dierence. In order to con®rm the mechanical property of the DEM model used in this study, an additional experiment was conducted using rigid balls which are made of alumina (Fig. 12). Because the alumina ball is perfectly discontinuous, such a mechanical condition seems to be similar to that of the simulation. The penetrating resistance in the case of the rigid ball is shown in Fig. 13. The tendency of ¯uctuation in the rigid ball qualitatively agreed with that obtained by the simulation. Therefore, it follows that the discontinuous property of the present DEM model is very strong and this result would justify the supposition in the previous section that the ¯uctuation is caused by the large diameter and no-tension joint. The elements with adequate diameter and cohesive eect would be essential to settle the ¯uctuation of force in order to predict the force between soil and machine parts.
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Fig. 11. The penetration resistance of the experiment and the simulation.
Fig. 12. The additional experiment using the rigid balls.
Fig. 13. The result of penetration resistance in the case of rigid balls.
55
56
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
The magnitude of resistance in the simulation is still smaller than that in the rigid ball. Therefore, it follows that the representation of the eect of interlocking in the simulation is still insucient even in simulating the behavior of rigid balls. In a further study, this problem should be discussed in order to improve the mechanical model in the DEM. 5. Conclusions The mechanical model of soil was constructed by the Distinct Element Method, and the bar penetration test was simulated using the soil model. As mentioned above, the behavior of the soil was able to be well simulated by the movement of the elements. This paper also shows the eect of one of the parameters (the coecient of friction between the element and wall) by comparing some results of the behavior of elements. This discussion could contribute to the establishment of the concrete method to determine the values of the DEM parameters based on the soil properties. As shown in the result of the penetrating resistance, there lie some problems for the DEM to predict the force which acts on machine parts. This paper con®rmed the characteristic property of the present DEM model by way of the experiment using the rigid balls. The discussion oers a key to further improvement of the soil model by the DEM. Furthermore, this paper also contributes to the eciency of simulation works using the DEM since the limitation of the value of the time interval to get the stability of solution was derived as the practical problem in the calculation. The DEM would improve its ability to analyze the mechanical interaction between soil and machine if this method is discussed more and becomes sophisticated. Acknowledgements We wish to express our warmest gratitude to A. Professor Dr. Eiichiro Sakaguchi (Tokyo University of Agriculture) for his assistance in developing the DEM program used in this study. References [1] Cundall PA, Strack ODL. A discrete numerical method for granular assemblies. Geotechnique 1979;29(1):47±65. [2] Kiyama H, Fujimura H. Application of Cundall's discrete block method to gravity ¯ow analysis of rock-like granular materials (in Japanese). Proc JSCE 1983;333:137±45. [3] Sakaguchi, E. et al. Simulation on ¯owing phenomena of grains by distinct element method. Proc CIGR Ag Eng Milano '94, 1994, Rep. No. 94-G-025.
www.elsevier.com/locate/jterra
Simulation of soil deformation and resistance at bar penetration by the Distinct Element Method$ Hiroaki Tanaka a,*, Masatoshi Momozu b, Akira Oida b, Minoru Yamazaki b a
Laboratory of Agricultural Machinery, Department of Hilly Land Agriculture, Shikoku National Agricultural Experiment Station, 2575, Ikano, Zentsuji, Kagawa 765, Japan b Agricultural Systems Engineering Section, Graduate School of Agriculture, Kyoto University, Kyoto 606-01, Japan
Abstract A mechanical model of soil is constructed using the Distinct Element Method (DEM) which makes it possible to analyze the discontinuous property of soil. To discuss the applicability of the soil model by the DEM, a bar penetration test was conducted and the result was compared with the simulation results. From the results of the behavior of elements, it could be said that the mechanical model by the DEM could well simulate the discontinuous behavior of soil and the parameters used in the simulation play important roles to make the soil model useful. As for the penetrating resistance, some problems which lie in the present DEM model are discussed and the key to solving these problems is indicated. Moreover, the method to determine the time interval used in the DEM simulations is mentioned in terms of the stability of the solution in the calculation. # 2000 ISTVS. All rights reserved.
1. Introduction There are many interactive phenomena between soil and farm machines. To investigate mechanical interactions between soil and machines is very important for designing farm machines or evaluating tillage methods. In terms of computational
$
This report was awarded ``best paper'' at the 12th International Conference of the ISTVS (Beijing, 1996). * Corresponding author. Tel.: +81-877-62-0800; fax: +81-877-62-1130. E-mail address: [email protected]rc.go.jp (H. Tanaka). 0022-4898/00/$20.00 # 2000 ISTVS. All rights reserved. PII: S0022-4898(99)00013-0
42
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Nomenclature dn t ; ds t en t ; en tÿt es t ; es tÿt Et ; Etÿt fn t ; fs t g Ii knp ; ksp knw ; ksw Mi t mi ri ; rj rij t ui ; uj un ; us ui t ; ui tÿt u_ i t ; u_ i tÿt u i t vi ; vj vi t ; vi tÿt v_i t ; v_i tÿt vi t Xi t xi t ; xi tÿt Yi t yi t ; yi tÿt ij np ; sp nw ; sw 'i ; 'j 'i t ; 'i tÿt 'it ;'itÿt '_ i t ; '_ i tÿt ' i t p ; w
normal and tangential damping forces on an element at time t normal spring forces on an element at time t and t ÿ t tangential spring forces on an element at time t and t ÿ t energies of a system at time t and t ÿ t normal and tangential contact forces on an element at time t gravitative acceleration moment of inertia of the element i normal and tangential spring constants between elements normal and tangential spring constants between element and wall resultant moment on the element i at time t mass of the element i radii of the element i and j distance between the centers of the elements i and j time interval horizontal displacements of the elements i and j during t normal and tangential relative displacements of the element i to the element j during t horizontal displacements of the element i at time t and t ÿ t horizontal velocities of the element i at time t and t ÿ t horizontal acceleration of the element i at time t vertical displacements of the elements i and j during t vertical displacements of the element i at time t and t ÿ t vertical velocities of the element i at time t and t ÿ t vertical acceleration of the element i at time t horizontal resultant force on the element i at time t horizontal positions of the element i at time t and t ÿ t vertical resultant force on the element i at time t vertical positions of the element i at time t and t ÿ t angle that the line connecting the centers of the elements i and j makes with the X axis normal and tangential coecients of viscous damping between elements normal and tangential coecients of viscous damping between element and wall rotational displacements of the elements i and j during t rotational positions of the element i at time t and t ÿ t rotational displacements of the element i at time t and t ÿ t rotational velocities of the element i at time t and t ÿ t rotational acceleration of the element i at time t coecients of friction between elements and between element and wall
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
43
mechanics, there have been many reports which have dealt with soil±machine mechanical interactions but almost all of them treated soil analyzed as a continuum. Soil originally consists of granular assemblies which can be cut and separated by the action of machine parts. In such a case, the continuum model is hard to apply. In 1979, Cundall and Strack [1] developed a discrete model to analyze behaviors of granular assemblies. It is called the Distinct Element Method (DEM). The DEM is applied in several ®elds such as powder technology, civil engineering [2] and soil mechanics. The mechanical model by the DEM has the capability to analyze the discontinuous property of the materials which is usually dicult to be analyzed by the continuum model. However, at the same time, the DEM also has some problems that the concrete method to determine the parameter has not been established yet, much time should be taken for calculation and so on. This report will begin with the explanation of the theory for construction of the soil model by the DEM and lead to the discussion about the constructed DEM model of soil by comparing with the result of the experiment. If the soil model by the DEM is constructed with high accuracy, it could be applied to many mechanical and dynamic problems between soil and machines in the ®eld of Terramechanics. 2. Mechanical model of soil using the DEM 2.1. Expression of soil by circular elements The most characteristic feature in the DEM is to assume that the object analyzed is an assembly of elements. In this study, the soil is represented as the assembly of two dimensional circular elements as shown in Fig. 1. It becomes possible to analyze the behavior of the soil by tracing the motion of each element which is derived from the mechanical interaction between elements or between the element and wall of the container or penetrating bar.
Fig. 1. The soil model by circular elements.
44
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
2.2. Mechanical interaction between elements Each element receives contact forces from contacting elements and/or sides of the walls. The magnitude of the contact force is determined by the relative displacement and relative velocity of the element. In the DEM simulation, a spring with the spring constant k, a dashpot with the coecient of viscous damping and a slider with the coecient of friction are supposed to be between elements. The same assumption is put between the element and the wall. The mechanical relationships are shown in Fig. 2. 2.3. Judgement of contact The judgement of contact of elements is based on the geometrical relationship between the element and the element/wall. In other words, if the distance between the centers of two elements i and j is smaller than the sum of the radii of these two elements, it follows that the element i contacts with the element j. Also in the case of the contact between element i and the wall of the container or the bar, if the distance between the center of the element i and the side of the wall or bar is smaller than the radius of the element i, it follows that the element i contacts with the wall/bar. In this study, the side of the metal bar is treated as the ``wall'' in the calculation. 2.4. Principle of motion of elements In the DEM, the positions of each element are determined step-by-step at intervals of t s. The position and contact force on each element determine the motion of the element. Simulation is conducted by repeating the calculations of these two values. The concept of this algorithm is shown in Fig. 3. The details of these calculations are mentioned in later sections. 2.5. Calculation of the contact force on each element Contact forces which act on a contacting element are proportional to the relative displacement and velocity of the element to the other contacting element. In this section, the calculation of contact force is formulated. As shown in Fig. 4, two elements
Fig. 2. Mechanical relationships between the element and the element/wall.
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
45
Fig. 3. The concept of the simulation algorithm.
Fig. 4. The displacements of the contact elements (ui ; vi are positive, uj ; vj are negative in this case).
i and j which contact with each other are taken as an example. This ®gure shows the positions of the elements i and j at time t, respectively. The tracks of motion of each element during t are represented by arrows. The horizontal (the X direction) and vertical (the Y direction) displacements of the elements i and j during t are expressed by the notations ui ; vi ; uj ; vj , respectively. Here the angle and directions are de®ned, which concern the formulation in this section. As shown in Fig. 5(a), considering a line which connects the centers of the elements i, j, and taking the center of the element i as the vertex, the angle that the line makes with the X axis is de®ned as ij . The angle ij is taken as positive in the counterclockwise direction. The sign rule for relative displacement of the element i to the element j is shown in Fig. 5(b). In the normal direction, un is taken as positive for compression. In the tangential direction, us is taken as positive in the
46
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Fig. 5. De®nition of the angle ij [(a)] and the sign rule [(b), (c) and (d)].
clockwise direction around the element j. In Fig. 5(c), the sign rule for contacting force on the element i given by the element j is indicated. fn t and fs t are taken as positive in the directions opposite to un and us . As for the rotational displacement i and the moment Mi t , each value is taken as positive in the counterclockwise direction as shown in Fig. 5(d). In accordance with these de®nitions, the relative displacement of the element i to the element j in the normal direction is expressed by Eq. (1), and in the tangential direction by Eq. (2). ÿ ÿ un ui ÿ uj cos ij vi ÿ vj sin ij
1
ÿ ÿ us ÿ ui ÿ uj sin ij vi ÿ vj cos ij ri 'i rj 'j
2
If these two elements contact for the ®rst time at time t, relative displacements in each direction are expressed by Eqs. (3) and (4), respectively, which are equal to the overlap quantities of the two elements. ÿ u0n ri rj ÿ rij u0s us
u0n un
3 4
From the relative displacements during t, the contact forces on the element i at time t are obtained based on the mechanical relationships in Fig. 2. In the normal direction for example, the spring force at time t, which is proportional to the overlap quantities, is obtained by adding the increment of the spring force knp un to the spring force at the previous step en tÿt as follows: en t en tÿt knp un
5
Therefore, the normal contact force on the element i given by the element j is obtained by adding the force by the dashpot dn t , which is proportional to the
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
47
relative velocity of the element i to the element j during t, to the force by the spring en t as follows: ÿ np un fn t en t dn t en tÿt knp un t
6
In the same way, the contact force in the tangential direction is calculated as follows: ÿ sp us fs t es t ds t es tÿt ksp us t
7
As shown by the no-tension joint in Fig. 2, the tensile force between contacting elements is not considered in this model, therefore: If en t < 0;
dn t e s t ds t 0
8
The magnitude of contact force by the spring in the tangential direction can not exceed the maximum force caused by friction, therefore: If
9 jes t j > p en t ; fs t es t p en t sig es t and ds t 0 where, sig es t is the sign of es t . Transforming the contact forces in the normal and tangential directions into the forces in the horizontal, vertical and rotational directions and summing up the contact forces caused by all of the elements which contact with the element i, resultant forces which act on the element i in the horizontal, vertical and rotational directions are obtained as follows: Xi t Yi t
Xÿ j
Xÿ j
Mi t ÿri
ÿ fn t cos ij fs t sin ij ÿ fn t sin ij ÿ fs t cos ij
Xÿ j
fs t
10 11 12
The same procedure is operated in the calculation between the element and wall. 2.6. Calculation of the position of each element The accelerations of the element i in the horizontal, vertical and rotational directions are obtained by dividing the contact forces in each direction [Eqs. (10)±(12)] by the mass or moment of inertia of the element i as follows:
48
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
u i t
Xi t mi
13
vi t
Yi t mi
14
' i t
Mi t Ii
15
Integrating the accelerations over the time interval t and adding the values to the previous velocities, the velocities of the element i in each direction at time t are obtained as follows: u_ i t u_ i tÿt u i t t
16
ÿ v_i t v_i tÿt vi t ÿg t
17
'_ i t '_ i tÿt ' i t t
18
The displacements of the element i in each direction at time t are given by integrating the velocities over the time interval t. However, the average of the displacements at time t and t ÿ t was adopted as the displacement at time t in actual calculation in order to get the stability of solution [2] as follows: ui t
1ÿ ui tÿt u_ i t t 2
19
vi t
1ÿ vi tÿt v_i t t 2
20
'i t
1ÿ 'i tÿt '_ i t t 2
21
Finally, adding the displacements of the element i at time t to the positions of the element i at time t ÿ t, the renewed positions of the element i at time t is obtained as follows: xi t xi tÿt ui t
22
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
49
yi t yi tÿt vi t
23
'i t 'i tÿt 'i t
24
Adopting these procedures to all of the elements in turn, the contacting forces, the displacements and the positions of all elements at time t are calculated. The values of the displacements ui t , vi t , 'i t calculated by Eqs. (19)±(21) are used as ui ; vi ; 'i in Eqs. (1) and (2) in the next step of the calculation to progress the simulation. 2.7. Stability of solution in the DEM simulation It is generally said that one of the DEM parameters eects so complicatedly the simulation result that the method to determine the DEM parameters is not clear. It is also said that the stability of solution is aected especially by the normal spring constant and the mass of element. In this study, the spring constant was set prior to the time interval and the value of the time interval was decided subsequently in order p to retain the stability of solution. The coecient p of viscous damping was set as 1= 2 times of the critical damping coecient 2 mk [3]. In this section, the interrelationship between normal spring constant k, mass of element m and time interval t is discussed and the limitation of the time interval is proposed. To simplify the problem, a certain system is assumed as shown in Fig. 6. This ®gure means that an element on the boundary of the wall at time t ÿ t moves perpendicularly to the wall during the time interval t, and overlaps with the wall at time t. The limitation of the time interval is to be considered in terms of the energy of the system at these two points in time t ÿ t and t. Assuming the displacement of element during t as un , the velocity of the element at t ÿ t is: u_ tÿt
un t
25
Fig. 6. Collision of an element to the wall.
50
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
At the point in time t ÿ t, as the element does not contact with the wall, the energy of the system at t ÿ t equals the kinetic energy of the element as follows: 1 Etÿt mu_ 2tÿt 2
26
After the time interval t, the element contacts with the wall and contact force generates. From the Eqs. (6),(10),(13) and (16), the velocity of the element at time t is:
un t u_ t u_ tÿt ÿ kun m t
27
At this point in time, the energy of the system equals the sum of the kinetic energy of the element and the potential energy of the spring as follows: 1 1 un un t 2 2 1 2 1 2 ÿ kun Et k un mu_ i k un m 2 2 2 2 m t t
28
Here, it is plain that the energy of the system is consumed by the dashpot during the time interval t. Accordingly, a certain condition is given to the energy of these two points in time as follows: Et < Etÿt
29
namely, 1 1 un un t 2 1 un 2 ÿ kun k un 2 m < m t t t 2 2 m 2
30
p In Eq. (30), the coecient of viscous damping can be replaced by 2mk as mentioned above. And Eq. (30) leads to a cubic inequality of t as follows: r pm32 p m 2 m <0 t t ÿ 2 2 t 2 2 k k k 3
31
The solution of this inequality within the range of 0 < t is: t < 0:75
r m k
32
In this way, the limitation of the time interval becomes clear as shown in Eq. (32). In actual calculation, some kinds of values are used as k and m because two kinds of normal spring constant (knp and knw ) are assumed in Fig. 2 and two kinds of
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
51
diameters of the elements are assumed in the simulation. Therefore, the least value of the limitations calculated by Eq. (32) should be taken as the limitation of the time interval in the calculation. 3. Experiment and simulation 3.1. Experimental procedure A laboratory experiment was conducted to con®rm the soil model by the DEM. As shown in Fig. 7, the soil was put into a container with the volume of 15151 cm3. The metal bar was penetrated vertically into the soil at the speed of 1 cm/s, and the resistant force which acts on the bar and the penetration depth of the bar were measured by the sensors. The behavior of the soil was photographed continuously. The properties of the soil are shown in Table 1. 3.2. Simulation procedure In the simulation of this experiment, the sizes of the container and the bar were set up the same as those in the above experiment. The diameters of the elements were determined to be 0.4 and 0.5 cm. The density of the element was determined to be 1.53 g/cm3 which is equal to that of soil in the experiment. These elements were located at random as the initial condition of the simulation. The number of elements
Fig. 7. The apparatus of the experiment.
52
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
was counted as 944. According to the parameters used in this simulation as shown in Table 2, the time interval t was set to 110ÿ4 s considering the stability of solution mentioned above. The outputs of calculation (the positions of elements and resistant force) were obtained every 1000 cycles of the calculation which equals 0.1 s of the simulation time. 4. Results and discussion 4.1. Result of the experiment (behavior of soil) The soil behavior at several points in time is shown in Fig. 8 as an example. As the soil was marked by powdered chalk every 3 cm from the bottom of the Table 1 Physical properties of soil used in the experiment Composition
Moisture content
Liquid limit
Plastic limit
Sand 57% Silt 31% Clay 12%
8.82%
37%
19.2%
Table 2 Parameters used in the simulation (cf. Fig. 2) Element±element knp ksp np sp p
1000 250 0.774 0.387 0.0, 0.25, 0.5, 0.75
Element±wall [N/m] [N/m] [N.s/m] [N.s/m] [ÿ]
knw ksw nw sw w
1000 250 0.774 0.387 0.0, 0.25, 0.50, 0.75
Fig. 8. The soil behavior in the experiment.
[N/m] [N/m] [N.s/m] [N.s/m] [ÿ]
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
53
container, it can be seen that the soil near the bar moved downward following the movement of the bar. 4.2. Result of the simulation (behavior of elements) Some results of the simulation are shown in Fig. 9. The elements near the bar moved downward following the movement of the bar. It can be said that the tendency of the behavior of elements was the same as that of the soil in the experiment. 4.3. The eect of the coecient of friction between the element and wall In this study, as the ®rst step to build up the concrete method to determine the parameters, the eect of one of the parameters (i.e. the coecient of friction) was examined. Fig. 10 shows some results of the experiment. These results show the positions of the elements at the same elapsed time of 9.0 s when the coecient of friction between element±element was assumed to be 0.75. But the coecient of friction between the element and wall was varied as shown in the ®gure. The eect of the coecient of friction clearly appeared in the behavior of elements. It can be seen that as the coecient of friction becomes higher, the movement
Fig. 9. Some results of the simulation.
Fig. 10. The eect of the coecient of friction between the element and wall.
54
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
of the elements near the bar comes larger. It also can be seen that as the coecient of friction becomes higher, a surface bump appeared clearer. Comparing these results with those of the experiment, the simulation with the coecient of friction between the element and wall of 0.25 can best simulate the behavior of the soil in the experiment. Therefore, it could be con®rmed that the value of parameters have essential eects on the behavior of the elements and the accuracy of the simulation becomes higher if the value of each parameter is selected appropriately. 4.4. The result of penetration resistance (soil and simulation) The penetration resistance in the experiment and the simulation are shown in Fig. 11. Each graph shows the relationship between the penetration depth and the resistance. The tendency that the penetration resistance becomes larger according to the penetration depth is shown in both results of the experiment and simulation, but the dierence between experiment and simulation appeared. The penetration resistance in the simulation considerably ¯uctuated. One reason could be that the diameter of the elements is large, and the other reason could be that the model is lacking the eect of cohesion because of the no-tension joint (cf. Fig. 2) which is put in the mechanical model, so that the frequent change of the resistance given by one element would cause a big change of the penetration resistance. It also can be seen that the value of penetration resistance in the simulation is apparently smaller than that of the experiment. One reason could be that the value of normal spring constant was not appropriate, and the other reason is possibly that the interlocking between elements is dicult to be achieved by using the circular elements which are unstable and are apt to be stirred by the penetrating bar through the calculation. These characteristic behaviors are the problem with establishing the accurate soil model by the DEM. These problems are to be discussed more in the next section. 4.5. Additional experiment using rigid balls and discussion As mentioned in the previous section, there is a large dierence between the result of penetration resistance of soil and that of simulation. Therefore, it is necessary to discuss the reason for such a dierence. In order to con®rm the mechanical property of the DEM model used in this study, an additional experiment was conducted using rigid balls which are made of alumina (Fig. 12). Because the alumina ball is perfectly discontinuous, such a mechanical condition seems to be similar to that of the simulation. The penetrating resistance in the case of the rigid ball is shown in Fig. 13. The tendency of ¯uctuation in the rigid ball qualitatively agreed with that obtained by the simulation. Therefore, it follows that the discontinuous property of the present DEM model is very strong and this result would justify the supposition in the previous section that the ¯uctuation is caused by the large diameter and no-tension joint. The elements with adequate diameter and cohesive eect would be essential to settle the ¯uctuation of force in order to predict the force between soil and machine parts.
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
Fig. 11. The penetration resistance of the experiment and the simulation.
Fig. 12. The additional experiment using the rigid balls.
Fig. 13. The result of penetration resistance in the case of rigid balls.
55
56
H. Tanaka et al. / Journal of Terramechanics 37 (2000) 41±56
The magnitude of resistance in the simulation is still smaller than that in the rigid ball. Therefore, it follows that the representation of the eect of interlocking in the simulation is still insucient even in simulating the behavior of rigid balls. In a further study, this problem should be discussed in order to improve the mechanical model in the DEM. 5. Conclusions The mechanical model of soil was constructed by the Distinct Element Method, and the bar penetration test was simulated using the soil model. As mentioned above, the behavior of the soil was able to be well simulated by the movement of the elements. This paper also shows the eect of one of the parameters (the coecient of friction between the element and wall) by comparing some results of the behavior of elements. This discussion could contribute to the establishment of the concrete method to determine the values of the DEM parameters based on the soil properties. As shown in the result of the penetrating resistance, there lie some problems for the DEM to predict the force which acts on machine parts. This paper con®rmed the characteristic property of the present DEM model by way of the experiment using the rigid balls. The discussion oers a key to further improvement of the soil model by the DEM. Furthermore, this paper also contributes to the eciency of simulation works using the DEM since the limitation of the value of the time interval to get the stability of solution was derived as the practical problem in the calculation. The DEM would improve its ability to analyze the mechanical interaction between soil and machine if this method is discussed more and becomes sophisticated. Acknowledgements We wish to express our warmest gratitude to A. Professor Dr. Eiichiro Sakaguchi (Tokyo University of Agriculture) for his assistance in developing the DEM program used in this study. References [1] Cundall PA, Strack ODL. A discrete numerical method for granular assemblies. Geotechnique 1979;29(1):47±65. [2] Kiyama H, Fujimura H. Application of Cundall's discrete block method to gravity ¯ow analysis of rock-like granular materials (in Japanese). Proc JSCE 1983;333:137±45. [3] Sakaguchi, E. et al. Simulation on ¯owing phenomena of grains by distinct element method. Proc CIGR Ag Eng Milano '94, 1994, Rep. No. 94-G-025.