Investigation of normal contact interaction between two bonded spherical particles with interface layer/Dvieju per tarpini sluoksni kontaktuojanqu sferiniu daleliu normalinio kontakto saveikos tyrimas.

Pilkavicius, S. ; Kacianauskas, R. ; Norkus, A. 等

1. Introduction

Characterisation of mechanical behavior of heterogeneous solids and evaluation of specified parameters is important issue to be solved in earlier stage of technological and structural design. Many of natural solids like soils or rocks as well as industrially manufactured materials like concrete and resilient materials are actually particulate solids. They are composed by grains of different size, shape and mechanical properties.

One of the main building construction materials concrete presents characteristic sample of particulate solid, where the relevant structural formation of the material relies on particles packing and structure voids. Moreover concrete is aggregate which consists of several nature products such as granite and/or limestone, sand particles and matrix, which in turn consist of even smaller structures. Consequently, concrete is a bonded particulate material in which both the grains and the cement bonds are deformable and can break.

It is obvious that mechanical properties of the heterogeneous solids, i.e. macroscopic properties, are result of grain properties and their interactions. However, it is extremely difficult to measure the strength of bonds and especially tracks the performance of each bond in tests.

Recently numerical simulating presents feasible support to reduce experimental effort. Among the many proposed numerical simulation methods, the Discrete Element Method (DEM), introduced by Cundall [1] has become the most useful tool. In DEM simulation, the dynamic motion of each particle of media is monitored during whole process. This allows simulation of motion and interaction between the particles taking into account the microscopic geometry and various constitutive contact models. Evolution of the contact behaviour is decisive in the DEM. To save computational time, the DEM can operate by simplified description of the contact. Review of the earlier developments is presented by Dziugys and Peters [2]. Various models of the normal contact are presented by Kruggel-Emdem et al. [3], Tomas [4], Maknickas et al. [5].

The theoretical frame of normal contact of homogeneous spheres rises from the classical work of Hertz dating back 1881, where an analytical solution for the frictionless contact of two elastic spheres was derived. Details of the Hertz method can be found in the book of Johnson [6]. In this method the elastic contact behavior is clearly characterized by the force-displacement relationship containing the effective radius and effective elasticity modulus. Generally, even homogeneous spheres may be of different radii and of different materials. The influence of the differences in particle properties is illustrated in [7].

An extensive review of the literature on spherical and cylindrical contacts under normal load was made by Adams and Nosonovsky [8]. It was shown the most of the existing works on spherical contact were restricted to a perfect slip contact. Elastic solution of the normal contact under stick and sliding with the wide range of elastic constants are reviewed by Brizmer et al. [9]. Continuum-based solutions are extensively explored in development of DEM models.

Many studies have been investigating the micromechanical behaviour of particulate, porous and heterogeneous materials. For example, studies on cemented particulate materials by Dvorkin et al. [10] and Zhu et al. [11] provided information on the normal and tangential load transfer between cemented particles. Applications of such contact-based micromechanical analysis for asphalt mixture behaviour have been reported by Chang and Gao [12], Cheung et al. [13] and Zhu and Nodes [14].

Originally, most of the DEM applications are aimed to simulate non-cohesive granular materials with unilateral repulsive normal contact while real deformation of particles is replaced by their overlap. When considering bonded particles the approach applied is non-unique. A micro model for loosely packed bonded granule material with cemented by a finite-sized piece of cementitous material bonds was developed by Jiang et al. [15-17]. Here the bond is regarded as a rod (or system of rods), the section of which depend on particular problem.

The approach developed by Patyondy and Cundall [18], assumes normal bond stiffness implemented as stiffness of contacting particles stiffened by additional parallel link. This calculation model is employed and widely used in commercial DEM codes EDEM [19], PFC2D/3D discontinuum program [20] and in another numerous applications.

The paper presents investigation of normal contact interaction between two bonded spherical particles. Analytical model for normal contact of two spherical particles connected by interface bond of finite size is suggested. Various models were applied for examination of the bond having different radius and material properties. Validity of various analytical models was examined by applying the Finite Element method (FEM).

The paper is organised in the following way. DEM methodology and normal contact of bonded particles are presented in section 2. Investigation of the bond properties was done in section 3. FE simulation of normal interaction is described in section 4. Evaluation of analytical and FE model results was described in section 5, while concluding remarks are given in section 6.

