Indeterminacy and the onset of motion in a simple granular packing

Indeterminacy and the onset of motion in a simple granular packing Sean McNamara, Ramón García-Rojo, and Hans Herrmann Institut für Computerphysik, Universität Stuttgart, 70569 Stuttgart, Germany 共Received 25 August 2004; revised manuscript received 19 May 2005; published 19 August 2005兲 We examine the relatively simple problem of a disk placed in a symmetric V-shaped channel, and subjected to gravity and a torque. We obtain analytic predictions of the contact forces and disk motion using two different models. In the first model, the disk is assumed to be perfectly rigid, leading to force indeterminacy. In the second model, the disk is assumed to interact with the walls of the channel via linear springs, leading to a unique solution for the contact forces. The results of these two models are compared. It is shown that there are two possible ways motion can occur—through the appearance of a null eigenvector of the stiffness matrix, or through an instability. When motion occurs through an instability, the first model cannot predict when the disk will rotate; it is necessary to know the undetermined forces in order to predict the motion of the disk. It is also shown how indeterminacy in the first model is linked to memory in the second. The analytical results are also compared with numerical simulations using two different methods, each related to one of the models. DOI: 10.1103/PhysRevE.72.021304 PACS number共s兲: 45.70.⫺n I. INTRODUCTION A. Motivation The quasistatic behavior of granular materials has been studied for many years by engineers seeking to provide stable foundations for buildings and roads. They typically use continuum equations, with appropriate constitutive relations, to predict the deformation of soils when a load is applied to them. Recently, there has been a great deal of work to relate these continuum equations to the micromechanics— what is happening at the level of individual grains or contacts. One hopes to obtain a deeper physical understanding of “granular solids,” as well as to suggest new continuum approaches. One powerful tool employed in this quest is numerical simulation, which gives access both to the macroscopic, continuumlike behavior of the granular material, and also to very detailed microscopic information. However, such simulations always involve approximations and various mathematical subtleties. One of the most intriguing issues raised by numerical simulations is the indeterminacy of packings of rough, perfectly rigid particles. Since the Young’s modulus of rock, glass, and metal is much greater than the stresses applied to granular materials, it is tempting to idealize the grains as infinitely hard bodies. When this is done, the problem of determining contact forces in mechanical equilibrium is indeterminate: there are usually many possible solutions 关1–5兴. Several questions arise: How are the different solutions related to one another? Can indeterminacy cause wrong predictions of particle motion to be given? In a model without indeterminacy, are all possible solutions accessible? What determines which one chosen? The goal of this paper is to study the mathematical properties of two different models of granular material. One model assumes perfectly rigid particles and exhibits indeterminacy, and the other assumes the particles are deformable and produces a unique solution. Both models can be solved analytically when applied to a very simple problem. The results of the two models are then compared, giving precise answers to the questions posed in the proceeding paragraph. Furthermore, our calculations are done using formalisms then can be generalized to many particle packings. We anticipate therefore that our findings in this paper are suggestive of what happens in larger, more complicated systems. B. Synopsis The granular “packing” shown in Fig. 1 is perhaps the simplest system that exhibits some effects that we wish to study. The disk is subjected to a gravitational acceleration g and a torque ␶ ⬎0. The contacts between the disk and the wall are assumed to be noncohesive and frictional, with Coulomb friction ratio ␮ . We investigate two questions: How do the contact forces at contacts ␣ and ␤ change in response to the imposed forces? At what value of the torque does the particle begin to rotate? Similar problems have been investigated in at least two other cases. Ref. 关5兴investigated a rod wedged between two converging walls. Indeed, this work can be considered as a deeper and more precise investigation of the ideas presented in that paper. The system we study is also very close to the disk placed into an asymmetric groove considered in Ref. FIG. 1. A disk supported by two walls through contacts ␣ and ␤ . The angle ␾ suffices to characterize the geometry. PHYSICAL REVIEW E 72, 021304 共2005兲 1539-3755/2005/72共2兲/021304共11兲/$23.00 ©2005 The American Physical Society021304-1 关6兴; however, the perspective taken in this paper is different. Instead of investigating a simple physical problem as realistically as possible, as in Ref. 关6兴, this paper is concerned with the mathematical properties of two different models in the hope of gaining insight into the behavior of many particle systems. These models involve idealizations 共for example, a linear force law is used at contacts, the disk is perfectly round兲whose effects are not analyzed here. We believe it is useful first to understand these simple models before turning to the complexities introduced by lifting these idealizations. We first examine some general considerations: the equations of static equilibrium, and the status of the contacts. Then we attempt to deduce the contact forces from these considerations. We show that for certain geometries and torques, one cannot decide whether the disk rotates or not. We then present a series of contact dynamics 共CD兲关7兴and molecular dynamics 共MD兲关8兴simulations, and point out the difference between them. Then we formulate a second method of calculating the contact forces, assuming they arise from small deformations of the particles. This second method explains all the features of the MD simulations, and sheds light on the questions mentioned above. II. RIGID PARTICLES A. Static equilibrium The equations of static equilibrium for the disk in Fig. 1 are R ␣ sin ␾ +T ␣ cos ␾ −R ␤ sin ␾ +T ␤ cos ␾ =0, R ␣ cos ␾ −T ␣ sin ␾ +R ␤ cos ␾ +T ␤ sin ␾ =mg, rT ␣ +rT ␤ =− ␶ .共1兲 Here, R ␣ indicates the normal force at contact ␣ , with R ␣ ⬎0 for repulsion. T ␣ is the tangential force, with T ␣ ⬎0 when this force exerts a positive torque on the disk. R ␤ and T ␤ are defined in the same way. These equations can be written in matrix form: cF +fext =0.共2兲 We call cthe contact matrix: c= 冢 sin ␾ cos ␾ − sin ␾ cos ␾ cos ␾ − sin ␾ cos ␾ sin ␾ 0r0r 冣 .共3兲 Note that the dimensions of crequire that it have at least one null eigenvalue. The contact forces Fand the external forces fext are F= 冢 R ␣ T ␣ R ␤ T ␤ 冣 ,fext = 冢 0 −mg ␶ 冣 .共4兲 It is convenient to use the following orthogonal basis of R4: Fx=1 2 冢 sin ␾ cos ␾ − sin ␾ cos ␾ 冣 ,Fy=1 2 冢 cos ␾ − sin ␾ cos ␾ sin ␾ 冣 , F ␪ =1 2 冢 − cos ␾ sin ␾ cos ␾ sin ␾ 冣 ,F0=1 2 冢 sin ␾ cos ␾ sin ␾ − cos ␾ 冣 .共5兲 If Fis written in this basis: F=axFx+ayFy+a ␪ F ␪ +a0F0,共6兲 then one can verify that cF = 冢 ax ay axrcos ␾ +a ␪ rsin ␾ 冣 =−fext.共7兲 This equation enables one to easily construct solutions to Eq. 共2兲, for one has immediately ax=0, ay=mg, and a ␪ = − ␶ /共rsin ␾ 兲. However, a0cannot be determined in this way because F0is a null eigenvector of c. This is an expression of the indeterminacy in the problem. B. Contact status Of course, one does not have absolute freedom to choose a0, because the following conditions must be satisfied at the contacts: R艌0, ␮ R艌兩T兩,共8兲 where the constant ␮ is the Coulomb friction ratio. The first inequality excludes cohesive forces, and the second requires that the tangential forces not exceed a certain threshold. At this point, it is useful to introduce the concept of contact status. If we have equality in the second condition in Eq. 共8兲, the contact is said to be sliding, otherwise it is nonsliding. The tangential relative motion vtmust be consistent with the contact status. If the contact is nonsliding, no tangential motion is allowed, but if the contact is sliding, then the tangential force must oppose the motion. Note that vt=0 is allowed under all circumstances, so that contacts in static equilibrium can be sliding. To see how this limits the possible values of a0, one can first use the definition of Fin Eq. 共4兲and the basis in Eq. 共5兲 to compute the forces in terms of ␮ , ␾ , and the coefficients in the expansion of Eq. 共6兲. Then one can use Eq. 共7兲to eliminate all the coefficients except for a0. Finally, one forms the inequalities ␮ R ␣ 艌−T ␣ , ␮ R ␣ 艌T ␣ , ␮ R ␤ 艌−T ␤ , and ␮ R ␤ 艌T ␤ . These yield MCNAMARA, GARCÍA-ROJO, AND HERRMANN PHYSICAL REVIEW E 72, 021304 共2005兲 021304-2 a0共 ␮ tan ␾ +1兲+ 冉 mg + ␶ rsin ␾ 冊 共 ␮ − tan ␾ 兲艌0, 共9a兲 a0共 ␮ tan ␾ −1兲+ 冉 mg + ␶ rsin ␾ 冊 共 ␮ + tan ␾ 兲艌0, 共9b兲 a0共 ␮ tan ␾ −1兲+ 冉 mg − ␶ rsin ␾ 冊 共 ␮ + tan ␾ 兲艌0, 共9c兲 a0共 ␮ tan ␾ +1兲+ 冉 mg − ␶ rsin ␾ 冊 共 ␮ − tan ␾ 兲艌0. 共9d兲 Some of these conditions are redundant. For example, Eq. 共9c兲implies Eq. 共9b兲. To see this, note that due to our choice ␶ ⬎0, 冉 2 ␶ rsin ␾ 冊 共 ␮ + tan ␾ 兲艌0. 共10兲 Adding this to Eq. 共9c兲gives Eq. 共9b兲. Thus Eq. 共9b兲never needs to be considered; it will be satisfied as long as Eq. 共9c兲 is. Similar reasoning can be used to show that Eq. 共9a兲is redundant when ␮ −tan ␾ ⬎0, and that Eq. 共9d兲is redundant when ␮ −tan ␾ ⬍0. When ␮ =tan ␾ , these conditions are equivalent. The conditions in Eqs. 共9兲put limits on the value that a0 can attain 关4兴. Furthermore, when a contact is sliding, we have equality in one of the above cases, and a0is determined: indeterminacy disappears. C. Application We now attempt to deduce the contact forces and particle motion, using only the considerations presented above. No relation between particle deformation and contact forces is assumed. Note that the biggest difficulty will be determining the coefficient a0; the other three coefficients in the expansion in Eq. 共6兲can be read off from Eq. 共7兲. The unknown value of a0is an expression of the indeterminacy of the problem. We must consider three separate cases, depending on the relation of ␮ to the slope of the sides of the groove, tan ␾ . 1. Shallow slopes: tan ␾ ⬍ ␮ Let us first restrict ourselves to tan ␾ ⬍ ␮ and ␮ ⬍1. Under these assumptions, we have ␮ −tan ␾ ⬎0, so the a0is constrained by the conditions Eqs. 共9c兲and 共9d兲. Furthermore, ␮ tan ␾ −1⬍0, so that Eq. 共9c兲sets a lower, and Eq. 共9d兲an upper bound on a0: amin 艋a0艋amax,共11兲 where amin =− 冋 mg − ␶ rsin ␾ 册 冉 ␮ − tan ␾ 1+ ␮ tan ␾ 冊 , amax = 冋 mg − ␶ rsin ␾ 册 冉 ␮ + tan ␾ 1− ␮ tan ␾ 冊 .共12兲 The quantities in the curved brackets are always positive, while the quantity in the square brackets is positive for ␶ =0, and decreases as ␶ increases. It vanishes when ␶ = ␶ 1⬅mgr sin ␾ ,共13兲 and becomes negative when ␶ ⬎ ␶ 1. When it is negative, no solution for a0exists, for amin⬎amax. Thus there are no solutions of the equations of static equilibrium, when ␶ ⬎ ␶ 1 that also satisfy Eq. 共8兲, and the disk must move. On the other hand, when ␶ ⬍ ␶ 1, an infinite number of such solutions exist, one for each value of a0satisfying Eq. 共11兲. Finally, at ␶ = ␶ 1, we have amin=a0=amax=0; hence there is a unique solution. We have shown that the disk must rotate when ␶ ⬎ ␶ 1, and that static solutions exist for ␶ 艋 ␶ 1, but we have not yet shown that the disk will not move when ␶ 艋 ␶ 1. After all, for some choice of a0, it might be possible to put the disk in motion. But this can be shown to be impossible by considering what happens at ␶ = ␶ 1. In this situation, we have a0 =0, so the contact forces are R ␣ =mg cos ␾ ,T ␣ =−mg sin ␾ ,R ␤ =T ␤ =0. 共14兲 Thus the disk is entirely supported by contact ␣ .R ␣ and T ␣ cancel the gravitational force, and ␶ balances the torque rT ␣ =− ␶ 1exerted on the disk by the contact force. A small increase in ␶ and T ␣ will cause the disk to roll upward to the left, out of the groove. On the other hand, if ␶ ⬍ ␶ 1, then the torque rT ␣ =− ␶ 1will not be balanced, and the disk will try to roll down the slope. Therefore, if ␶ ⬍ ␶ 1, the disk cannot rotate. Another possibility is that the disk rotate in place. But this is impossible, when tan ␾ ⬍ ␮ , because this requires T ␣ = − ␮ R ␣ and T ␤ =− ␮ R ␤ , i.e., contacts ␣ and ␤ must be sliding. This occurs when we have equality in both Eqs. 共9a兲and 共9c兲. But this cannot happen because when equality holds in Eq. 共9a兲, then Eq. 共9d兲is violated. Therefore, the value of a0has no effect on the motion of the disk. If ␶ ⬍ ␶ 1, the disk remains in the groove, although the forces cannot be uniquely determined. At ␶ = ␶ 1, there is a unique solution for the forces: contact ␤ opens and the particle’s weight is supported by contact ␣ . For ␶ ⬎ ␶ 1, the disk rolls out of the groove. Note that there is a relationship between the onset of motion and the disappearance of indeterminacy: we have motion for ␶ ⬎ ␶ 1and indeterminacy for ␶ ⬍ ␶ 1. Indeed, we could have anticipated this when we wrote down Eqs. 共9兲. For a disk to move, at least one contact must be sliding 共or open兲, leading to an equality in at least one of Eqs. 共9兲. Then this equality can be used to determine a0, eliminating the indeterminacy. 2. Intermediate slopes: ␮ ⬍tan ␾ ⬍1/ ␮ Now let us consider intermediate angles ␮ 艋tan ␾ 艋1/ ␮ . At these values of ␾ , the relevant conditions from Eqs. 共9兲 are Eqs. 共9a兲and 共9c兲. Equation 共9c兲sets an upper, and Eq. INDETERMINACY AND THE ONSET OF MOTION IN A…PHYSICAL REVIEW E 72, 021304 共2005兲 021304-3 共9a兲a lower, bound on a0. Thus we have Eq. 共11兲again, except now amin = 冋 mg + ␶ rsin ␾ 册 冉 tan ␾ − ␮ 1+ ␮ tan ␾ 冊 , amax = 冋 mg − ␶ rsin ␾ 册 冉 tan ␾ + ␮ 1− ␮ tan ␾ 冊 .共15兲 Again, the quantities in the curved brackets are positive for ␮ ⬍tan ␾ ⬍1/ ␮ . The quantities in the square brackets are equal and positive for ␶ =0. But as ␶ increases, amin increases while amax decreases. When ␶ = ␶ 2⬅mgr ␮ 1+ ␮ 2sec ␾ ,共16兲 we have amin=a0=amax=a2, where a2=mg 冋 1+ ␮ tan ␾ 1 + tan2 ␾ 1+ ␮ 2 册 冉 tan ␾ − ␮ 1+ ␮ tan ␾ 冊 ,共17兲 and there is a unique solution for the contact forces. When ␶ ⬎ ␶ 2, no static solution satisfying the contact conditions exists, as amin⬎amax. Furthermore, we can show that the disk cannot move when ␶ ⬍ ␶ 2. The disk cannot roll out of the groove, because the forces in Eq. 共14兲violate 兩T ␣ 兩艋 ␮ R ␣ . If the disk is to rotate in place, we must have T ␣ =− ␮ R ␣ and T ␤ =− ␮ R ␤ , i.e., we must have equality in Eqs. 共9a兲and 共9c兲. But we have just showed that requiring these equalities leads to ␶ = ␶ 2, amin=a0=amax. Therefore, we have exactly the same situation as in Sec. II C 1. For ␶ ⬍ ␶ 2, there is indeterminacy without motion, and for ␶ ⬎ ␶ 2there is motion. There is a unique solution for the forces precisely at ␶ = ␶ 2. Furthermore, note that ␶ 1= ␶ 2when tan ␾ = ␮ and when tan ␾ =1/ ␮ . 3. Large slopes: tan ␾ ⬎1/ ␮ At large angles 共tan ␾ ⬎1/ ␮ 兲, the relevant conditions are still Eqs. 共9a兲and 共9c兲as in the previous section. But now the factor 1− ␮ tan ␾ becomes negative, transforming Eq. 共9c兲into a lower bound on a0. Equation 共9a兲remains a lower bound, meaning that there is no upper bound on a0. It can be made arbitrarily large. Therefore, for any ␶ ⬎0, it is possible to satisfy all conditions in Eq. 共9兲simply by making a0very large. 共Note that all the factors multiplying a0in these conditions are positive when tan ␾ ⬎1/ ␮ .兲This means that an infinite torque can be put on the disk without causing it to rotate. On the other hand, if ␶ = ␶ 2, one can obtain equality in Eqs. 共9a兲and 共9c兲, meaning contacts ␣ and ␤ are sliding. If ␶ is increased slightly, then the disk can start to rotate in place. Therefore, both static and moving solutions coexist for ␶ ⬎ ␶ 2, and we cannot decide on the basis of Eqs. 共9兲whether the disk rotates or not. This remarkable situation poses two questions: First, how will the CD algorithm handle this case? This algorithm tries to use the information presented in Sec. II to deduce the motion of the particles, but we have shown that this is insufficient. Second, it would be very easy to construct an experiment to measure the torque needed to make the disk rotate. One could construct a groove with tan ␾ ⬎1/ ␮ , and place a cylinder in the groove, and try to turn the cylinder. What determines whether the cylinder rotate? These questions are investigated in the next sections. III. SIMULATIONS To verify the conclusions of the previous section, and to clarify the situation for tan ␾ ⬎1/ ␮ , we carried out CD and MD simulations. The disk was placed in the groove without any forces acting on it. Then gravity was turned on, and increased, until it reached a maximum of gmax at t=tA. Then it was decreased, reaching g*at t=tB, and thereafter was held constant. A linearly increasing torque ␶ was then applied until the particle began to rotate. Figure 2 shows the gravity force and torque as functions of time. In the following, we will vary gmax, but keep g*fixed. In all cases, ␮ =0.5. The experiment can be considered as a very simple soil mechanics experiment. First, the sample, consisting of one disk, is “prepared” by pushing the disk into the groove, and “loaded” by the torque until it “yields.” We are especially interested in knowing if and how the preparation affects the yielding torque. A typical result of a simulation is shown in Fig. 3. The three nonzero coefficients ay,a ␪ , and a0of the expansion in Eq. 共6兲are shown. In accordance with Eq. 共7兲,ayand a ␪ are identical for the CD and MD solution methods, and given by FIG. 2. The forces applied to the disk. Time is given in units of t*=冑d/g. The gravity attains a maximum at t=tA, and the torque begins to increase at t=tB. FIG. 3. The coefficients of the expansion Eq. 共6兲as a function of time. Here ␾ =20° and ␮ =0.5 so that tan ␾ ⬍ ␮ . MCNAMARA, GARCÍA-ROJO, AND HERRMANN PHYSICAL REVIEW E 72, 021304 共2005兲 021304-4 ay=mg,a ␪ =− ␶ /共rsin ␾ 兲. The coefficient a0, however, differs between the CD and MD solutions. In the following, we will plot only a0, as the other coefficients are always uniquely determined by the imposed forces. Three different values of ␾ will be considered, corresponding to the three subsections of Sec. II C. A. Shallow slopes: tan ␾ ⬍ ␮ In Fig. 4, we plot a0as a function of time when ␾ =20° 共tan ␾ ⬍ ␮ 兲, for the MD and CD simulations. In the MD simulations, a0=0, independent of time and of gmax.Inthe CD simulation, a0is proportional to the gravity for t⬍tB. Note that the values of a0are significant—almost twice the weight of the particle for gmax/g*=8. However, at t=tB, all CD simulations have the same value of a0. As the torque is applied, a0moves linearly towards 0. This occurs because when ␶ 艋 ␶ 1, the disk cannot move, but a0is constrained between the two values given in Eq. 共12兲, which approach each other as ␶ approaches ␶ 1. Indeed, the diagonal line segment near 50⬍t/t*⬍60 traced out by the CD simulations corresponds to requiring equality in Eq. 共9d兲. The conditions Eqs. 共9c兲and 共9d兲act as a “funnel” that guides a0toward 0 as the torque increases, so that a0=0 when ␶ = ␶ 1. Thus the MD and CD methods both predict that the disk moves at ␶ = ␶ 1. However, the two methods predict different routes to failure. MD predicts that both contacts remain nonsliding until the force at ␤ vanishes, and the disk starts to roll. On the other hand, CD predicts that the contact at ␤ first becomes sliding, and then later starts to roll. B. Intermediate slopes: ␮ ⬍tan ␾ ⬍1/ ␮ In Fig. 5, we show a0for the case of ␾ =50° 共 ␮ ⬍tan ␾ ⬍1/ ␮ 兲. The behavior of the CD simulation is quite similar to the preceding case. For t⬍tB,a0is again proportional to gravity, but attaining even larger values than before. At t =tB, all CD simulations are identical, and they evolve together during the loading. They again encounter the “funnel” arising from the lower and upper bounds on a0given in Eq. 共15兲, and move toward a0=a2, and the disk starts to rotate when ␶ = ␶ 2. The MD simulation has a quite different behavior. For t ⬍tA,a0is also proportional to the gravity, although lower than the CD value. But for t⬎tA,a0remains constant, so that a0is proportional to gmax, even at t=tB. In this case, therefore, the preparation does make a difference. The MD simulations remember the value of gmax, and this memory is stored in a0. This simple example shows how the indeterminacy of forces in perfectly rigid particles is related to history dependency in packings of deformable particles. The coefficient a0, which is undetermined in the perfectly rigid case, here contains the memory of the packing, and is simply proportional to the greatest downward force that the disk has experienced in the past. The history dependency of the MD simulation, however, does not affect the yielding torque. The different MD simulations encounter the funnel at different times, but are all guided towards the point ␶ = ␶ 2,a0=a2, where the disk begins to rotate. Therefore, the memory in the MD simulation can only be detected by inspecting the contact forces. C. Large slopes: tan ␾ ⬎1/ ␮ Finally, in Fig. 6, we show the case where ␾ =80° 共tan ␾ ⬎1/ ␮ 兲. This case has some important differences from the preceding cases. First of all, both the CD as well as MD exhibits history dependence, as the contact forces at the end of the sample preparation and the yielding torque depend on gmax. Second, only half of the funnel operates as before. One may ask how the CD simulation can exhibit memory, because it deduces the contact forces from the principles stated in Sec. II but the history of the packing does not enter into those considerations. The memory of the CD algorithm comes from the way it chooses the contact forces. It begins with an initial guess, which it then refines through an iterative process until it arrives at a satisfactory solution. One usually takes the initial guess to be the solution of the previous time step because convergence is faster. But in Fig. 6, it is clear that this choice of initial guess also serves to make FIG. 4. The coefficient a0as a function of time for various values of gmax/g*. Here ␾ =20° and ␮ =0.5, so that tan ␾ ⬍ ␮ . The thin lines are the CD simulations, with gmax/g=2 共lowest curve兲,4, 6, and 8 共highest curve兲. The thick lines are the corresponding MD simulations, which fall on top of each other in this case. The large dot indicates the point where the disk begins to rotate. FIG. 5. The undetermined coefficient a0from Eq. 共6兲as a function of time. Here, ␾ =50° and ␮ =0.5, so that ␮ ⬍tan ␾ ⬍1/ ␮ . The thin lines are the CD simulations, with gmax/g=2 共lowest curve兲,4, 6, and 8 共highest curve兲. The thick lines are the MD simulations with the same values of gmax/g. The large dot indicates the point where the disk begins to rotate. INDETERMINACY AND THE ONSET OF MOTION IN A…PHYSICAL REVIEW E 72, 021304 共2005兲 021304-5 the CD algorithm history dependent, as the solution chosen depends on the initial guess 关5兴. But this memory is not the same as in the MD simulations. Nor is it clear why the preparation of the sample makes a difference at ␾ =80° but not at ␾ =50°. But the most important difference with the previous examples is that only half of the funnel works as before. As was noted earlier, there is no longer an upper bound on a0; Eqs. 共9a兲and 共9c兲are both lower bounds. From the figure, one can see that the lower bound set by Eq. 共9a兲functions as before. When a0reaches this lower bound, it then increases linearly as the torque is increased, until equality is obtained in both Eqs. 共9a兲and 共9c兲. Then the disk begins to rotate. The behavior of the lower bound set by Eq. 共9c兲is quite different. When equality is obtained in that condition, the value of a0 jumps discontinuously, and the particle begins to rotate. As a result, the yielding torque is determined by when the system first satisfies equality in Eq. 共9c兲, and this in turn depends on the value of a0set by the preparation phase of the experiment. Thus in this case, the memory does affect the yielding torque. In Fig. 7, we show the yielding torque as a function of gmax. It is initially independent of gmax, corresponding to the case where the system first meets the lower edge of the funnel. Then, CD and MD give different yielding torques. The MD values are well predicted by a result that will be obtained in Sec. IV D 4. The CD values for the yielding torque are equal to or below the MD ones, because a0decreases as the torque is increased, while in the MD case, a0remains constant. This is plainly visible in Fig. 6. The different CD points correspond to different iteration procedures. The CD algorithm searches for a solution by adjusting the contact forces one by one. It considers a given contact, and calculates the change in its contact forces needed to prevent interpenetration and minimize sliding. It then adds this to change to its current guess for the forces. After passing over all the contacts a certain number of times, the changes in the force needed fall below a certain threshold, and the solution is accepted. But this procedure can be altered by multiplying the calculated change in forces by a number ␭, with ␭=1 corresponding to the usual iteration procedure. Changing ␭changes the yielding torque, as shown above. Furthermore, the memory in the CD algorithm can be erased at each time step by always using F=0as the initial guess, instead of the solution of the last time step. When this is done, one obtains the points labeled “RESET” in Fig. 7 and the yielding torque is independent of gmax. This result emphasizes the important role of memory in this experiment. IV. DEFORMABLE PARTICLES A. Formalism We now present a second method of calculating the contact forces that is able to explain all features of the MD simulations presented in the previous section. The cost of this additional information is that more assumptions must be made. First of all, one must specify how the forces depend on the deformations, but more significantly, one must specify the past history of the packing. This method shows how indeterminacy is replaced by memory. More precisely, when tan ␾ ⬎ ␮ , we show that a0in Eq. 共6兲stores information about the largest downward force that has been exerted on the disk in the past. Indeed, a0is a quantification of the “degree of wedging” discussed in Ref. 关5兴. We model the deformations in a very simple way, assuming that contact forces are generated by springs that are stretched by the motion of the disk. Our model is inspired by the MD simulation method 关8兴, but we are able to apply it analytically to our simple problem. When two bodies touch, a normal and a tangential spring be created at the instant of contact. The contact forces are simply proportional to the spring lengths: R=−k ␦ n,T=−k ␦ t,共18兲 where ␦ nand ␦ tare the normal and tangential spring lengths, and kis the spring constant, here assumed to be equal for both the tangential and normal springs. One normally includes damping in Eq. 共18兲, but we will assume that the motion is quasistatic, i.e., a sequence of equilibrium states. FIG. 6. Coefficients of the expansion Eq. 共6兲as a function of time. Here, ␾ =80° and ␮ =0.5, so that tan ␾ ⬎1/ ␮ . The thin lines are the CD simulations, with gmax/g=2 共lowest curve兲,4,6,and8 共highest curve兲. The thick lines are the MD simulations with the same values of gmax/g. The large dots indicates the points where the disk begins to rotate. The diagonal dashed line is obtained by requiring equality in Eq. 共9c兲. FIG. 7. The torque at which the disk begins to rotate at ␾ =80° and ␮ =0.5, for different values of gmax. The four different sets of points for CD correspond to the different iteration algorithms described in the text. The line shows the theoretical prediction given in Sec. IV D 4. MCNAMARA, GARCÍA-ROJO, AND HERRMANN PHYSICAL REVIEW E 72, 021304 共2005兲 021304-6 Under this assumption, the particle velocities are vanishingly small, and so are the damping forces. Equation 共18兲leads to a vector equation for Fin terms of the spring lengths D: F=−kD,共19兲 where ␦ n, ␣ , ␦ t, ␣ , ␦ n, ␤ , and ␦ t, ␤ are gathered into Djust as the contact forces are arranged in F关see Eq. 共4兲兴. The spring lengths obey ␦ ˙n=vn, ␦ ˙t= 再 vt, ± ␮ vn,兵共20兲 where vnand vtare the normal and tangential components of the relative velocity. There are two choices for ␦ ˙tbecause ␮ R艌兩T兩means that the tangential spring has a maximum allowable length: 兩 ␦ t兩艋 ␮ ␦ n. If applying ␦ ˙t=vtwould lead to a violation of this condition, the second choice is taken. Eq. 共20兲relates the time derivative of Dto the velocity of the disk. If all contacts are nonsliding, we have ␦ ˙n, ␣ =vxsin ␾ +vycos ␾ , ␦ ˙t, ␣ =vxcos ␾ −vysin ␾ +r ␻ , ␦ ˙n, ␤ =−vxsin ␾ +vysin ␾ , ␦ ˙t, ␤ =vxcos ␾ +vysin ␾ +r ␻ .共21兲 Gathering the velocities of the disk into a single vector v= 冢 vx vy ␻ 冣 ,共22兲 Eqs. 共21兲becomes D ˙=cTv,共23兲 where cTis the transpose of cin Eq. 共3兲. The appearance of cTis not a coincidence, but a general property valid for all granular packings 关9兴. It remains to incorporate the status of the contacts into our formalism. This can be be done by inserting a matrix Sthat depends on the contact status into Eq. 共23兲: D ˙=ScTv,共24兲 where S= 冉 S ␣ 0 0S ␤ 冊 .共25兲 Here, S ␣ is the identity matrix if contact ␣ is nonsliding, S ␣ =0 if it is open, and S ␣ = 冉 10 ± ␮ 0 冊 共26兲 if it is sliding. Finally, let us not forget Newton’s equation of motion: Mv ˙=cF +fext.共27兲 Here, the matrix Mcontains the mass mand moment of inertia Iof the disk on the diagonal: M= 冢 m00 0m0 00I 冣 .共28兲 Equations 共19兲,共24兲, and 共27兲form a set of equations that can be solved for the motion of the disk. Note that cremains constant because all contacts are between the disk and the straight walls. The status of the contacts can change, however, so it is useful to divide time up into segments 关t0,t1兴,关t1,t2兴,..., where the changes of contact status occur at the times t1,t2, ... . In the interior of a time interval, the matrices c,S, and Mare constant, so Eqs. 共19兲,共24兲, and 共27兲can be combined into Mv ¨=−Qv +f ˙ext,共29兲 where the matrix Q=kcScTis called the stiffness matrix and gives the changes in the forces that arise from a small displacement of the disk. The MD simulations used to obtain the results in Sec. III are numerical solutions of Eq. 共29兲. Those simulations concerned very slow motion of the disk, and thus we make the quasistatic approximation, and assume that the disk is always in an equilibrium state. Thus the left hand side of Eq. 共29兲vanishes: Qv =f ˙ext.共30兲 This is the equation that we will analyze in the rest of this paper. It is capable of explaining the results obtained in Sec. III. B. Relation to other theories Before proceeding, several remarks on Eq. 共30兲are in order. Let us begin by pointing out that the derivation of Eq. 共30兲can be generalized to packings with a large number of particles. In this case, the matrices cand Scan be constructed 关4兴. It is also possible to generalize to the case where tangential and normal motions have differing stiffness, or even where each contact has different stiffness. We therefore expect that the results presented here can be generalized to packings of many particles. Let us also note that a numerical method, sometimes called the “granular element method” 关10兴, is based on Eq. 共30兲. It is less common than MD and CD because one must rebuild the stiffness matrix Qevery time a contact opens, closes, or becomes sliding. Nevertheless, it has been compared with MD simulations and used to investigate the quasistatic deformation of granular materials 关10,11兴. Next, we point out that Eq. 共30兲can be written in terms of increments instead of time derivatives: Q ␦ u= ␦ fext.共31兲 Here, ␦ uis a small increment of displacement caused by an incremental change in the applied load ␦ fext. The difference between this equation and Eq. 共30兲is that all reference to INDETERMINACY AND THE ONSET OF MOTION IN A…PHYSICAL REVIEW E 72, 021304 共2005兲 021304-7 time has been removed from Eq. 共31兲. This is a consequence of the quasistatic approximation: since the system is always in mechanical equilibrium, it does not matter how it moves from one state to another 共as long as it is done slowly enough兲. Finally let us point out several parallels between Eqs. 共30兲 and 共31兲and the theories of elasticity and elastoplasticity in continuum mechanics. Equation 共31兲resembles Hooke’s law in elasticity, with Qplaying the role of the elasticity tensor, and giving a linear relation between the forces and the displacements. However, in continuum mechanics, the forces and displacements are given by the stress and strain tensors, which are second rank tensor functions of position. In Eq. 共31兲, the analogous quantities are ␦ uand ␦ fext. Furthermore, as long as there are no sliding contacts, all deformations are reversible—they can be removed simply by removing the applied force, just as in elasticity. When a contact becomes sliding, however, the analogy to elasticity breaks down, energy will be dissipated at the sliding contact, and permanent deformation will occur. However, one can then consider an analogy to elastoplasticity 关12,13兴, which was developed to describe the irreversible deformation of metals. In this theory, the stress is considered as a point in a space whose coordinates are the components of the stress tensor. As the stress changes, this point moves correspondingly through stress space. Elastoplasticity posits that the origin of stress space is enclosed by a “yield surface.” As long as the system remains within the yield surface, the material behaves elastically—all deformation vanishes if the stress returns to 0. But when the system arrives at the yield surface, the material deforms “plastically” 共irreversibly兲as well as elastically. As the stress increases further, the yield surface expands as well, so that the system remains always on the yield surface. If the stress is decreased, the system moves inside the yield surface, and the deformations are again elastic. Although not obvious from Eq. 共30兲, the description put forth here behaves in a similar way. Requiring equality in Eq. 共9兲defines a “yield surface” in the four-dimensional space 共R ␣ ,T ␣ ,R ␤ ,T ␤ 兲. As long as the system remains inside this surface, the disk’s motions are reversible and elastic. Once the system arrives at this surface, irreversible motion occurs. Another parallel to elastoplasticity occurs if the motion is reversed: in this case, the sliding contact becomes nonsliding, and the deformation is again elastic. This analogy could be studied further, but this is beyond the scope of this paper. C. The onset of motion In deriving Eq. 共30兲, we assumed that the motion is quasistatic, so that the left hand side of Eq. 共29兲containing the inertia of the disk could be neglected. However, we are also interested in the onset of motion, i.e., when a static configuration yields and begins to move. There are two different ways motion can begin. The first way occurs when the stiffness matrix Qhas a null eigenvector, that is, there is a vector v*⫽0 such that Qv*=0. The motion given by v*is called a “mechanism” 关9兴. When there is a mechanism, Qmaps R3into a twodimensional subspace, and if f ˙ext is not in this subspace, Eq. 共30兲has no solution. This means that the inertia terms must be included, i.e., Eq. 共29兲must be used instead. Equation 共29兲always has a solution, since Mis a diagonal matrix with positive entries, and thus maps R3onto itself. If t*is the time when Qbecomes singular 共perhaps due to a change in contact status兲. Then for times tjust after t*, the solution to Eq. 共29兲has the form v⬃共t−t*兲2M−1f ˙*,共32兲 where f ˙*is the part of f ˙ext for which no solution to Eq. 共30兲 exists. We will call this “motion through a mechanism.” It is characterized by accelerations of order the imposed external force divided by the inertia. The second way to motion occurs when the quadratic form vTQv becomes negative. To see this, we dot Eq. 共29兲 through by vand obtain vTMv ¨=−vTQv +v·f ˙ext.共33兲 We must have v·f ˙ext⬎0, because when the external forces are changed, the disk must move in the direction of the change in the force. In the quasistatic approximation, the left hand side of Eq. 共33兲is set to zero. Thus to obtain equality, we must have vTQv ⬎0.共34兲 This inequality must hold, otherwise the quasistatic approximation cannot be made. In the context of elastoplasticity, this condition is called Hill’s stability criterion 关14兴. Let us now investigate the motion just after Eq. 共34兲has been violated. As before, we let t*be the time when the motion begins. Just after t*, the motion will be given by Eq. 共32兲, since the term Qv requires a finite time to change. Soon, however, this term dominates the imposed forces, so that Eq. 共29兲becomes Mv ¨=−Qv. The solution is then v⬃v*et冑Q/M,共35兲 where v*is the vector for which Q=−共v* TQv*兲/共v*·v*兲is a maximum, and M=共v* TMv*兲/共v*·v*兲. Thus, when Eq. 共34兲is violated, Qacts like a negative number in Eq. 共29兲, and the contact forces then amplify the velocities, which then grow exponentially. In this case, the packing fails catastrophically, with the contact forces changing much more quickly than fext. We call this “motion through an instability.” It is characterized by a time scale that is much shorter than motion by a mechanism. D. Application 1. Calculation of displacements Let us now show how the contact forces can be calculated without indeterminacy. We calculate the coefficients of the expansion in Eq. 共6兲. We can construct an equation for these coefficients by differentiating Eq. 共19兲with respect to time, and combining it with Eq. 共24兲, leading to F ˙=ScTv. Then we project this equation onto the basis given in Eq. 共5兲by leftMCNAMARA, GARCÍA-ROJO, AND HERRMANN PHYSICAL REVIEW E 72, 021304 共2005兲 021304-8 multiplying it by a,a4⫻4 matrix whose rows are the basis vectors in Eq. 共5兲, but without the factor of 1/2. After doing this, we obtain 冢 a ˙x a ˙y a ˙ ␪ a ˙0 冣 =kaScTv.共36兲 Equation 共36兲represents four equations, one corresponding to each of ax,ay,a ␪ , and a0. The first three coefficients are known from Eq. 共7兲; thus the first three equations above can be used to find the three components of v. Then the last equation can be used to determine a ˙0. The displacement of the disk is obtained by integrating v. 2. Shallow angles: tan ␾ ⬍ ␮ When tan ␾ ⬍ ␮ , the contacts remain nonsliding until the disk begins to rotate. Therefore, we can consider the entire experiment to take place during one time interval. When the contacts are nonsliding, kaScT=2k 冢 10rcos ␾ 010 00rsin ␾ 000 冣 .共37兲 Using this result in Eq. 共36兲, we have a ˙x=2kvx+2kr ␻ cos ␾ , a ˙y=2kvy, a ˙ ␪ =2kr ␻ sin ␾ , a ˙0=0. 共38兲 From the last line, we can see that a0remains constant as long as the contacts are nonsliding. This explains why a0 =0 for all the MD simulations in Fig. 4: a0=0 at the beginning of the simulation, and thus remains so until the disk rolls. This also explains why a0is constant when gis decreased in Figs. 5 and 6. Next let us calculate the disk’s motion. Equations 共38兲can be solved for the velocities and integrated. In this way, the displacement uof the disk can be shown to be u=1 2k 冢 0 −mg ␶ /共r2sin2 ␾ 兲 冣 .共39兲 Because no contact changes its status, the motion is perfectly reversible. This is another way to see that the packing has no memory when tan ␾ ⬍ ␮ . Finally, let us analyze the onset of motion. For tan ␾ ⬍ ␮ , motion occurs through a mechanism, as can be seen by constructing the stiffness matrix Qat the time when the motion begins. It was shown in Sec. II C 1, Eq. 共14兲that the contact ␤ opens when the motion begins, while contact ␣ remains nonsliding. The submatrices making up Sin Eq. 共25兲 are thus S ␣ = 共 10 01 兲 and S ␤ =0. This leads to Q=k 冢 10rcos ␾ 01−rsin ␾ rcos ␾ −rsin ␾ r2 冣 .共40兲 There is a vector v*⫽0such that Qv*=0, namely, v*= 冢 −rcos ␾ rsin ␾ 1 冣 .共41兲 The motion given by v*is the particle rolling upward to the left out of the groove, as described in Sec. II C 1. 3. Intermediate slopes: ␮ ⬍tan ␾ ⬍1/ ␮ When tan ␾ ⬎ ␮ and t⬍tA, the contacts are sliding with T ␣ =− ␮ R ␣ and T ␤ = ␮ R ␤ : kaScT=2kcos2 ␾ 冢 tan2 ␾ − ␮ tan ␾ 00 01+ ␮ tan ␾ 0 − tan ␾ − ␮ tan2 ␾ 00 0tan ␾ − ␮ 0 冣 . 共42兲 Changes in contact status will occur, so it is necessary to subdivide the experiment into time intervals. The first change of status occurs at t=tA, when the derivative of the gravity changes sign, so we define our first interval to end at t=tA. Calculating the particle displacements using the same method as before yields uy共tA兲=−mgmax 2k sec2 ␾ 1+ ␮ tan ␾ .共43兲 as before, ux共tA兲=u ␪ 共tA兲=0. We also have a0共tA兲=amem =mgmax tan ␾ − ␮ 1+ ␮ tan ␾ .共44兲 When the gravitational force begins to decrease, all contacts are nonsliding. The displacements must be calculated using the matrix given in Eq. 共37兲. The particle rises by a distance m共gmax−g*兲/共2k兲, while a0remains constant, as in Eq. 共38兲. The particle’s position is now uy=−mg*+amem 2k.共45兲 It is important to realize that the particle is amem/共2k兲lower than it would be if the gravity had been increased monotonically from 0 to g*. This displacement thus represents a “memory” of the forces exerted on disk in the past. But it is more convenient to think of the memory as being stored in a0, for the displacement is proportional to 1/kand very small if the particle is stiff, whereas a0is proportional to gmax and independent of the current value of g. When the torque is applied, the contacts are still nonsliding, and they remain so until one of the contacts becomes INDETERMINACY AND THE ONSET OF MOTION IN A…PHYSICAL REVIEW E 72, 021304 共2005兲 021304-9