J. Indones. Math. Soc. (MIHMI) Vol. No. , pp. 173Ae189. A MODEL FOR STATIC AND DYNAMIC PHENOMENA IN DEPOSITION PROCESSES Hasan and J. van Opheusden Abstract. We have developed a model for friction in a dry granular material, that allows multiple reversible transitions between stick and slip contacts. The ultimate purpose of the model is to simulate deposition processes. During these processes a pile structure is formed that is repeatedly destabilised by addition of new material. The present model does not yet include rotation. The model is tested for a few simple systems. Finally we perform more extensive granular dynamics simulations using the model for 2D and 3D In each system, a steady flow of material is dropped onto a rough surface. INTRODUCTION When dry granular materials are dropped onto a surface or into a container a pile structure is formed. The characteristics of the pile, its structure and stability, are largely determined by the mechanical properties of the material, and the details of the deposition process. In the early stages grains form a small pile, which would be stable if left undisturbed. The further pile production depends on the reaction of the existing structure to disturbances caused by newly dropped grains. Grains previously at rest can start sliding before coming to rest again. This implies there can be many transitions between static . or dynamic . contact during the whole process. Several experiments on real systems that investigate pile structure during deposition have been reported in the literature. Avalanches are present during the pile production, and at the end a stable pile is formed which is almost triangular . in shape with tails near the base . , . and a flat curve at the top . a dropped grain having large impact energy hits a stable pile, a crater is formed . Received 8 November 2006, revised 4 January 2007, accepted 15 November 2007. 2000 Mathematics Subject Classification: 74E20 Key words and Phrases: deposition processes, granular dynamics Hasan and J. van Opheusden In order to establish further insight into the relation between quantities at microscopic level and overall system behaviour, numerical simulations using approximate models are very useful. Various simulation methods have been developed, including Cellular Automata . , . , and Discrete Element Methods (Granular Dynamic. We will use the latter method. The Granular Dynamics method traces the trajectory of each grain, using a numerical integration of NewtonAos equation of motion. The force on each grain is calculated at every time step. This force is then used to determine the motion of the grain, giving a dynamic simulation. The most important aspect of our model is how we describe the friction force between the particles. During the deposition process there may be many transitions between slip and stick modes for all contacts involved. The Coulomb law of friction states that once the tangential force during the static mode exceeds a certain threshold, the dynamic mode engages, during which the friction force is just proportional to the normal force. Unfortunately this simple law can not be implemented directly in the dynamic simulation due to the unknown force during the static mode. Some models have been used to avoid this problem in investigating the pile formation. One of the models applies a friction force proportional to the relative tangential velocity at small velocities . In this case at rest there is no friction, and the system gradually yields to any external force. Hence, this model is not suited to fully describe deposition, as in the end any pile will flatten out due to gravity. Other models introduce a virtual spring at contact . , . This gives a more realistic description of the static regime. However, during the dynamic mode this model may generate an extra force due to this virtual spring when the total tangential displacement and the relative tangential velocity are in opposite directions . In our model we introduce a virtual spring whenever the relative tangential velocity at a contact point decreases below a certain threshold. This gives a transition from a slip to a stick mode. When the virtual spring is stretched beyond another threshold, similarly to the Coulomb law, a transition occurs to the slip mode, and the virtual spring is removed. Thus both transitions from slip to stick mode and vice versa are possible within our model. During the deposition process each contact may undergo many of these reversible transitions. We start with a detailed example explaining why the standard Coulomb law of friction leads to possible difficulties in a fully dynamical multi-particle model. Next we specify our model with the critically damped virtual springs that break upon overstretching to have a quasi-static regime within a dynamical model. The model is first applied to a single contact point of a tethered sphere on a conveyor For this simple system the equation of motion can be solved analytically, allowing direct comparison with simulation results. Also the Coulomb law can be used, as there is only a single contact involved. To investigate whether our model is suited to describe stable pile formation we consider a three dimensional system with five identical soft spheres. Four spheres rest on a rough horizontal surface, the fifth falls onto the array due to gravity. Finally we apply our model to two and three dimensional pile formations of soft discs and spheres respectively. Upon collision there is the tangential force as specified above, in the normal direction there is a Model for static and dynamic phenomena compression force and a dissipative term. Rotation is excluded. Spherical particles are employed to model the shape of the grains. Thus our model will not produce an adequate description of actual deposition processes. Spherical particles will roll, and aspherical particles will have different packing structures from the ones observed in this model study. The main purpose of this study is to test the multiple reversible transitions between stick and slip contacts. Rotation could in principle be added to the model straightforwardly, albeit with considerable effort, mainly with respect to coding. Without rotation also the aspect of rolling friction can be MAIN RESULTS The Coulomb law of friction The law of Coulomb is commonly applied to describe the friction force between colliding dry materials. It relates the normal and the tangential forces as given by the following equation [OeAAs . n |. AAs . n |], vt = 0 . tatic frictio. OevCt AAd . n |, vt 6= 0 . ynamic frictio. where vt and vCt are the relative tangential velocity and a corresponding unit vector, fn is the normal force and ft is the tangential . force, while AAs and AAd represent the coefficients of static and dynamic friction. It can be seen from equation . that the Coulomb law is simple in the sense that the coefficients of friction are constant and the velocity dependence is only at vt = 0. Due to its simplicity, this law has been widely implemented to model the tangential force in simulating the dynamic behaviour of frictional granular systems as discussed in the previous section. On the other hand, it can also be seen that there are some difficulties in dynamic simulations, caused by the Coulomb One is that the tangential force does not change gradually during a staticdynamic transition. Moreover during the static regime, the tangential force can not be determined directly from the state variables, i. , all is known is that the magnitude of the tangential force is inside the so called friction cone . s indicated by the interval in the equation abov. All the other forces acting on the two colliding bodies determine its actual value. This value assures that vt remains Depending on the number of spheres and the number of contacts there may or may not be unique solution for the tangential forces to maintain a static state. If every sphere has only a single contact, either with another sphere or with the boundary, the Coulomb law can be directly implemented since the friction force at the contact point can always be determined. The only realisation of such a system in our study is the one where each individual sphere is just resting on the bottom surface. For pile formation there must be at least two contacts per sphere, to have a single vertical array of grains, but for a realistic pile three or Hasan and J. van Opheusden more contacts per sphere are needed. Without rotation the number of equations to determine whether a pile of N particles is stable . ll accelerations zer. is 2N in two and 3N in three dimensions. To describe the size of the tangential forces M unknowns must be solved for, and another M for their directions . n 3D). It implies that the system of equations becomes indeterminate in 2D when there are more than 2 contacts per grain, against 1. 5 contacts in 3D. As stated, two contacts are the minimum needed to have a realistic pile, and in general the number of contacts in a 3D pile will be higher than that in a 2D one. Indeterminacy of the tangential forces is related to existing internal stresses within the pile. In reality these stresses will build up during the pile formation, so history is important. A full model to describe these internal stresses hence must include the history of the pile. Another complication is that it is not clear, if the static state cannot be maintained . hat is when the system of equations has no solutio. which contact starts slipping. Again historical information is needed to resolve that problem in practice. Before we try to deal with the complexity of multi-particle contacts, let us first consider a system in which a sphere is wedged into two fixed spheres, as shown in Fig. The free sphere 2 has two contacts with the fixed spheres 1 and 3. Here we show that even for such a simple system, equation . can not fully describe the dynamics. The only moving sphere is sphere 2, the others are considered fixed, and we have three dimensional system. There are two contacts, so there are four unknowns and just three equations for zero acceleration of sphere 2 in the three Apparently one degree of freedom is left, as we will show in more detail With rotation another three equations would be given by zero rotational acceleration of sphere 2, and we would not have the problem of indeterminacy, but rather too many constraints to be met. Figure 1: The collision geometry . , and dependency of collision modes on c . The spheres within our model are considered to be elastic and compressible, which is represented by an overlap in this figure. The compression of the spheres is associated with the . compression force at the contact. At each contact Model for static and dynamic phenomena point, the related normal direction and tangential plane can always be identified. At each contact the colliding spheres exert a force onto each other with opposite The force consists of a normal compression force, acting in the normal direction, e. nC21 , and a tangential friction force, acting in the tangential plane spanned by tC1 and tC2 . Suppose the line of intersection between the two tangent planes is in the direction of uC2 , one of the unit vectors forming the tangent plane where sphere 2 collides with sphere 3, then we can establish a co-ordinate system consisting of the unit vectors tC1 , uC1 , and uC2 ( tC2 and uC2 are collinea. Let ft21 and ft23 be the friction forces at the contact points of colliding sphere 2 with spheres 1 and 3 respectively. As derived in appendix A1, there is one free parameter, c. This parameter c can have any value in the range cl O c O cr , to give a static In figure 1. , ft21 . and ft23 . are the friction forces, and ft21 ft23 denote the maximum static friction forces at the related contact points. which a transition to slipping occurs within the Coulomb model. These maximum forces determine the critical values of the free parameter. If we choose the value of c in the interval ccrit then the contact between spheres 1 and 2 is in the dynamic mode, cd21 , while the contact between spheres 2 and 3 is in the static mode, cs23 . Assigning a value of c in the other intervals yields contact modes as shown in figure 1. , where the superscript AusAy and AudAy stand for AustaticAy and AudynamicAy. Between the critical values there is a regime where both contacts can remain static. When the static contact mode is being maintained, the tangential force is evaluated using equations . On the other hand, when the colliding spheres start sliding the tangential force is given by the dynamic friction of equation . Depending on the value of c, figure 1. gives a tool to classify the contact modes and hence the tangential forces. Unfortunately, this value can not be determined directly from the state variables. Thus, it is clearly shown that the Coulomb law can induces force indeterminacy in dynamic simulations where a sphere has two or more simultaneous collisions in three dimensions, as well as indeterminacy as to where a stick-slip transition will occur if stability is compromised. 2 Force models 1 Tangential force To avoid the difficulties of Coulomb friction, fully dynamical models, where the static mode is approximated with little or no motion. Elperin and Golshtein . and Matuttis et al . introduce a tangential virtual spring once two spheres collide at time t0 , as is pioneered by Cundall and Strack . The model is represented by OeCt kt . | if kt . | O AA. n |, . OeCt AA. n | where kt denotes the tangential spring constant, t = t0 vt (E )dE is the total tangential displacement built up during the time duration t Oe t0 , and the cap means an associated unit vector, indicating the direction only. Hasan and J. van Opheusden The first term within the bracket in equation . represents the spring force during the AustaticAy mode. The relative velocity may not be zero, but if the spring constant is sufficiently large, the tangential displacement will be small. When the force exceeds the threshold, given by the constraint, the spheres start sliding, and the friction force is given by the second term opposing the tangential displacement. Thus, there is a transition criterion from static to dynamic mode. This model gives a good agreement with experimental data for acetate spheres . However, further investigations show that this model may lead to unphysical behaviour in the tangential force . , . since in the model the friction force opposes the tangential displacement rather than the relative tangential velocity. Another shortcoming of this model is that it does not allow for a transition back from the dynamic to the static mode. Finally in the static mode the system may oscillate forever. Figure 2: The transition of motion . Elperin and Golshtein . and Brendel and Dippel . have made various suggestions to remedy the shortcomings mentioned. We have developed a model that includes a dynamic and a . static mode, and transitions between these, to describe the deposition process. In the static regime we use a . harmonic spring, as in the Cundall and Strack model, but with a damping term included. Furthermore, rather than engaging static mode at contact, an AuA - criterionAy as shown in figure 2 is introduced: once the magnitude of relative tangential velocity is less than A, the static regime engages, otherwise, the dynamic regime applies. When the spring force exceeds the Coulomb threshold, the spring is removed and the dynamic mode engages. During the quasi-static regime, the friction force is given by the damped spring, while beyond this regime the friction force is just proportional to the normal force. In addition, the model applies two different coefficients of friction. AAs and AAd . The tangential force is then modelled by Oe. t t t vt ), . t | O A . tatic mod. OevCt AAd . n |, . t | > A . ynamic mod. where t is the tangential damping coefficient, and the remaining notations are the same as above. In our case t0 is not the time of contact, as applied in . , but is Model for static and dynamic phenomena the time when the quasi-static regime engages. Oscillations in the static regime are avoided by applying critical damping, i. Oo t = 2 mkt . At the onset of the dynamic regime the velocity is small. To avoid the system engaging static mode again, we make the additional requirement that the tangential force opposes the velocity when applying the A - criterium. As long as the velocity is increasing the dynamic mode will pertain. Hence the discontinuity at the transition from dynamic to static mode is replaced by a . the transition from static to dynamic the virtual spring is completely removed, all information about the direction or size on breaking is lost. The simplicity of this model is that, unless the A - criterion is matched, it only applies dynamic friction, which is velocity independent. A similar velocity independent model is also applied in . , . , but the virtual springs are implemented Experimental results, using spherical glass particles . showed complicated behaviour of the stick - slip transition, with a continuously changing force, instead of a discrete step. It would be difficult to include these effects in the model. 1 Normal force The second component of the contact force in our model is the compression force, perpendicular to the friction plane. In our model, it consists of an elastic and a dissipation force, modelled by a Hookean spring and damper fn = Oe. n n n vn ), . where n and vn are the displacement . pheres may overla. and the relative velocity in the normal direction, and kn and n are the normal elastic and damping The displacement is determined in practice by treating the grains as perfect spheres with a constant radius, and the displacement is simply the overlap between the . The tangential plane hence is always simply perpendicular to the vector connecting the centres of the colliding spheres. As the normal force, fn , is not constant, any quantities determined by the normal force, , dynamic friction and maximum static friction, may vary over time, inducing transitions between static and dynamic modes. The normal force . is simple in the sense that the elastic force is a linear function of the displacement, and the dissipation force is a linear function of the relative normal velocity. The simplicity leads to a non-zero force at the onset of contact due to the damping term . , . , which may not be appropriate. the other hand, this simple model results in a coefficient of restitution. O, which is velocity independent. By applying the condition suggested by Luding . , i. c ) = 0, where tc is the collision time. O can be related to the damping coefficient A1 n n O = eOe O (AOe2 arctan( O )) , with O = 4mkt Oe n2 2 . Here we used a collision with zero tangential velocity. During actual contacts energy is not only dissipated in the normal direction, but also in the tangential direction through friction. The advantage of this formulation is that it directly Hasan and J. van Opheusden Figure 3: The conveyor belt relates a model parameter, the normal damping coefficient, to an experimental system parameter, the coefficient of restitution. 3 A single sphere on the conveyor belt To test our model with transitions between slip and stick motion we study a simple system of a single incompressible sphere on a conveyor belt. The sphere is connected to a fixed support by a spring having an elasticity constant kr , as shown in figure 3. The sphere is initially at rest with respect to the conveyor belt, and the spring is in the rest position. As the conveyor belt is driven with a constant velocity ve , the sphere moves with the same velocity, vs = ve , causing the spring to stretch, generating the spring force, fsp = Oekx. As long as the magnitude of spring force does not exceed the maximum static friction force, fsmax , the sphere remains stationary with respect to the belt. We say the system is in the static mode. When the spring force exceeds the maximal static friction the sphere starts sliding along the belt, and the dynamic friction, fd , is applied. This we call the dynamic mode. When the relative velocity of the sphere drops below the critical value, the sphere re-sticks to the belt, and we have static mode again. The equation of motion of the sphere is then given by xN = ve , . sp | O fsmax . tatic mod. , . mxO = fsp Oe fd , . sp | > fsmax . ynamic mod. Equations . shows that during the static mode the sphere is displaced until x = fsmax /k. Furthermore, the dynamic mode always starts when x = fsmax /k and xN = ve , and the system returns to static mode when xN = ve or equivalently x = . AAd Oe AAs )mg. Hence, the system imposes the criteria to switch from static to dynamic modes and vice versa. Since the relative velocity, vr = vs Oe ve , is always negative, the second equation of . can be rewritten in the form xN = y yN = Oe m . Oe . kd | ), with the initial conditions x. d0 ) = fsmax /k, and y. d0 ) = ve , where td0 is the time when the sphere starts to slide. The system equation . has an equilibrium point Model for static and dynamic phenomena d |/k, . , and its solution is given by . d | 2 ) y 2 = A2 where A2 = . (OIAA)2 g 2 ve2 , and OIAA = AAs Oe AAd . Substituting the first equation of . into equation . gives x1 = AAs mg/k and x2 = . AAd Oe AAs )mg/k. Both x1 and x2 denote the positions where the transitions The first value is related to the transition from stick to slip, the second one is from slip to stick. The sphere undergoes a periodic stick-slip motion. In the phase plane, it is represented by a straight line during stick mode, and by a part of ellipse during the dynamic mode. Note the equation for x2 shows that if the roughness of the materials is such that AAd O 12 AAs then x2 O 0, i. the spring is compressed prior to the transition. We want a numerical procedure to solve equation . to replicate this, especially we want to avoid oscillatory behaviour near the transition from slip to As a test, we apply our model of tangential friction . that introduces a virtual spring once the quasi-static regime is valid. In the first leg of the curve, the stick regime, we apply the Coulomb law. When the spring force exceeds the maximal static friction force, we use our model with the onset of the dynamic The used parameters were a mass m = 0. 05 kg, the gravitational acceleration g = 9. 81 msOe2 , stiffness of the real spring kr = 105 kgsOe2 , stiffness of the virtual spring kv = 106 kgsOe2 , coefficients of friction AAs = 0. 6 and AAd = 0. 3, stick velocity A = 10Oe4 msOe1 , and belt velocities ve = 5 y 10Oe3 msOe1 . Oe3 Oe3 ve = 0. analytical solution numerical solution Oe2 equilibrium point analytical solution kv = 5x106 Oe4 kv = 10 kv = 10 Oe6 Oe3 Oe2 Oe1 . Oe6 Oe3 . Figure 4: Phase plane of sphereAos motion . and vel. at stick-slip transition . The analytical and numerical solutions of the sphereAos motion in the phase plane are shown in figure 4. At the first loop, both solutions are the same, the Coulomb law and the dynamical model copy the analytical solution within rounding Hasan and J. van Opheusden of error of the algorithm. The deviation occurs when the numerical procedure enters the quasi-static regime. When the dynamic friction disappears, the virtual spring and damping force build up the quasi-static regime. Initially these are less than the dynamic friction, so the velocity drops below the analytical curve. As can be seen from figure 4. , the sphere actually never catches up with the belt. The reason is that in the static equilibrium the real spring force and the virtual spring force = Oekr x = kr ve t, so to match this The real spring force increases as fsp the virtual spring force must be fsp = kv = kv vv t, with vv = ve kr /kv a remaining velocity difference between the belt and the sphere. The more rigid the virtual spring, the smaller the deviation from the analytical solution. 3 A small pile z5 = 0. z = 0. diagonal distance z5 = 0. z5 = 0. O OO 0. s = 0. O OO 0. O OO 0. s = 0. diagonal distance diagonal distance s = 0. s = 0. O OO 0. Figure 5: Five spheres construction . and behaviour of diagonal distance for varying initial height . , init. mutual dist. of spheres . , and coeff. of rest. Figure 5. shows a system of five identical spheres with the following initial four non-colliding spheres are placed at rest on a rough surface in a Model for static and dynamic phenomena square array of side length s, and a fifth sphere is dropped from a vertical position exactly above the centre of the square, hitting the spheres on the floor. We use our model to investigate whether a stable pile structure results, with the upper sphere resting on the slightly extended array of four. In appendix A2 we have investigated under what conditions this system allows a stable static pile structure. Figure 8 shows the critical static friction coefficient for a system consisting of 4, 5 and 6 spheres, as a function of the distance between the lower spheres. If the actual static friction coefficient is less than critical, the pile is unstable. Here we assumed the upper sphere is placed gently on the lower Dropping it from a distance has a destabilising effect. We apply our quasi-static model once the A - criterion is satisfied. Throughout these simulations the used parameters values were diameter E = 0. 05 m, mass m = 0. 05 kg, elasticity constants in the normal and tangential directions kn = kt = 105 kgsOe2 , friction coefficients AAs = 0. 6 and AAd = 0. 3, stick velocity A = 10Oe3 msOe1 , and time step is chosen sufficiently smaller than the smallest time scale occuring within our model system. We have varied the initial positions of the spheres . pper and lower sphere. and the normal damping constant . r the coefficient of restitution via equation . In the simulations we only varied a single parameter, as a reference we used a normal damping constant n = 0. 5ncrit or the coefficient of restitution O OO 0. 3, where ncrit is the normal critical damping. The diagonal distances of the lower spheres are shown in figures 5. Ae. As shown in figures 5. Ae. , when the upper sphere hits the lower spheres, they move apart. Depending on the applied parameters, the final diagonal distances may be less than twice the sphere diameter, in which case a stable pile is observed. If they exceed this value all spheres are on the surface. Figure 5. describes the effect of varying the initial vertical position of the upper sphere, z5 . In this case the initial mutual distance of the lower spheres is 0. 058 m. It can be seen that the higher z5 , the flatter the pile structure. At z5 = 0. 3 m there is no stable pile formed. For this value of z5 a pile is formed if the initial gap is smaller, e. s = 0. 054 m, as shown in figure 5. nitial height of the upper sphere = 0. Figure 5. nitial height of upper sphere = 0. 1 m, and mutual distance of lower spheres = 058 . shows the effect of the coefficient of restitution O. In figure 5. , it can be seen that the system with s = 0. 058 m, z5 = 0. 1 m and O OO 0. 3 forms a stable pile. However, the same system with O OO 0. 7 produces an unstable pile. As can be seen from figure 5. , the upper sphere bounces back and hits the lower spheres again until they give way. In order for a larger pile of many particles to be stable, it is not completely necessary that the small system is stable. Once a somewhat dense layer of particles on the surface has formed, a particle on top of the layer must displace more than just the particle it is touching, and the friction forces involved may allow the formation of a larger pile. 3 Hourglass simulations The previous two examples were simple test cases. The conveyor belt simulation indicates for what values of the calculational parameters a stable algorithm Hasan and J. van Opheusden is obtained, the small pile simulation tells us what system parameters we must choose to obtain a stable pile. This first actual application of the model involves a system where particles fall one by one, from a certain vertical position, and with a small random horizontal velocity, onto a surface which has the same material properties as the objects. We have named these calculations hourglass simulations, as it mimics the free falling grains from the opening in an hourglass. The small horizontal velocity is needed to introduce some randomness in the system, and to avoid having all grains stack exactly on top of each other. We simulate two and three dimensional . D and 3D) systems consisting of 600 identical discs and 1200 spheres respectively. The parameter values were chosen similar to the previous problem, i. E = 0. 05 m, m = 0. 05 kg, kn = kt = 105 kgsOe2 . AAs = 0. AAd = 0. 3, n = 0. 5ncrit , t = tcrit , and A = 10Oe3 msOe1 . Two time scales can be derived from the force equations 3 and 4. The oscillations of a particle of mass p m in the normal direction with a force described by 4 occur on a time scale of m/kn , which is 7 y 10Oe4 Oo s for our values, provided the normal friction constant is small enough . O mkn ). For larger values of the dissipation in the normal direction the motion is critically or overdamped, leading to longer, rather than smaller time Oo As the dissipation in the tangential direction critically damped, t = 2 mkt , the time scale for the tangential force is m/kt , 5 y 10Oe4 s. The time scale in the transition from dynamic to static mode is related to the critical velocity at which the transition occurs. Combination with the dynamic friction force yields a time scale of mA/(AAd . n |), which may be interpreted as the time over which the velocity drops to zero given the force. The problem with this time scale is that it depends on the current value of the force. For the previous applications this time scale is not important, as the normal force is at about the order of the weight of a single grain. For real deposition with a large number of grains the normal force can be large. For a friction force of 10 kgmsOe2 we obtain a time scale of 5 y 10Oe5 s. We set the time step to OIt = 10Oe5 s. Figure 6: Stable piles of 600 discs . , and 1200 spheres . The obtained stable piles are displayed in figures 6. The model generates stable piles for both systems. For the 2D and 3D systems, the different colours of discs and spheres represent, respectively, the first, second and last 200 discs and 400 spheres dropped. It appears that the pile has almost an onion-like layering, as has been observed in experiments . and numerical simulations . Model for static and dynamic phenomena The pile almost forms a triangle with an irregular surface. Furthermore, it has an almost regular close packing structure with dislocation lines. In the dynamical study it was found that along these dislocation lines landslides occur, during which a large part of the pile shifts horizontally. 3D piles must have a much larger number of particles before a substantial pile, suited for analysis, is formed. To shorten the simulation time, grains are dropped in small batches of five. On the other aspects the model and procedure are analogous to the 2D case. Figure 6. shows the front view of the half part of the 3D pile cut by plane y = 0. The pile formed has a somewhat conical shape with spheres dropped earlier covered by the ones dropped Unlike the 2D pile, the structure of 3D pile, both inside and at the surface, is rather irregular. CONCLUDING REMARKS It has been shown that via a virtual spring with damping our friction model is able to describe the phenomena during the contact. from stick to slip, during slipping, and from slip to stick. The model solves the inherent problem of the Coulomb friction law, and has not shown any systematic error. As is shown in the phase plane of the conveyor belt problem, the model produces the same behaviour as the analytical solution, with slight deviation. Furthermore, the hourglass simulations show the model capable of generating stable piles, convincing it adequately captures this aspect of the deposition problem for which it was developed. In its present form the model may be used to study deposition of model materials or pile formations for which rolling is not important, if such materials do exist. REFERENCES