2. Normal contact of bonded particles

2.1. Problem description

Particulate solids are regarded as a system of the finite number n of particles, characterised by the prescribed material properties and the constitutive interaction laws. The dynamical behaviour and deformation of the solid is considered by the motion of an individual particle i, (i = 1, 2, ... , n) by applying the Newton's second law. Translational motion of the center of gravity of particle i can be fully described by the system of ordinary differential equations

[m.sub.i] [d.sup.2][x.sub.i]/[dt.sup.2] = [m.sub.i][a.sub.i] = [F.sub.i] (1)

where [a.sub.i] and [x.sub.i] are vectors of acceleration and the position of the centre of gravity [m.sub.i] of the particle i, respectively. Vector [F.sub.i] is the sum of all forces acting on the centre of gravity of the particle.

If we restrict ourselves to contacts with j (j = 1, 2, ... , nc) neighbor particles, resultant force can be described as the sum of direct contact forces between the particle i and another particle j

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

The inter-particle contact force [F.sub.ij] can be expressed as the sum of normal and tangential components

[F.sub.ij] = [F.sub.n,ij] + [F.sub.t,ij] (3)

Regarding the contact, we focus to the normal component [F.sub.n,ij]. The normal contact force between particles can be expressed as the resultant force reflecting various interaction effects.

Normal contact of bonded particles comprises contact of two spherical particles bonded by the interface of finite size.

Geometry of two bonded deformable spheres i and j with the radii [R.sub.i] and [R.sub.j] is shown in Fig. 1, a.

[FIGURE 1 OMITTED]

The location of spheres is characterised by the central points [O.sub.i] and [O.sub.j] referring to local Cartesian coordinates x, y, z where axis Oz points the normal direction. The centres of spheres are defined by the coordinates [z.sub.i] and [Z.sub.j], respectively. The initial distance between surfaces of interacting spheres is denoted by [L.sub.c].

The spheres are assumed to be connected by the interface material of finite size termed usually a bond. Geometry of the bond is circular cylinder bounded by spheres surfaces. The bond radius [R.sub.b] is defined by the minimal radius of contacting particles and may be characterized by fraction factor [lambda] (0 [less than or equal to] [lambda] [less than or equal to] 1)

[R.sub.b] = [lambda]min ([R.sub.i], [R.sub.j]) (4)

The cylinder length is denoted by [L.sub.b]. Estimation it of requires broader discussion which will be presented below.

The material of each particle and apparent bond is assumed isotropic and elastic. Elasticity properties of the particles and bond are characterised basically by the elasticity modulus [E.sub.i], [E.sub.j] and [E.sub.b]. The shear moduli and Poisson's ratios are specified in the same manner.

Motion of contacting spheres in time t is defined by normal displacements of spheres centres [u.sub.i](t) and [u.sub.j](t). Parameters [u.sub.i] and [u.sub.j] are assumed to be much smaller than the radii [R.sub.i] and [R.sub.j] of the spheres and radius [R.sub.b] of the bond.

Deformation behaviour of normal particle-particle contact is characterised by the normal inter-particle displacement (particles overlap) [h.sub.ij](t) which depends on time and is expressed as

[h.sub.ij] (t)= [u.sub.i] (t) + [u.sub.j] (t) (5)

Finally, the linear constitutive relationship relating forces and displacements

[F.sub.n,ij] = [K.sub.n,ij][h.sub.ij] (t) (6)

here, [F.sub.n,ij] is the normal force component discussed in Eq. (3). [K.sub.n,ij] is the resultant normal stiffness reflecting various properties of contacting particles and the bond per unit length. It may be obtained analytically or numerically by conducting various tests. Thereby, test loading can be imposed by applying time-driven (Fig. 1, b) or force-driven approaches.

2.2. Normal contact model with parallel bond

Most of developments consider the normal contact model of bonded spherical particles where influence of the interface material is evaluated by applying concept of the parallel bond. This model was proposed by Potyondy and Cundall [18].

Detailed analysis shows that it considers the single point contact of spherical particles when the space between the spheres in the vicinity of contact point is fulfilled by the interface material.

Geometry of the particles and the bond defined by this model is shown in Fig. 2, a. The material bond is presented by the cylinder variable length of which [L.sub.b] is predefined by its intersection with the sphere. Most of the considerations [15-18] are restricted to the limit case of porousless solid, where cylinder is defined by the maximum cylinder radius [R.sub.b] = R covering entire sphere of radius R, while the bond length [L.sub.b] is equal to sum of contacting particles radii [L.sub.b] = 2R.

[FIGURE 2 OMITTED]

The contact model is assumed as direct normal frictionless contact of two spheres while influence of the interface material is considered as parallel spring (Fig. 2 b). Consequently, the resultant normal force [F.sub.n,ij] acting in the bonded contact is comprised of two normal force components

[F.sub.n,ij] = [F.sub.pn,ij] + [F.sub.bn,ij] (7)

n,ij pn,ij bn,ij V /

The first term [F.sub.pn,ij] in Eq. (7) arises from particle-particle contact, while the second reflects the influence of the parallel bond. Separate particle-particle contact is presented in Fig. 3, a. Here, initial positions of particles i and j are marked by dashed line. New positions after deformation are denoted by [[O.sup.'].sub.i] and [[O.sup.'].sub.j], respectively. Due to rotational symmetry about z axis, contact surface between particles is a plane, which is described by the area of the circle in the Cartesian coordinate system x, y, z.

Formally, expression (7) may be simply transformed to the constitutive relationship (6) relating forces and displacements, where the bond properties may be characterised by the resultant stiffness comprising contribution of the particles and the bond

[K.sub.n,ij] = [K.sub.pn,ij] + [K.sub.bn,ij] (8)

here [K.sub.pnij] is the resultant normal stiffness of contacting particles given by expression

[K.sub.pn,ij] = [k.sub.i][k.sub.j]/[k.sub.i] + [k.sub.j] (9)

Individual stiffness [k.sub.i] and [k.sub.j] of individual contacting particles are defined by Potyondy and Cundall [18] as follows

[k.sub.i] = [alpha][E.sub.i][R.sub.i] (10 a)

[k.sub.j] = [alpha][E.sub.j][R.sub.j] (10 b)

In Eqs. (10, a) and (10, b) equations, parameters [E.sub.i], [R.sub.i] and [E.sub.j], [R.sub.j] is the single particles elasticity moduli and radii, respectively. Coefficient [alpha] depends on the amount of the parts of equal division of the sphere. For full 3D spherical particle [alpha] = 4 (Fig. 2, a).

[FIGURE 3 OMITTED]

Second term [F.sub.bn,ij] in Eq. (7) arise from apparent parallel bond and can be denoted as the cement behavior. Illustration of parallel bond is given in Fig. 4, a.

[FIGURE 4 OMITTED]

In the Fig. 4, a denoted parameters are described previously in Fig. 2, a, while new positions of particles i and j described in Fig. 3, a.

The second stiffness term [K.sub.bn,ij] in Eq. (8) reflects behaviour of the parallel bond. It is considered as the stiffness of a simple axially loaded rod, and is given by expression

[K.sub.bn,ij] = [E.sub.b] - [A.sub.b]/[L.sub.b] (11)

here [E.sub.b], [A.sub.b] and [L.sub.b] is the apparent elasticity modulus, the area of the circular cross-section and the length of the bond, respectively.

The main drawback of the previous model (7) relies on the assumption of particle-particle contact. In reality this situation is rather exceptional, while particles interaction is transferred via bond of the finite length.

Modified model of normal contact for spherical particles contacting via interface of the finite length [L.sub.k] = [L.sub.c] is newly elaborated in this paper. Modification was carried out by setting up supplementary sequential spring for bonding material imposing real distance between contacting particles. Modified model is illustrated in Fig. 5.

[FIGURE 5 OMITTED]

The new model presents model of two parallel springs which was already described by the above expressions (7) and (8), where the stiffness of the parallel bond [K.sub.bn,ij] is calculated by Eq. (11). When contacting spheres radii are equal, the length of parallel bond is reduced by extracting the interface length

[L.sub.b] = 2R - [L.sub.c] (12)

The essential difference of the new model is in evaluation of the inter-particle stiffness. Thereby, deformation of interface layer of thickness [L.sub.c] shaded by the darkest gray colour in Fig. 5, a is evaluated by introducing additional sequentially connected interface spring kc. This stiffness is evaluated by the linear model and can be calculated by the following equation

[k.sub.c] = [E.sub.c][A.sub.c]/[L.sub.c] (13)

where [E.sub.c] = [E.sub.b], [A.sub.c] = [A.sub.b] and [L.sub.c] is the interface elasticity module, the cross-section area and the length (distance between contacting particles), respectively. Cross-section area of interface can be calculated by circle area formula.

Finally, the resultant sequential stiffness [K.sub.pn,ij] in Eq. (7) comprises three terms and is obtained by formula

[K.sub.pn,ijc] = [k.sub.i][k.sub.j][k.sub.c]/[k.sub.j][k.sub.c] + [K.sub.i][k.sub.c] + [k.sub.i][k.sub.j] (14)

where [k.sub.i] and [k.sub.j] is the single particles stiffness described by the Eqs. (10, a) and (10, b) expressions respectively.

3. Investigation of the bond properties

Influence of the normal contact interaction between two bonded spherical particles with interface layer was examined by comparing two different analytical contact models. Here contact without interface is considered using Eqs. (7)-(11) and the newly developed model with interface using Eqs. (12)-(14). Two equal spherical particles were assumed for investigation purposes.

[FIGURE 6 OMITTED]

To preserve reality, the model could be useful to imitate concrete. The assumed radius of spherical particles [R.sub.p] = 7.5x[10.sup.-3] m corresponds to the average radius of real granite particles being frequently used as a filler material. Initial gap interface adopted as real matrix of cement stone with varying radius by the ratio [L.sub.b]/[R.sub.p].

Thereby, various parameters of the bond were investigated. The initial distance between the spherical particles is assumed to be fraction of the radius [L.sub.c] = 0.07[R.sub.p]. Variation of the normal contact force against relative bond radius of various bond properties defined by ratios of elasticity moduli [E.sub.p]/[E.sub.b] = 0.2, [E.sub.p]/[E.sub.b] = 0.5, [E.sub.p]/[E.sub.b] = 1.0, [E.sub.p]/[E.sub.b] = 2.0 and [E.sub.p]/[E.sub.b] = 5.0 is considered. Simulation results are presented in terms of normal contact force expressed by Eq. (6). Assuming contact overlap [h.sub.ij] = 0.04[R.sub.p], this force reflects stiffness inter-particle contact.

Comparison of bond models by variation of the normal contact force against relative bond radius of various bond properties is presented in Fig. 6. Here the force is given in the dimensionless form by the ratio F = [F.sub.n]/([E.sub.p][R.sub.p]x[10.sup.-3]), while bond radius is presented by the ratio with respect to the particle radius [R.sub.b]/[R.sub.p]. As the result each part of the graph corresponds to the contact stiffness.

Numerical results show that the models are quite identical for the higher values of elasticity ratio ([E.sub.p]/[E.sub.b] < 1.0). Weakening of the bond above [E.sub.p]/[E.sub.b] > 1.0 leads to considerable differences for small radius of the bond.

4. Finite element simulation of normal interaction

4.1. FE model

Evaluation of the analytical models was performed by conducting of the finite element analysis. Two bonded spherical particles with the interface layer are investigated numerically by using the FEM ANSYS 12.1 Workbench software [21].

The computational domain was identical to those used in analytical calculation as shown in Fig. 5. Geometry of the FE model is defined parametrically and controlled by the basic parameters such as the radius of the particles [R.sub.i] = [R.sub.j] = [R.sub.p], the radius and the length of the cylinder bond [R.sub.b] and [L.sub.b], respectively, and the initial distance between contacting particles [L.sub.c]. Because of the axial symmetry the quarter of the half spheres will be considered in computational model which is shown in Fig. 7.

[FIGURE 7 OMITTED]

Elasticity properties of the spheres and the bond material are described by the elasticity modules [E.sub.i], [E.sub.j], and [E.sub.b], respectively. Interface elasticity properties adopted by the ratio [E.sub.p]/[E.sub.b], which is defined by values 0.2, 0.5, 1.0, 2.0 and 5.0. Poisson's ratios for the spheres and interface for the simplification are constant, [v.sub.i] = [V.sub.j] = [v.sub.b] = v = 0.25.

The symmetric boundary conditions are specified on the vertical [O.sub.i]yz and the horizontal [O.sub.i]xz coordinate planes for both spherical particles and the bond. Actually, their motion was restricted into perpendicular directions. The curved surfaces of the particles and interface cylinder are free.

Loading was imposed by the motion of the central sections y[O.sub.i]x and y[O.sub.j]x of the spherical particles, which in Fig. 7 shaded in gray, and controlled by the displacements. Loading was restricted to maximum displacement value [u.sub.i](t) = [u.sub.j](t) = u(t) = 0.02[R.sub.p].

The volume domain of the particles and the bond are discretised by the tetrahedron 10 nodes elements with six degree of freedom, which are interconnected only by nodes.

Generation of the FE mesh was performed on a basis of the following strategy. Two regions were distinguished in the solution domain. The particles interface zone, which was conditionally limited by the fictive sphere, is covered by the finer mesh. Mesh is characterised by the average element size equal to 0.007[R.sub.p]. Remainder region is covered by the course mesh having average element size equal to 0.07[R.sub.p]. The mesh for the model is relatively small bond radius [R.sub.b]/[R.sub.p] = 0.2 is shown in Fig. 8, a, while refine zone is shown in Fig. 8, b. This model consists of 1223504 nodes and 891960 elements. It should be noted that by increase of the bond radius ([R.sub.b] > 0.5[R.sub.p]) the size of finer elements zone was increased.

Connection type between spherical particles and interface is assigned as bonded. This contact approach is applied to all contact regions (surfaces, solids, lines, faces, edges). In this contact no sliding or separation between faces or edges is allowed.

[FIGURE 8 OMITTED]

Numerical tests with coarser meshes were performed to check influence of the element size. Results confirmed that the above illustrated mesh was quite sufficient to satisfy convergence criteria. The density of the accepted relatively dense mesh with the characteristic element size was expected to be satisfactory for evaluation of exact solution. Normal contact force [F.sub.n,ij] is relatively insensible to mesh refinement.

4.2. Performance of the FE model

Several numerical tests were conducted to study particles interaction behaviour. To limit the number of possible situations some of the model parameters have been fixed for the all samples, i.e. the radius of the particles [R.sub.i] = [R.sub.j] = [R.sub.p] = 7.5x[10.sup.-3] m and the initial distance between the spherical particles is assumed to be fraction of the radius [L.sub.c] = 0.07[R.sub.p].

On the basis of analytical results, see Fig. 6, only the weakest bond was examined. The radius of the cylinder bond is assumed equal to 0.2[R.sub.p]. The bond material properties are defined by ratios of elasticity moduli [E.sub.p]/[E.sub.b] = 0.2.

The geometrically linear and nonlinear (large strain and large displacements) models using the same previously described FE mesh were compared.

Calculation results are presented in Fig. 9. The normal force is presented in the dimensionless form by the ratio F = [F.sub.n]/([E.sub.p][R.sub.p] x [10.sup.-3]). The displacement is presented by the ratio with respect to particle radius [h.sub.ij]/[R.sub.p].

[FIGURE 9 OMITTED]

Obtained results in Fig. 9 seem to be illustrating well-known tendency observed by compering linear and nonlinear contact models. The nonlinear solution exhibits nonlinear geometric stiffening and reflects Hertz solution. In the range of the small overlap differences are insignificant. Consequently, usage of geometrically nonlinear models is redundant. Because loading (maximum displacement value) of particles is such a small that assumption of small strain and small displacement is satisfied. Therefore, further normal contact calculations of FE models were performed by linear principle.

The most of simulation effort was put, however, to examine influence of the bond parameters. The obtained force-displacement relationship for elasticity module ratio [E.sub.p]/[E.sub.b] = 2.0 and various radii ratios between the particles and the bond [R.sub.b]/[R.sub.p] equal to 0.2, 0.6 and 1,0 is shown in Fig. 10.

From graphs shown in Fig. 10 it could be observed, that for fully bonded equal particles, when [R.sub.b]/[R.sub.p] = 1, the obtained 5.5% difference between the FE and analytical models may be considered as insignificantly. Normal forces at the maximum particles displacements by the ratio [R.sub.b]/[R.sub.p] = 0.6 - differ about 11%. For quite small bonds ([R.sub.b]/[R.sub.p] < 0.2) difference between normal forces reach about 40%.

[FIGURE 10 OMITTED]

The character of obtained results may be explained on qualitative analysis of deformation behaviour. Two-dimensional plots of directional deformation fields by the z axis in the central section for various radii ratios [R.sub.b]/[R.sub.p] are shown in Fig. 11. It could be noticed that interaction behaviour highly depends on the bond radius, while the same tendencies observed in previously shown normal force-displacement relationship (Fig. 10) are confirmed by the three-dimensional continuum model. Firstly, transition from Hertz contact (Fig. 11, a) to continuous cylinder (Fig. 11, c) was clearly indicted. Secondly, assumption (14) about uni-directional uniform deformation of the interface layer is held through all of three plots.

[FIGURE 11 OMITTED]

In the calculations of this problem by the FE program it have been encountered some accuracy problems. Therefore, can be said, that not only the mesh size, but also several other factors, such as the size of load increment, definition of initial contact or solution algorithm contribute to simulation results.

5. Evaluation of analytical and FE model results

The finite element results are explored to evaluate approximate analytical models. Regarding very fine mesh and three dimensional models we assume that FE results are very accurate and maybe used for evaluation of approximate analytical expressions (4-14). The FE calculations are restricted to the mostly heterogeneous case, where bond between particles is much lower than particles.

Detailed comparison of analytically and numerically obtained results of contact normal forces varying radii ratio and elasticity moduli ratio of the particles and the bond is given in Fig. 12. Here the dashed lines indicate results obtained by two analytical models, while solid curve indicate FE solution. Additionally the horizontal dotted lines indicate interaction of two spheres without any bond.

[FIGURE 12 OMITTED]

From Fig. 12 can be observed that analytical models exhibit quite different behaviour in the range of the small bond radii. Moreover neglecting of the interface size yields solution bounded by contact model of unbounded particles. Decreasing of the bond radius of the finite size leads to decrease of contact stiffness. This fact obtained by analytical model of the particles-bond with the finite size and is confirmed by the FE calculations.

Comparison calculation results of the newly developed analytical model and FE numerical model gave better coincidence, because newly developed model for normal bonded contact estimates real distance between contacting particles.

Finally, it could be stated, that evaluation of the interface layer is necessary. It could be evaluated by additional sequential spring by expressions (13-14).

7. Conclusions

Basing on the calculation results of two analytical models and the FE simulation the following conclusions were drawn:

1. The numerical simulation of the normal contact interaction between two spherical particles obtained by the FEM using relatively fine mesh is assumed to be exact.

2. It was found that the model of contacting particles with the parallel bond proposed by Potyondy and Cundall [18] gives overestimates contact stiffness, because the bond volume overlaps the spheres.

3. The above model of contacting particles with the parallel bond is suitable for the evaluation of contact with interface, for the relatively strong bonds ([E.sub.p]/[E.sub.b] < 1.0), but for weaker bonds it is suitable for thick bonds with relative radii ([R.sub.b]/[R.sub.p] > 0.6). This model fails however for the case relatively weaker slender interface bonds.

4. For the evaluation of interface, the new analytical model was suggested, where influence of interface of the finite size is evaluated by additional sequential bond. Accuracy of this model was confirmed by the finite element analysis.

5. Further investigations are still required to evaluate not only compression, but also shearing and bending of the bond.

10.5755/j01.mech.18.6.3160

Received August 22, 2011 Accepted November 15, 2012

References

[1.] Cundall, P.A; Strack, O.D.L. 1979. A discrete numerical model for granular assemblies, Geotechnique 29(1): 47-65. http://dx.doi.org/10.1680/geot.1979.29.L47.

[2.] Dziugys, B.; Peters, J. 2001. An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers, Granular Material 3(4): 231-277. http://dx.doi.org/10.1007/PL00010918.

[3.] Kruggel-Emden, H.; Simsek, E.; Rickelt, S.; Wirtz, S.; Scherer, V. 2007. Review and extension of normal force models for the Discrete Element Method. Powder Technology 171(3): 157-173. http://dx.doi.org/10.1016/j.powtec.2006.10.004.

[4.] Tomas, J. 2004. Fundamentals of cohesive powder consolidation and flow, Granular Matter (7): 75-87. http://dx.doi.org/10.1007/s10035-004-0167-9.

[5.] Maknickas, A.; Kacianiauskas, A.; Kacianauskas, R.; Balevicius, R.; Dziugys, A. 2006. DEM software for simulation of granular media, Informatica 17(2): 207-24.

[6.] Johnson, K.L. 1985. Contact Mechanics. Cambridge University Press. Cambridge: 472p.

[7.] Jasevicius, R.; Kacianauskas, R. 2007. Modeling deformable boundary by spherical particle for normal contact, Mechanika 6(68): 5-13.

[8.] Adams, G.G.; Nosonovsky, M. 2000. Contact modeling forces, Trib. Int. 33: 431-442. http://dx.doi.org/10.1016/S0301-679X(00)00063-3.

[9.] Brizmer, V.; Kligerman, Y.; Etsion, I. 2007. The effect of contact conditions and material properties on the elasticity terminus of a spherical contact, International Journal of Solids and Structures 43: 5737-5749.

[10.] Dvorkin, J.; Nur, A.; Yin, H. 1994. Effective properties of cemented granular materials, Mechanics of Materials 18: 351-366. http://dx.doi.org/10.1016/0167-6636(94)90044-2.

[11.] Zhu, H.; Chang, C.S.; Rish, J.W. 1996. Normal and tangential compliance for conforming binder contact I: Elastic binder, International Journal of Solids and Structures 33: 4337-4349. http://dx.doi.org/10.1016/0020-7683(95)00238-3.

[12.] Chang, C.S.; Gao, J. 1997. Rheological modeling of randomly packed granules with visco-elastic binders of Maxwell type, Computers and Geotechnics 21: 41-63. http://dx.doi.org/10.1016/S0266-352X(97)00012-8.

[13.] Cheung, C.Y.; Cocks, A.C. F.; Cebon, D. 1999. Isolated contact model of an idealized asphalt mix, International Journal of Mechanical Sciences 41: 767-792. http://dx.doi.org/10.1016/S0020-7403(98)00046-0.

[14.] Zhu, H.; Nodes, J.E. 2000. Contact based analysis of asphalt pavement with the effect of aggregate angularity, Mechanics of Materials 32: 193-202. http://dx.doi.org/10.1016/S0167-6636(99)00054-X.

[15.] Jiang, M.; Leroueil, S.; Konrad, J.M. 2005. Yielding of microstructured geomaterial by distinct element method analysis, Journal of Engineering Mechanics 131(11): 1209-1213. http://dx.doi.org/10.1061/(ASCE)07339399(2005)131:11(1209).

[16.] Jiang, M.J.; Yu, H.S.; Harris, D. 2005. A novel discrete model for granular material incorporating rolling resistance, Computers and Geotechnics 32(5): 340-357. http://dx.doi.org/10.1016/jxompgeo.2005.05.00L

[17.] Jiang, M.J.; Yu, H.S.; Harris, D. 2006. Bond rolling resistance and its effect on yielding of bonded granulates by DEM analyses, International Journal For Numerical And Analytical Methods In Geomechanics 30(8): 723-761. http://dx.doi.org/10.1002/nag.498.

[18.] Potyondy, D.O.; Cundall, P.A. 2004. A bonded-particle model for rock, International Journal of Rock Mechanics and Mining Sciences 41(8): 1329-1364. http://dx.doi.org/10.1016/jijrmms.2004.09.01L

[19.] http://www.dem-solutions.com/.

[20.] Itasca Consulting Group Inc. PFC2D/3D (Particle Flow Code in 2/3 Dimensions), Version 2.0. Minneapolis, MN: ICG; 1999.

[21.] http://www.ansys.com/.

S. Pilkavicius *, R. Kacianauskas **, A. Norkus ***

* Vilnius Gediminas Technical University, Sauletekio 11, 10223 Vilnius, Lithuania, E-mail: saulius.pilkavicius@vgtu.lt

** Vilnius Gediminas Technical University, Sauletekio 11, 10223 Vilnius, Lithuania, E-mail: rkac@vgtu.lt

*** Vilnius Gediminas Technical University, Sauletekio 11, 10223 Vilnius, Lithuania, E-mail: Arnoldas.Norkus@st.vgtu.lt