text:
# *

## 1 Introduction

Traditionally, the subgrade of railway tracks is modeled using continuum mechanics. These methods have been proven to yield reliable results in many applications and have been developed to standard methods. In certain applications, however, continuum models fail in describing the mechanical properties of the subgrade. This is the case when the ballast must not be considered as a continuous medium, but when the granularity of the material is important. Typical processes which cannot be sufficiently explained by continuum models are sedimentation of the ballast due to recurrent load, abrasion of the ballast particles which leads to less efficient damping properties and the formation of force chains inside the ballast material.

In the past decade mainly by physicists much work has been done in the field of Molecular Dynamics of granular material, i.e., the numerical simulation of granular material as many-particle systems. This technique was applied to many interesting systems and has contributed to the explanation of several exciting and technologically important effects, such as mixing and demixing of granular materials, avalanche statistics on sand heaps, milling processes, convection dynamics in shaken granular materials and others. Many examples of such simulations can be found, e.g., in .

The idea of Molecular Dynamics is to simulate the granular material as a many particle system and to determine the dynamics of the system by numerical integration of Newton’s equation of motion for each of the $`N`$ particles:

$$\begin{array}{cc}\hfill \ddot{\stackrel{}{r}}_i=& \stackrel{}{F}_i/m_i\hfill \\ \hfill \ddot{\varphi }_i=& \widehat{J}_{i}^{}{}_{}{}^{1}\stackrel{}{M}_i,\hfill \end{array}$$
(1)
where $`\stackrel{}{r}_i`$ and $`\stackrel{}{\varphi }_i`$ are the position and the orientation of the $`i`$-th particle of mass $`m_i`$ and moment of inertia $`\widehat{J}_i`$ while $`\stackrel{}{F}_i`$ and $`\stackrel{}{M}_i`$ are the force and the torque acting on this grain. In three dimensions Eqs. (1) establish a set of $`6N`$ coupled non-linear differential equations which have to be numerically integrated. The force $`\stackrel{}{F}_i`$ consists of gravity and of the interaction force of the particle $`i`$ with other particles $`j`$

$$\stackrel{}{F}_i=m_i\stackrel{}{g}+\underset{j}{}\stackrel{}{F}_{ij},$$
(2)
where $`\stackrel{}{g}`$ is the gravitational acceleration. Granular particles interact with each other only if they are in mechanical contact. For spheres of radii $`R_i`$ and $`R_j`$ we write

$$\stackrel{}{F}_{ij}=\{\begin{array}{cc}F^N_{ij}\stackrel{}{n}_{ij}+F^T\stackrel{}{t}\hfill & \text{if }|\stackrel{}{r}_i\stackrel{}{r}_j|<R_i+R_j\hfill \\ \text{0}\hfill & \text{else\hspace{0.17em},}\hfill \end{array}$$
(3)
where $`F^N`$ and $`F^T`$ are the components of the force in normal and tangential direction with respect to the inter-center vector $`\stackrel{}{r}_i\stackrel{}{r}_j`$ and $`\stackrel{}{n}`$ and $`\stackrel{}{t}`$ are the corresponding unit vectors. There exist several models for the interaction forces in normal and tangential direction $`F_{ij}^{N/T}(\stackrel{}{r}_i,\stackrel{}{r}_j,\dot{\stackrel{}{r}}_i,\dot{\stackrel{}{r}}_j)`$ (see, e.g. ) which shall not be discussed in detail here.

## 2 Molecular Dynamics fails for the simulation of railway ballast

There are many examples where granular systems have been simulated by Molecular Dynamics, however, except for few examples, realistic simulations have been achieved only for systems where the dynamical behavior of the grains dominates the system properties. When the static properties of the particle system become important, i.e., when the relative velocities are small or zero, the details of the interaction force become essential for understanding the system behavior. We are faced with two major problems:

* The interaction force of contacting particles must be known as a function of the particle positions and velocities. In the case of sharply edged grains such as railway ballast, this function is unknown.
* As soon as the realistic simulation of static properties matters for the system behavior, the simulation slows down extremely which implies that to achieve affordable computation times one has to make simplifying assumptions on the particle contact which are not justified from the point of view of mechanics and material science.

For these reasons we believe that the described Molecular Dynamics method is in principle unsuitable for the simulation of railway ballast. Below we list arguments which support this statement:

1. Elastic normal force: Whereas the normal elastic force of contacting spheres is given by the Hertz law

$$F_{ij}^{N,\mathrm{el}}=\frac{2}{3}\frac{Y\sqrt{R_{\mathrm{eff}}}}{1\nu ^2}\xi _{ij}^{3/2}$$
(4)
with $`Y`$, $`\nu `$, and $`R_{\mathrm{eff}}`$ being the Young modulus, the Poisson ratio, and the effective radius, respectively, this force is not known for more complicatedly shaped particles. For smooth particles (when the local radius of the contact area is large as compared to the compression $`\xi _{ij}R_i+R_j\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|`$) the function (4) is certainly a good approximation, for sharply edged particles, however, this function fails. Some authors assume that the normal force is proportional to what they call “overlap”, i.e., the volume of the compressed material or (yet more simple) in two dimensions the “compression length” (for spheres the value $`\xi `$) which is certainly incorrect even for the most simple case of contacting spheres and the more for more complicatedly shaped grains.
2. Dissipative normal force: The dissipative normal force of contacting bodies is unknown, in general. For viscoelastically colliding spheres and other smooth contacting surfaces it is given by

$$F_{ij}^{N,\mathrm{diss}}=A\dot{\xi }\frac{}{\xi }F^{N,\mathrm{el}}.$$
(5)
The prefactor $`A`$ is a complicated function (for details see ) which contains the viscous constants of the material which are unknown in general. For particles which are not smooth, such as railway ballast, even the functional form of the dissipative force is unknown. The frequently applied force law $`F_{ij}^{N,\mathrm{diss}}\dot{\xi }`$ is not justified and fails even for spheres.
3. Tangential force: Whereas for smooth bodies the normal force can be determined from bulk properties of the material, the tangential force is determined by the bulk and by surface properties. A natural (phenomenological) assumption is

$$F_{ij}^T\mu F_{ij}^N,$$
(6)
where $`\mu `$ is the Coulomb friction parameter. Unfortunately, this model is not sufficient for static systems: Assume a (non-spherical) particle which rests on an inclined plane (angle $`\alpha `$). Its normal force is $`F^N=mg\mathrm{cos}\alpha `$ and there is a corresponding tangential force too. To prevent the particle from sliding along the plane one has to assume an additional force which mimics static friction, e.g. . This force cannot be derived from material or surface properties, hence, it is arbitrary. Choices for this function which can be found in literatures are even in disagreement with basic mechanics.
4. Numerical integration of Newton’s law: Besides the fact that the numerical integration of the set (1) for a relevant number of particles $`N`$, e.g., $`N=1,000`$ which establish a quite small system of $`10\times 10\times 10`$ particles in three dimensions, requires serious numerical effort, there are principal problems when applying Molecular Dynamics to a particle system:

* For very rigid particles the gradient of the force is large and, therefore, requires a very small integration time step. The more rigid the particle material the slower progresses the simulation. Assuming a time step of $`\mathrm{\Delta }t=10^7`$ sec which is used in many MD simulations of granular material, one needs $`10^9`$ integration steps of the equations of motion (1) to achieve a real time of $`100`$ sec only. Hence, a simulation of long time behavior using MD seems to be unrealistic.
    * The accuracy of the algorithm cannot be increased too much by reducing the time step. Frequently, predictor-corrector algorithms are used, which require powers of the time step, $`(\mathrm{\Delta }t)^2=10^{14}`$, $`(\mathrm{\Delta }t)^3=10^{21}`$, etc. (again for $`\mathrm{\Delta }t=10^7`$). If we ignore the prefactor of this powers in the integration scheme to obtain a crude estimate, one has to sum numbers of the orders $`O(1)`$, $`𝒪\left(10^7\right)`$, $`𝒪\left(10^{14}\right)`$, etc. Double precision numbers have a mantissa of 15 digits, therefore, one cannot sum numbers which are different by more than $`10^{15}`$. This limits the minimal value of the integration time step. If one requires higher precision one needs real numbers in quadruple precision or higher which implies higher memory consumption and (more importantly) yet larger computation time since multiplication of two numbers in quadruple precision requires approximately quadruple time as a multiplication of double precision numbers.

    We are aware that these numbers are a very crude estimate which have been given to illustrate the problem. Of course, more sophisticated integration schemes suffer less critical from the addition of different numbers. In principle, however, the problem persists.

## 3 Rigid Body Dynamics

In Molecular Dynamics simulations the trajectories of the particles are determined by numerical integration of Newton’s equations of motion. As discussed above, this implies the “soft particle assumption”, i.e., the particle deform under load. Since the particles are very hard but not completely rigid we encounter serious numerical problems as described above.

The Rigid Body Dynamics originates from the opposite idea: The interaction forces are determined from the required behavior of the particles. This method is, therefore, suited to simulate perfectly rigid particles without the necessity to specify a certain force-deformation law. We consider the example of a rigid sphere which rests on a rigid flat surface (see Fig. 1).

There is a point-like contact of the sphere and the surface, hence, there is a contact force $`F^N`$ in vertical direction. Moreover gravity acts on the sphere causing a force $`mg`$. If one would chose the contact force $`F^N=0`$, the sphere would move downwards with an acceleration $`g`$, i.e., it would penetrate the surface. This unphysical behavior has to be avoided by the proper choice of the force $`F^N`$ due to the following conditions:

1. Contact forces have to be chosen in order to avoid penetration of contacting particles.
For our example $`F^Nmg`$ follows. On the other hand, a force $`F^N>mg`$ would lead to unphysical behavior since the sphere would move upwards. This is avoided by the second condition.
2. A contact force vanishes when the contact breaks.
A contact is said to break if the normal component of the relative acceleration of the concerned particles, or their normal velocity is larger than zero. (The relative velocity is counted positive if the particles separate from each other.)
3. There are no attractive normal forces.
In our example $`F^N>mg`$ causes the contact to break which implies that the contact force vanishes. Therefore, from these conditions follows $`F^N=mg`$. For this choice the total force is zero and the sphere rests on the plane. The conditions 1.-3. are sufficient to describe any particle system as long as there are no friction forces. For systems with friction we need one more condition:
4. Friction forces act in parallel with the contact plane, i.e., perpendicular to $`F^N`$. Given the tangential force $`F^{}`$ which is necessary to keep two particles from sliding. Then the acting tangential force is $`\left|F\right|=\mathrm{min}(\left|F^{}\right|,\left|\mu F^N\right|)`$. Its sign has to be chosen opposite to the tangential relative acceleration (or the tangential relative velocity).
In agreement with Coulomb’s friction law the particles slide only if $`\left|F^T\right|\left|F^N\right|`$. If this condition is fulfilled the friction force adopts its maximal value $`\pm \mu F^N`$.

To perform simulations one has to derive the forces of the particles in normal and tangential directions from these four conditions. The corresponding algorithm will be explained in sections 5 and 6.

In our simple example (Fig. 1) there exists only one contact between the sphere and the plane. In more complex situations, e.g. for a resting cube, there are contacts areas instead of points. These contacts may be always reduced to point contacts. It will be shown that the described conditions are sufficient to determine the forces and torques which act on the particles, provided there are not too many contacts in the system. If the number of contacts is too large only the total force and the total torque which act on a particle may be determined, but not each of the pairwise contact forces. For the computation of the particle trajectories, however, the total forces and torques are sufficient. We will return to this issue in Sec. 8.

Rigid Body Dynamics has been intensively studied in the past two decades. Descriptions of the algorithm can be found in . The core of the algorithm is the numerical computation of the contact forces which is a Linear Complementarity Problem . An efficient algorithm for this type of problems can be found in . Rigid Body Dynamics has also been applied to granular systems, e.g. , where frictionless smooth spheres have been simulated. Systems of granular particles subject to friction have been studied, e.g., in .

Due to our understanding the Rigid Body Dynamics is much better suited for the simulation of railway ballast for the following reasons:

* Ballast particles are irregularly shaped and sharply edged. Even if the bulk material properties were precisely known, the contact force law is unknown due to the complicated shape.
* Ballast particles are very stiff which implies that the gradient of the interaction force is very steep. In this regime the numerical integration of Newton’s equation is problematic. For the dynamics of the system the deformation of single particles is unimportant, i.e., the Rigid Body assumption is well justified.
* Static friction, whose treatment in MD-simulations is problematic too, is essential for the dynamics of the system. It is correctly modeled in Rigid Body Dynamics.
* The long time behavior of railway ballast is affected by abrasion and fragmentation of particles. Molecular Dynamics of fragmenting particles may cause artifacts for several reasons which cannot be discussed here in detail (see ). Rigid Body Dynamics is very well suited for this case.

## 4 Schedule of Rigid Body Simulations

The state of the granular system is described by the position and orientation of its particles and by the according time derivatives. Contacts between the particles may be classified as sticking and sliding contacts. In the due of time the contact network is modified by creation and termination of contacts as well as by transformation of sticking contacts into sliding ones and vice versa. Whenever the contact network is modified the state of the system changes qualitatively. The simulation proceeds in discrete time steps. Each of them consists of

1. Contact detection: all existing contacts are registered.
2. Treatment of collisions: A collision takes place if two contacting particles move with negative normal relative velocity. In this case one cannot determine a finite contact force which would avoid penetration of the particles since any force, however large it is, would need a short but finite path to retard the colliding particles. Hence, mutual penetration would be unavoidable. Therefore, we need a special treatment for collisions (see Sec. 7).
3. Formulation of the geometry equations: After a collision, in general, there is a number of contacts of particles which have a positive normal relative velocity, i.e., the particles lose contact. These contacts have to be erased from the list of contacts. The normal components of the relative velocities at all remaining contacts are zero. We have to establish the geometry equations which contain the information about the geometry of the system (see Sec. 5).
4. Computation of the forces: Section 6 deals with the computation of the relative accelerations by means of the geometry equations.
5. Integration of the equations of motion: Finally we have to integrate the equations of motion for all particles. During this operation it may be necessary to update the geometry equations and to repeat the computation of the forces according to the integration scheme used.

## 5 Mathematical description

For convenience at first we will restrict to frictionless particles. When the mathematical framework has been developed for this simplified case we will then introduce friction forces between particles. The rigidity of the particles is enforced by means of mathematical motion constraints of the form

$$g(q)0,$$
(7)
where $`q`$ is a vector containing the positions (center of mass position and orientation) of all particles of the system. The constraint function $`g`$ shall be zero if particles in the system are in contact and larger than zero otherwise. For spheres, which is the most simple case, the constraint function reads

$$g(q)=\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|R_iR_j.$$
(8)
If the spheres would deform each other ($`\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|<R_i+R_j`$) the function $`g(q)`$ would be negative, if the particles touch each other it would be zero. To prevent the deformation we require $`g(q)`$ to be positive or zero, i.e., it has to fulfill the condition (7). For sharp edged particles, such as particles which are described by polyhedrons or polygons, the motion constraints are

$$g(q)=\stackrel{}{n}_j\left(\stackrel{}{r}_i+\stackrel{}{x}_i\stackrel{}{r}_j\right)d,$$
(9)
where $`\stackrel{}{n}_j`$ is normal of the face of particle $`j`$ which is in contact with an edge of the particle $`i`$. The vectors $`\stackrel{}{r}_i`$ and $`\stackrel{}{r}_j`$ are the center of mass positions of the particles, the edge of particle $`i`$ that is in contact with the particle $`j`$ is at position $`\stackrel{}{r}_i+\stackrel{}{x}_i`$. The constant $`d`$ is the distance of the contacting face of $`j`$ from the center of the particle. The left picture in Fig. 2 shows a sketch of a typical contact of two particles. Face-face contacts can be described by two face-edge contacts (see Fig. 2).

Every motion constraint corresponds to a scalar contact force $`f`$. According to d’Alembert’s principle the direction of the contact force is given by the spacial derivative of $`g`$ with respect to all components of the coordinate vector $`q`$, namely $`g/q`$. The contact force that acts on a certain particle $`i`$ is $`fg/q_i`$, with $`q_i`$ being the coordinates of particle $`i`$. We can formulate the equation of motion for the particles111Particle indices are written in Latin letters, contact indices in Greek letters.

$`\widehat{M}_i\ddot{q}_i`$ $`=`$ $`Q_i+{\displaystyle \underset{\alpha }{}}f_\alpha {\displaystyle \frac{g_\alpha }{q_i}}`$ (10)
$`g_\alpha (q)`$ $``$ $`0.`$ (11)
$`\widehat{M}_i`$ is the mass matrix of the particle $`i`$, which has the form

$$\widehat{M}_i=\left(\begin{array}{cccc}\hfill m_i& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill m_i& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill m_i& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill \widehat{J}_i\end{array}\right),$$
(12)
where $`\widehat{J}_i`$ is the moment of inertia tensor. In two dimensional systems $`\widehat{J}_i`$ is only a scalar $`J_i`$ and there are only 2 entries of $`m_i`$. $`Q_i`$ finally is the external force (and torque) acting on particle $`i`$. This is usually gravitation, but other external forces can be incorporated at this point as well.

Although we formulated the motion constraints in the form $`g(q)0`$ to allow separation of particles, contact forces can only act if particles are actually in contact. Therefore, constraints which are strictly positive, i.e. $`g(q)>0`$ (the particles are separated), can be disregarded. These constraints are said to be inactive. The remaining constraints, the active ones, are thus satisfied by $`g(q)=0`$. Since the $`g(q)`$ have to remain non-negative their time derivatives have to be non-negative as well. We therefore have

$`\dot{g}_\alpha `$ $`=`$ $`{\displaystyle \frac{g_\alpha }{q_k}}\dot{q}_k0`$ (13)
$`\ddot{g}_\alpha `$ $`=`$ $`{\displaystyle \frac{g_\alpha }{q_k}}\ddot{q}_k+{\displaystyle \frac{^2g_\alpha }{q_kq_l}}\dot{q}_k\dot{q}_l0.`$ (14)
For simplicity of notation we used the Einstein convention, i.e. summation over doubly occurring indices $`k`$ and $`l`$ is implied. These time derivatives are the relative velocity and relative acceleration of the particles at their contact points. It is important to note that $`\dot{g}_\alpha `$ and $`\ddot{g}_\alpha `$ are not the relative velocity or acceleration of the centers of mass of the particles but of the points of both particles which are actually in contact. It can easily happen that the relative velocity or acceleration of the contact points are positive (the particles are about to separate) although their centers of mass approach each other.

We insert the equation of motion (10) into (14) and find

$`\ddot{g}_\alpha `$ $`=`$ $`{\displaystyle \frac{g_\alpha }{q_k}}\left[\widehat{M}_k^1\left(Q_k+{\displaystyle \underset{\beta }{}}f_\beta {\displaystyle \frac{g_\beta }{q_k}}\right)\right]+{\displaystyle \frac{^2g_\alpha }{q_kq_l}}\dot{q}_k\dot{q}_l`$ (15)
$`=`$ $`{\displaystyle \frac{g_\alpha }{q_k}}\left(\widehat{M}_k^1Q_k\right)+{\displaystyle \frac{^2g_\alpha }{q_kq_l}}\dot{q}_k\dot{q}_l+{\displaystyle \frac{g_\alpha }{q_k}}\widehat{M}_k^1\left({\displaystyle \underset{\beta }{}}f_\beta {\displaystyle \frac{g_\beta }{q_k}}\right).`$ (16)
The first term on the right hand side describes the action of the external forces, the second term describes the action of inertial forces as, e.g., centrifugal force and Coriolis force, the third term finally describes the action of the contact forces. We can rewrite this equation as

$$\ddot{g}_\alpha =b_\alpha +\underset{\beta }{}A_{\alpha \beta }f_\beta ,$$
(17)
with $`A_{\alpha \beta }`$ and $`b_\alpha `$ abbreviating

$$\begin{array}{cc}\hfill A_{\alpha \beta }& =\frac{g_\alpha }{q_k}\widehat{M}_k^1\frac{g_\beta }{q_k}\hfill \\ \hfill b_\alpha & =\frac{g_\alpha }{q_k}\left(\widehat{M}_k^1Q_k\right)+\frac{^2g_\alpha }{q_kq_l}\dot{q}_k\dot{q}_l.\hfill \end{array}$$
(18)
From now on we will denote the relative acceleration of the contacting particles at their contact points – the contact acceleration – as $`a_\alpha `$ instead of $`\ddot{g}_\alpha `$. Equation (17) turns into

$$a_\alpha =b_\alpha +A_{\alpha \beta }f_\beta ,$$
(19)
where again summation over $`\beta `$ is implied. We will call this equation the geometry equation. By means of this equation and the consistency conditions introduced in Sec. 3 we can now determine the contact forces $`f_\beta `$. The consistency conditions read

$$\begin{array}{cc}\hfill a_\alpha & 0\hfill \\ \hfill f_\alpha & 0\hfill \\ \hfill a_\alpha f_\alpha & =0.\hfill \end{array}$$
(20)
The first condition prevents deformation of particles, the second one excludes attractive forces and the third one requests that contact forces may only act if the particles stay in contact, i.e., if $`a_\alpha =0`$. These conditions together with the geometry equation (19) allows to determine the unknown contact forces $`f_\alpha `$. The whole system consisting of Eq. (19) and the conditions (20) is called a Linear Complementarity Problem. It can be solved by Dantzig’s algorithm .

To incorporate friction we introduce additional motion constraints which shall, if possible, impede a tangential motion of the contacting particles. For polygonal particles they are of the form

$$g(q)=\stackrel{}{t}_j\left(\stackrel{}{r}_i+\stackrel{}{x}_i\stackrel{}{r}_j\stackrel{}{x}_j\right).$$
(21)
This constraint has a similar form as the constraint of the normal motion of the particles (9) but instead of the normal unit vector of the contacting face of particle $`j`$ the tangential unit vector $`\stackrel{}{t}_j`$ appears, thus ensuring that the edge of particle $`i`$ at position $`\stackrel{}{r}_i+\stackrel{}{x}_i`$ does not move along the face of particle $`i`$ away from the point $`\stackrel{}{r}_j+\stackrel{}{x}_j`$ of initial contact.

These motion constraints are, however, of different nature than the normal motion constraints. Whereas in the case of the normal motion the constraints must never be violated, the constraints on the tangential motion may actually be violated, as it happens when the particles start to slide. This is due to the fact that the magnitude of the friction forces are limited by $`\mu f_N`$, with $`\mu `$ being the friction constant and $`f_N`$ the corresponding normal force. As reflected by the consistency condition 4 (ref. Sec. 3) the friction force must adopt its maximum value if particles actually slide, i.e. if the tangential motion constraint is inactive. Thus, in this case the value of the friction force is determined without need of further consideration. In contrast to the case of normal motion constraints we may not neglect the inactive constraints because now the corresponding contact forces are non-zero and thus the constraint has to be kept in consideration in order to determine the direction of the tangential force.

Since friction causes only further motion constraints there is, in principle, no need of further discussion of the problem. The geometry equation (19) can describe systems with friction as well. For simplicity of notation it is worth, however, to consider normal and friction forces and their corresponding motion constraints separately. We now have 2 classes of motion constraints

$$\begin{array}{cc}\hfill g^N(q)& 0\hfill \\ \hfill g^T(q)=& 0\hfill \end{array}$$
(22)
and the corresponding contact forces $`f^N`$ and $`f^T`$. Now the equation of motion reads

$$\widehat{M}_k\ddot{q}_k=Q_k+\underset{\alpha }{}\left(f_\alpha ^N\frac{g_\alpha ^N}{q_k}+f_\alpha ^T\frac{g_\alpha ^T}{q_k}\right).$$
(23)
For the second time derivative of the motion constraints we obtain

$$\begin{array}{cc}\hfill \ddot{g}_\alpha ^N=& \frac{g_\alpha ^N}{q_k}\left(\widehat{M}_k^1Q_k\right)+\frac{^2g_\alpha ^N}{q_kq_l}\dot{q}_k\dot{q}_l+\frac{g_\alpha ^N}{q_k}\widehat{M}_k^1\left(\underset{\beta }{}f_\beta ^N\frac{g_\beta ^N}{q_k}+f_\beta ^T\frac{g_\beta ^T}{q_k}\right)\hfill \\ \hfill \ddot{g}_\alpha ^T=& \frac{g_\alpha ^T}{q_k}\left(\widehat{M}_k^1Q_k\right)+\frac{^2g_\alpha ^T}{q_kq_l}\dot{q}_k\dot{q}_l+\frac{g_\alpha ^T}{q_k}\widehat{M}_k^1\left(\underset{\beta }{}f_\beta ^N\frac{g_\beta ^N}{q_k}+f_\beta ^T\frac{g_\beta ^T}{q_k}\right).\hfill \end{array}$$
(24)
Renaming again $`\ddot{g}_\alpha ^{\{N,T\}}`$ by $`a_\alpha ^{\{N,T\}}`$ and using the abbreviations

$$\begin{array}{cc}\hfill b_\alpha ^N=& \frac{g_\alpha ^N}{q_k}\left(\widehat{M}_k^1Q_k\right)+\frac{^2g_\alpha ^N}{q_kq_l}\dot{q}_k\dot{q}_l\hfill \\ \hfill b_\alpha ^T=& \frac{g_\alpha ^T}{q_k}\left(\widehat{M}_k^1Q_k\right)+\frac{^2g_\alpha ^T}{q_kq_l}\dot{q}_k\dot{q}_l\hfill \end{array}$$
(25)
and

| $`A_{\alpha \beta }^{NN}`$ | $`={\displaystyle \frac{g_\alpha ^N}{q_k}}\widehat{M}_k^1{\displaystyle \frac{g_\beta ^N}{q_k}}`$ | $`A_{\alpha \beta }^{NT}`$ | $`={\displaystyle \frac{g_\alpha ^N}{q_k}}\widehat{M}_k^1{\displaystyle \frac{g_\beta ^T}{q_k}}`$ |
| --- | --- | --- | --- |
| $`A_{\alpha \beta }^{TN}`$ | $`={\displaystyle \frac{g_\alpha ^T}{q_k}}\widehat{M}_k^1{\displaystyle \frac{g_\beta ^N}{q_k}}`$ | $`A_{\alpha \beta }^{TT}`$ | $`={\displaystyle \frac{g_\alpha ^T}{q_k}}\widehat{M}_k^1{\displaystyle \frac{g_\beta ^T}{q_k}}`$ |

(26)
we can write the modified geometry equations

$$\begin{array}{cc}\hfill a_\alpha ^N=& b_\alpha ^N+\underset{\beta }{}\left(A_{\alpha \beta }^{NN}f_\beta ^N+A_{\alpha \beta }^{NT}f_\beta ^T\right)\hfill \\ \hfill a_\alpha ^T=& b_\alpha ^T+\underset{\beta }{}\left(A_{\alpha \beta }^{TN}f_\beta ^N+A_{\alpha \beta }^{TT}f_\beta ^T\right).\hfill \end{array}$$
(27)
The full set of consistency conditions then reads

$$\begin{array}{cc}\hfill a_\alpha ^N& 0\hfill \\ \hfill f_\alpha ^N& 0\hfill \\ \hfill a_\alpha ^Nf_\alpha ^N=& 0\hfill \\ \hfill \left|f_\alpha ^T\right|& \mu f_\alpha ^N\hfill \\ \hfill a_\alpha ^T\left(\left|f_\alpha ^T\right|\mu f_\alpha ^N\right)=& 0.\hfill \end{array}$$
(28)
Equations (27) together with the conditions (28) can be solved with a modified Dantzig’s algorithm  which will be discussed in the next section. Note that some of the tangential forces may be directly determined by the respective normal forces, i.e., when sliding at this contact occurs. The consistency conditions for these forces have to be fulfilled, nevertheless.

## 6 Dantzig’s Algorithm

Contacts can be classified into breaking contacts ($`a^N=0`$ thus $`f^N=f^T=0`$), permanent static contacts, ($`a^N=0`$ and $`a^T=0`$) or permanent sliding contacts ($`a^N=0`$ but $`a^T0`$). If we knew à priory into which category each contact belongs the contact forces could be determined by solving an inhomogeneous system of linear equations which consists of all equations for which either $`a_\alpha ^N=0`$ or $`a_\alpha ^T=0`$, with the corresponding $`f_\alpha ^N`$ and $`f_\alpha ^T`$ as variables. All remaining normal forces are zero, the remaining tangential forces assume their maximum values. Unfortunately the contact classification is only known if we know the contact forces as well.

We apply Dantzig’s Algorithm to determine the forces together along with the corresponding contact classification. It starts with considering a certain contact, disregarding all others, i.e., their contact forces are set to zero. After having found a solution for this contact its classification is also known. The algorithm proceeds then with the next contact. Again the contact forces are determined preserving the consistency of all contacts considered before. In this process the contact classifications of the already consistent contacts may be changed if necessary. The process is repeated until the last contact has been classified.

All contacts are assigned to one of four lists:

* List $`NC`$ of breaking contacts
* List $`C_F`$ of permanent static contacts
* List $`C^\pm `$ of permanent sliding contacts. In list $`C^+`$ are all contacts where $`f^T=\mu f^N`$, in $`C^{}`$ all contacts where $`f^T=\mu f^N`$.

All contacts in the above lists are considered to be consistent, i.e., they satisfy the consistency conditions. The classification is done successively for all contacts $`\alpha `$ by:

1. Check if the normal force $`f_\alpha ^N=0`$ satisfies the consistency conditions $`a_\alpha ^N0`$. If this is the case the contact is consistent and belongs to $`NC`$.
2. If $`a_\alpha ^N<0`$ we have to increase the normal contact force $`f_\alpha ^N`$ to obtain a non-negative normal acceleration. However, increasing the normal force $`f_\alpha ^N`$ will change the contact accelerations of the already classified contacts as well. Since for the permanent contacts $`\beta `$ we need $`a_\beta ^N=0`$ (and for the static contacts also $`a_\beta ^T=0`$) this would invalidate the classification of these contacts. To preserve the consistency of the classified contacts we will have to change the contact forces of the the permanent contacts as well. We now have to determine how much we have to change these contact forces for a given increase $`s`$ of the new force $`f_\alpha ^N`$ to keep $`a_\beta ^N=0`$ and, if necessary $`a_\beta ^T=0`$. To calculate the necessary changes we formulate a reduced set of geometry equations:

$$\begin{array}{cc}\hfill 0=a_\beta ^N=& b_\beta ^N+A_{\beta \gamma }^{NN,\mathrm{red}}f_\gamma ^N+A_{\beta \gamma }^{NT,\mathrm{red}}f_\gamma ^T+A_{\beta \alpha }^{NN}f_\alpha ^N\hfill \\ \hfill 0=a_\beta ^T=& b_\beta ^T+A_{\beta \gamma }^{TN,\mathrm{red}}f_\gamma ^N+A_{\beta \gamma }^{TT,\mathrm{red}}f_\gamma ^T+A_{\beta \alpha }^{TN}f_\alpha ^N.\hfill \end{array}$$
(29)
The reduced set of geometry equations can be obtained from the original geometry equations (27) disregarding all breaking contacts and replacing $`f_\beta ^T=\pm \mu f_\beta ^N`$ for all contacts $`\beta `$ from $`C^\pm `$ (sliding contacts). Since only permanent contacts remain all contact accelerations in the reduced geometry equations are thus zero. If we now change the new contact force $`f_\alpha ^Nf_\alpha ^N+s`$ we have to vary the previously known contact forces $`f_\beta ^{\{N,T\}}`$ in the reduced geometry equations by an unknown amount $`\mathrm{\Delta }f_\beta ^{\{N,T\}}`$ in order to keep the contact accelerations at their required value of zero:

$$\begin{array}{cc}\hfill 0=& b_\beta ^N+A_{\beta \gamma }^{NN,\mathrm{red}}\left(f_\gamma ^N+\mathrm{\Delta }f_\gamma ^N\right)+A_{\beta \gamma }^{NT,\mathrm{red}}\left(f_\gamma ^T+\mathrm{\Delta }f_\gamma ^T\right)\hfill \\ & +A_{\beta \alpha }^{NN}\left(f_\alpha ^N+s\right)\hfill \\ \hfill 0=& b_\beta ^T+A_{\beta \gamma }^{TN,\mathrm{red}}\left(f_\gamma ^N+\mathrm{\Delta }f_\gamma ^N\right)+A_{\beta \gamma }^{TT,\mathrm{red}}\left(f_\gamma ^T+\mathrm{\Delta }f_\gamma ^T\right)\hfill \\ & +A_{\beta \alpha }^{TN}\left(f_\alpha ^N+s\right).\hfill \end{array}$$
(30)
Inserting Eq. (29) we find

$$\begin{array}{cc}\hfill 0=& A_{\beta \gamma }^{NN,\mathrm{red}}\mathrm{\Delta }f_\gamma ^N+A_{\beta \gamma }^{NT,\mathrm{red}}\mathrm{\Delta }f_\gamma ^T+A_{\beta \alpha }^{NN}s\hfill \\ \hfill 0=& A_{\beta \gamma }^{TN,\mathrm{red}}\mathrm{\Delta }f_\gamma ^N+A_{\beta \gamma }^{TT,\mathrm{red}}\mathrm{\Delta }f_\gamma ^T+A_{\beta \alpha }^{TN}s.\hfill \end{array}$$
(31)
This is linear system of equations for the unknown $`\mathrm{\Delta }f_\gamma ^{\{N,T\}}`$ with the step size $`s`$ being a parameter. The necessary variations $`\mathrm{\Delta }f_\gamma ^{\{N,T\}}`$ are proportional to the step size $`s`$, hence the solution is of the form

$$\mathrm{\Delta }f_\gamma ^{\{N,T\}}=F_\gamma ^{\{N,T\}}s,$$
(32)
where $`F_\gamma ^{\{N,T\}}`$ is the necessary variation if $`s=1`$. By inserting the changed values back into the original geometry equations we now know as well by how much the accelerations of the breaking and sliding contacts change.
3. We increase the value of the new force $`f_\alpha ^N`$ until either the acceleration $`a_\alpha ^N=0`$ (the contact is now consistent) or until the classification of any already consistent contact changes. The classification changes if either

1. a normal acceleration $`a_\beta ^N>0`$ becomes zero: The contact becomes permanent and has to be moved from $`NC`$ to either $`C_F`$ or $`C^\pm `$ according to its present value of the corresponding tangential acceleration.
    2. a normal force $`f_\beta ^N>0`$ (permanent contact) becomes zero: the contact is now a breaking contact and has to be moved from $`C_F`$ or $`C^\pm `$ to $`NC`$.
    3. a tangential acceleration $`a_\beta ^T0`$ becomes zero and the corresponding tangential velocity is zero: the contact now is static and has to be moved from $`C^\pm `$ to $`C_F`$.
    4. a friction force previously of smaller magnitude than its allowed maximum reaches the maximum value $`\pm \mu f_\beta ^N`$: the contact becomes sliding and has to be moved from $`C_F`$ to $`C^\pm `$.
4. if we have not yet found consistent values for $`a_\alpha ^N`$ and $`f_\alpha ^N`$ we have to proceed with step 2.

If the contact $`\alpha `$ is permanent we have to consider the tangential component $`a_\alpha ^T`$ of contact $`\alpha `$ too. The procedure is very similar to the calculation of $`a_\alpha ^N`$. The only difference is that if $`a_\alpha ^T>0`$ we have to decrease the friction force until it assumes its negative maximum value and if $`a_\alpha ^T<0`$ we have to increase the friction until it adopts its positive maximum value.

## 7 Collisions

In the framework of Rigid Body Dynamics collisions occur if two contacting particles have a negative normal relative velocity at their contact point. We can easily convince ourselves that no finite contact force can prevent a deformation of the particles. No matter how large the force is, it will always take a finite, however small, time to stop the approaching particles, hence they will deform each other. Thus, to prevent deformation of the particles an infinite repulsive force of infinitesimal duration is necessary. It turns out that the total momentum transfer $`\mathrm{\Delta }p`$ between the two particles is finite222Note that the duration $`t_c`$ of the collision is only formally the parameter of the limit in the equation above. This limit is actually achieved by starting with deformable particles and increasing their stiffness to infinity. In this limiting process the duration of collision approaches zero.:

$$\mathrm{\Delta }p=\underset{t_c0}{lim}\underset{0}{\overset{t_c}{}}f𝑑t=\mathrm{finite}.$$
(33)
Resolving multi-particle collisions it turns out that the resulting state of the particles after the collision is not unique. This is due to the infinite stiffness of the particles, or equivalently, the infinite speed of sound in the particle material. In the limit of infinite speed of sound all information on the exact collision mechanism, e.g., the sequence of individual pair collisions, is lost, since any of the collisions is of vanishing duration. Colliding particles of finite stiffness do not exhibit this feature, since here all processes take finite time.

According to these arguments the necessary information on the detailed collision mechanism is not accessible. The following set of assumptions turned out to yield realistic results, although they cannot be uniquely derived:

1. All individual pair collisions occur at once.
2. The transfer of momentum at contact points is finite. There is no momentum transfer which corresponds to attractive forces.
3. The relative velocity after a collision can never be smaller than $`ϵv`$, where $`ϵ`$ is the coefficient of restitution and $`v`$ is the impact velocity at that contact ($`v<0`$ !).
4. If the velocity after the collision is strictly larger than $`ϵv`$ there is no momentum transfer at this contact.

The first two assumptions have been introduced already. The remaining two assumptions deserve further discussion. Two-particle collisions can be described by means of the coefficient of restitution which relates the precollisional relative velocity $`v`$ and the final velocity $`v^{}`$ after the collision:

$$v^{}=ϵv.$$
(34)
In the case of multi-particle collisions, however, two particles which initially rest relatively to each other may separate after a collision. Therefore, the final velocity may indeed be larger than the value $`ϵv`$. Contrary it may no happen that two particle which collide with a finite impact velocity are at rest relative to each other afterwards. Therefore, $`v^{}ϵv`$.

The fourth assumption simply means, that if two particles separate from each other with higher velocity than $`ϵv`$ their aftercollisional velocity may not be increased further by an additional momentum transfer.

For convenience we define the excess velocity $`\mathrm{\Delta }v`$

$$\mathrm{\Delta }v=v^{}+ϵv,$$
(35)
which is zero if $`v^{}=ϵv`$. The source of any velocity change is a momentum transfer between the colliding particles. We can relate the excess velocities at the contacts with the momentum transfers by means of the collisional geometry equation, which is derived in a similar way as the geometry equation of the force algorithm:

$$\mathrm{\Delta }v_\alpha =\left(1+ϵ\right)v_\alpha +A_{\alpha \beta }^{\mathrm{col}}\mathrm{\Delta }p_\beta .$$
(36)
Note that the $`v_\alpha `$ are relative velocities of the contact points of two contacting particles, but not the velocities of the particles themselves. In mathematical terms the above discussed assumptions read

$$\begin{array}{cc}\hfill \mathrm{\Delta }v_\alpha & 0\hfill \\ \hfill \mathrm{\Delta }p_\alpha & 0\hfill \\ \hfill \mathrm{\Delta }v_\alpha \mathrm{\Delta }p_\alpha & =0.\hfill \end{array}$$
(37)
The first condition prevents the final velocity from being smaller than $`ϵv`$. The second condition excludes attractive interaction between particles. Finally, the third condition means that there may be a finite momentum transfer only if $`\mathrm{\Delta }v=0`$. These conditions together with the geometry equation (36) form a Linear Complementarity Problem which is already familiar from Sec. 6 and can be solved by Dantzig’s Algorithm.

## 8 Resolution of static indeterminacy

If the number of contacts in the system is too large, the contact forces cannot be uniquely determined by the force algorithm. This occurs, if the number of free variables in the system, i.e. the number of contact forces, is larger than the number of mechanical degrees of freedom, $`3N`$ in 2d or $`6N`$ in 3d, with $`N`$ being the number of particles. However, the total forces and torques acting on the particles and, hence, their trajectories are unique. This drawback restricts the applicability of Rigid Body Dynamics for the simulation of railway ballast, since this system is of mainly static nature, the exact knowledge of the contact forces is crucial for understanding its behavior.

So far we have considered the contact forces as independent of each other. This assumption is the reason for the force indeterminacy. In realistic systems, however, the forces are not independent (Fig. 3). If we apply an external force directed to the right on the central particle, we increase the contact force with the particle to its right while at the same time decreasing the contact force with the particle to its left. In this example both contact forces in reality depend on a single parameter, which is the applied external force.

We can mimic this behavior by introducing small displacements of the particles. Each particle has a set of macroscopic coordinates $`\stackrel{}{r}`$ and $`\stackrel{}{\varphi }`$ and a set of microscopic coordinates $`\delta \stackrel{}{r}`$ and $`\delta \stackrel{}{\varphi }`$. In analogy to the vector $`q`$ of all (macroscopic) coordinates of all particles, we define the vector $`\delta q`$ of all microscopic (infinitesimal) coordinates. The contact network and the kinematic state of the particles is determined solely by the macroscopic coordinates, hence they can be considered as the actual coordinates of the particles. The microscopic coordinates lead to a deformation of contacting particles. By definition they are of infinitesimal size, which allows us to restrict ourselves to a linear approximation in $`\delta q`$ for the computation of the deformations. The vector $`\xi `$ of all deformations in the system

$$\xi =\widehat{D}\delta q$$
(38)
is defined by the deformation matrix $`\widehat{D}`$ and the microscopic coordinates. The dependence of $`\widehat{D}`$ on the geometric properties of the systems is straightforward but lengthy, therefore we will not give explicit expressions here (for the full derivation see ). We define a force law to relate the deformations at the contact points with the contact forces:

$$f=f(\xi ).$$
(39)
Now the contact forces are functions of the displacements $`\delta q`$. To calculate the forces we use the geometry equation (27) together with the consistency conditions (28). This set of equations is to be solved for the microscopic displacements $`\delta q`$. Hence, there are as many variables as degrees of freedom, i.e., the system is unique.

To determine the microscopic coordinates we use an overdamped relaxation method. We start with a set of inconsistent coordinates $`\delta q`$. Inconsistent means that the consistency conditions for the resulting forces and accelerations are not fulfilled. Now we let the system relax. This way we find new microscopic coordinates such that $`\xi ^{}=\xi ha`$, with $`h`$ being the step size and $`a`$ the vector of all contact accelerations. If the contact acceleration is negative the displacement will be larger, yielding a larger force to stop the approaching motion of the contacting particles. Since the adjustment step for the displacements is proportional to the acceleration the sequence of microscopic coordinates of the particles can be understood as overdamped motion.

The adjusted microscopic coordinates are the solution of the linear system of equations

$$\xi ha=\widehat{D}\delta q^{}$$
(40)
or, equivalently

$$\widehat{D}\left(\delta q^{}\delta q\right)=ha.$$
(41)
In cases where there are less degrees of freedom than contact accelerations the system of equations is overdetermined. There may be vectors $`ha`$ which are not representable by any vector $`\delta q^{}\delta q`$. In this case we have to project the vector $`ha`$ into the image space of the operator $`\widehat{D}`$ before solving the system of equations.

We repeat the adjustment of microscopic coordinates until the consistency conditions are met. To further improve the speed of this method we can save the microscopic coordinates which yielded a consistent system in the previous time step. If the system did not change too much this set of microscopic coordinates is very close to the new solution and we need only few iteration steps to arrive at the new solution.

This method combines the advantage of Rigid Body Dynamics, namely the ability to simulate very stiff particles, with the advantage of Molecular Dynamics, namely uniquely defined contact forces. An additional advantage of the method of small displacements is that we can now easily simulate certain degradation mechanisms. If we, for example, want to simulate the effects of abrasion of edges of the particles we can do this by gradually changing the force law (39), which describes changing edge properties of the particles.

## 9 Step size control

The integration scheme used in our implementation is a Runge-Kutta method of fourth order. During one time step there are four force computations necessary. Since we use discrete time steps we are frequently faced with the problem that after a given time step some particles do deform each other. In this case, obviously, the chosen time step was too large. This problem could be solved by predicting the time of next contact (collision) from the present state of the particles. Since the particles are subject to forces which vary in time this prediction cannot be accurate, as we approach a collision we would have to update it repeatedly. Since this prediction method is quite complicated we chose a simpler method. We advance the system by a certain time step. When determining the contacts in the next time step (see Sec. 4) we check for deformations of the particles. If any deformations occur we restore the state of the system before this time step and take a time step of half its original value. We then repeat this process, i.e. advance the system by the new time step and check for deformations. If we again encounter deformations we again divide the time step by two and repeat the computation until a state without deformations is reached. The new state of the system is accepted and the other steps of the algorithm are performed. This procedure ensures that there will be no particle deformations in the system at the beginning of any accepted time step.

Since this method can only decrease the time step, however, we need a procedure to increase the time step again333There is, of course, an upper limit to the time step which is dictated by the accuracy of the chosen integration scheme. in order to avoid permanently slowing down the simulation. We can always increase the time step if there are no collisions in the systems for a certain time. Hence, if we successfully performed a number of time steps without encountering collisions we can increase the time step again by a factor of two. We cannot increase the time step after only one successful time step since this may lead to frequently alternating increase and decrease steps. A requirement of three successful time steps before increasing the time step has shown to yield good computational efficiency.

With this choice of a step size control method we have completed the discussion of the simulation algorithm.

## 10 Conclusions

Although Molecular Dynamics methods have recently been applied successfully to the simulation of the dynamics of railway ballast  we have doubts that this technique is suitable to the simulation of almost rigid, sharply edged particles such as railway ballast, in particular if the system behavior is governed mainly by static properties of the system. Moreover, it seems to be unsuitable for the simulation of long time effects such as densification and wear of ballast.

The Rigid Body Dynamics is a method which is intended to describe the motion of systems of very stiff particles. In a natural way it avoids the problems of Molecular Dynamics simulations which have been discussed in detail in Sec. 2. These problems are unavoidable within the concept of Molecular Dynamics. Therefore we believe that Rigid Body Dynamics is much better suited to the simulation of railway ballast than Molecular Dynamics. From the presentation of the algorithm we have seen that a time step in Rigid Body Dynamics is a lot more complicated than a time step in Molecular Dynamics. Hence the according implementation requires by far more computing time than an implementation of a Molecular Dynamics code. Fortunately this disadvantage is offset by the fact that in Rigid Body Dynamics we can chose a larger time step than in Molecular Dynamics. Whereas in Molecular Dynamics the time step is determined by the critical deformation (the length on which the force changes by one unit) divided by the characteristic velocity in the system, which usually yields a very small time step (e.g. $`10^7`$ sec), the time step in Rigid Body Dynamics is determined by the characteristic time in which the forces in the system change significantly. The latter quantity is usually much larger ($`10^3`$ sec). Thus, the computing time requirements for both methods are comparable when simulating mainly static problems as it is the case for railway ballast.

An additional advantage of Rigid Body Dynamics is that fracture of particle can be introduced in a very natural way. When dealing with deformable particles in MD simulations the fracture of particles can lead to numerical difficulties. After the fracture of a particle, i.e., when a large particle is replaced by two or more fragments, there will be a finite time in which the system has to readjust and to reach a new stable state. There are situations when the new stable state deviates significantly from the realistic (experimental) situation. This problem is not encountered in Rigid Body Dynamics. After a fragmentation (which is algorithmically done by constructing two or more new polygons from the previous particle simply by introducing another edge into the old particle and splitting along this line) the system remains always stable.

It is our firm believe that Rigid Body Dynamics provides a very promising alternative to Molecular Dynamics for the simulation of railway ballast.


format: text

text:
# 1 Introduction

## 1 Introduction

Despite the impressive progress of these last years, 2d quantum field theories (QFT) still provide many interesting open problems. One of these is the finite temperature behaviour of correlators. Besides its experimental and theoretical relevance, this issue is particularly important when these theories are studied by means of numerical simulations. In fact choosing a finite temperature setting corresponds, in the euclidean formulation of the theory, to compactify the (imaginary) time direction on a circle whose circumference $`R`$ coincides with the inverse temperature. Thus the finite temperature corrections become, in this framework, finite size corrections in the cylinder geometry which, as it is well known, play a crucial role in the extraction of continuum limit expectation values from numerical simulations.

While for a generic 2d QFT finding a functional form for the finite temperature correction seems a hopeless task, recently two interesting proposals by Delfino  and LeClair and Mussardo  (based on a previous work by LeClair and collaborators ) appeared in the literature to address this problem in the simpler case of integrable models. In particular, in , the authors studied the finite temperature behaviour of one-point correlation functions. General arguments suggest that, outside the critical point, any one-point function evaluated on a finite size lattice $`\mathrm{\Phi }_R`$ should approach its infinite volume limit$`\mathrm{\Phi }_l_{R=\mathrm{}}`$ with an exponential decay of the type

$$\mathrm{\Phi }_l_R\mathrm{\Phi }_l_{R=\mathrm{}}+Ce^{\frac{R}{\xi }}+\mathrm{}$$
(1)
$`\xi `$ being the correlation length (i.e. the inverse of the lowest mass of the theory). However with these considerations nothing can be said on the constant $`C`$ and the higher order terms in the above equation.

The main achievement of  was to show that the constant $`C`$ is indeed universal and to predict its value. At the same time they were able to give an explicit expression for higher order corrections.

The two proposals  and , have different theoretical starting points and indeed the predictions for these higher order corrections turn out to be different in the two approaches (for a critical comparison and a discussion of these differences see ). However, if one looks at the first few orders in a low temperature expansion (those which do not involve multi-particle form factors) the results coincide.

The aim of this paper is to test these proposals in a particular integrable model: the 2d Ising model in a magnetic field. This model has recently been the subject of several theoretical and numerical studies -. In particular, in the numerical analyses (performed using transfer matrix techniques  or montecarlo simulations ) a great attention was devoted to the analysis of finite size correction. However, despite this careful treatment these correction were the major source of uncertainty in the final results. Our goal in this paper is thus twofold: first we shall use the numerical results for finite size lattices to test the proposals of . Second, we shall see if, assuming the functional form for the finite size corrections proposed in  we can improve the precision of the continuum limit results for the one -point correlation functions.

This paper is organized as follows. In sect.2 we shall briefly discuss the model, the observables and the proposal of . Sect.3 will be devoted to a discussion of our numerical test and to a comparison with . Finally sect.4 will be devoted to some concluding remarks.

## 2 Ising model in a magnetic field

The Ising model in a magnetic field is defined by the partition function

$$Z=\underset{\sigma _i=\pm 1}{}e^{\beta (_{n,m}\sigma _n\sigma _m+H_n\sigma _n)}$$
(2)
where the field variable $`\sigma _n`$ takes the values $`\{\pm 1\}`$; $`n(n_0,n_1)`$ labels the sites of a square lattice size $`L_0`$ and $`L_1`$ in the two directions and lattice spacing $`a`$111Since the lattice spacing will play no role in the following we shall set $`a=1`$ in the rest of the paper.. $`n,m`$ denotes nearest neighbour sites on the lattice. In the following we shall treat asymmetrically the two directions. We shall denote $`n_0`$ as the compactified “time” coordinate and $`n_1`$ as the space one. The number of sites of the lattice will be denoted by $`NL_0L_1`$. The lattice extent in the transverse direction will be denoted as $`RL_0`$. This length will play a major role in the following. In fact our goal will be to describe the $`R`$ dependence of the expectation values of the spin and energy operators of the theory.

In order to select only the magnetic perturbation, the coupling $`\beta `$ must be fixed to its critical value

$$\beta =\beta _c=\frac{1}{2}\mathrm{log}(\sqrt{2}+1)=0.4406868\mathrm{}$$

by defining $`h_l=\beta _cH`$ we end up with

$$Z(h_l)=\underset{\sigma _i=\pm 1}{}e^{\beta _c_{n,m}\sigma _m\sigma _m+h_l_n\sigma _n}.$$
(3)
In the continuum limit the model is described by the action:

$$𝒜=𝒜_0+hd^2x\sigma (x).$$
(4)
where $`\sigma (x)`$ is the perturbing operator and the perturbing field is the magnetic field $`h`$.

In this paper we want to use results obtained in the continuum theory to describe the finite size scaling behaviour of the lattice model. In general this would require a precise definition of the relations between lattice and continuum quantities like for instance that between $`h_l`$ which appears in eq.(3) and the perturbing field $`h`$. However, as we shall see below, we shall only be interested in adimensional ratios in which all these normalizations cancel out. Thus we shall neglect this problem in the following. The interested reader can find a careful treatment of this issue in .

### 2.1 Lattice operators

As a first step of our analysis let us define the lattice analogous of the spin and energy operators of the continuum theory222Strictly speaking these lattice operators do not correspond exactly to the continuum ones but are instead linear combinations of all possible relevant and irrelevant operators of the continuum theory with compatible symmetry properties with respect to the $`Z_2`$ symmetry of the model (odd for the spin operator and even for the energy one). However near the critical point this linear combination is dominated by the corresponding relevant operator ($`\sigma `$ for $`\sigma _l`$ and $`ϵ`$ for $`ϵ_l`$) and the only remaining freedom will be a conversion constant relating the continuum and lattice versions of the two operators. As mentioned above we shall neglect in the following these normalization constants.. The simplest choices for these lattice analogous are

* Spin operator

$$\sigma _l(x)\sigma _x$$
(5)
i.e. the operator which associates to each site of the lattice the value of the spin in that site.
* Energy operator

$$ϵ_l(x)\frac{1}{4}\sigma _x\left(\underset{yn.n.x}{}\sigma _y\right)ϵ_b$$
(6)
where the sum runs over the four nearest neighbour sites $`y`$ of $`x`$. $`ϵ_b`$ represents a constant “bulk” term which we shall discuss below 333As a matter of fact, for technical reasons, in our Transfer Matrix calculations we evaluated only the “time-like” part of the action i.e. $`_{n_0,n_1}\sigma _{(n_0,n_1)}\sigma _{(n_0+1,n_1)}`$. While, obviously, this choice makes no difference in the thermodynamic limit, it might have some effect for finite values of $`R`$. We shall further discuss this point in sect.3.1 below..

The index $`l`$ indicates that these are the lattice discretizations of the continuous operators. We shall denote in the following the normalized sum over all the sites of these operators simply as

$$\sigma _l\frac{1}{N}\underset{x}{}\sigma _l(x)ϵ_l\frac{1}{N}\underset{x}{}ϵ_l(x).$$
(7)
In the following we shall be interested in the expectation value of these lattice operators. More precisely we shall be interested in:

* Magnetization

The magnetization per site $`M(h_l)`$ defined as

$$M(h_l)\frac{1}{N}\frac{}{h_l}(\mathrm{log}Z(h_l))|_{\beta =\beta _c}=\frac{1}{N}\underset{i}{}\sigma _i.$$
(8)
which implies

$$M(h_l)=\sigma _l.$$
(9)
* Internal Energy

The internal energy density $`\widehat{E}(h_l)`$ defined as

$$\widehat{E}(h_l)\frac{1}{2N}\underset{n,m}{}\sigma _n\sigma _m.$$
(10)
It is important to stress that in this case one also has to take into account a bulk analytic contribution (as it happens also for the free energy itself) $`E_b(h_l)`$ which is an even function of $`h_l`$. Let us define $`ϵ_bE_b(0)`$. The value of $`E_b(0)`$ can be easily evaluated (for instance by using Kramers-Wannier duality) to be $`ϵ_b=\frac{1}{\sqrt{2}}`$. Let us define $`E(h_l)\widehat{E}(h_l)ϵ_b`$, we have

$$E(h_l)=\frac{1}{2N}\underset{n,m}{}\sigma _n\sigma _m\frac{1}{\sqrt{2}}.$$
(11)
Hence we have

$$E(h_l)=ϵ_l.$$
(12)

In the following we shall in particular be interested in the ratios

$$\frac{\sigma _l_R}{\sigma _l_{R=\mathrm{}}}\frac{ϵ_l_R}{ϵ_l_{R=\mathrm{}}}$$
(13)
$`\mathrm{\Phi }_l_R`$ being the mean value on a lattice with transverse extent $`R`$ of the lattice operator $`\mathrm{\Phi }_l`$. Since in the ratio all the normalization constants cancel out, we can identify the two ratios with the analogous ones evaluated in the continuum theory, with the continuum operators.

$$\frac{\mathrm{\Phi }_l_R}{\mathrm{\Phi }_l_{R=\mathrm{}}}=\frac{\mathrm{\Phi }_R}{\mathrm{\Phi }_{R=\mathrm{}}}.$$
(14)
Notice that to perform this identification it is mandatory to eliminate the bulk correction as we did in eq.(11)

### 2.2 Critical behaviour

The critical behaviour of magnetization and internal energy can be easily obtained by means of standard renormalization group methods. One finds

$$M|h|^{d/y_h1}$$
(15)
$$E|h|^{(dy_t)/y_h}.$$
(16)
where $`y_h`$, $`y_t`$ are the are the magnetic and thermal RG-exponents respectively.

From the exact solution of the Ising model at the critical point we know that $`y_h=\frac{15}{8}`$ and $`y_t=1`$. Inserting these values in the above expressions we find

$$M|h|^{\frac{1}{15}}$$
(17)
$$E|h|^{\frac{8}{15}}.$$
(18)
From the conformal field theory (CFT) description of the Ising model at the critical point it is possible to construct also the $`h`$ dependence of the following terms in the two scaling functions (17) and (18). A detailed analysis can be found in . We shall make use of this result in the following.

### 2.3 S-matrix results

In 1989 A. Zamolodchikov  suggested that the scaling limit of the Ising Model in a magnetic field could be described by a a scattering theory which contains eight different species of self-conjugated particles $`A_a`$, $`a=1,\mathrm{},8`$ with masses

$`m_2`$ $`=`$ $`2m_1\mathrm{cos}{\displaystyle \frac{\pi }{5}}=(1.6180339887..)m_1,`$
$`m_3`$ $`=`$ $`2m_1\mathrm{cos}{\displaystyle \frac{\pi }{30}}=(1.9890437907..)m_1,`$
$`m_4`$ $`=`$ $`2m_2\mathrm{cos}{\displaystyle \frac{7\pi }{30}}=(2.4048671724..)m_1,`$
$`m_5`$ $`=`$ $`2m_2\mathrm{cos}{\displaystyle \frac{2\pi }{15}}=(2.9562952015..)m_1,`$ (19)
$`m_6`$ $`=`$ $`2m_2\mathrm{cos}{\displaystyle \frac{\pi }{30}}=(3.2183404585..)m_1,`$
$`m_7`$ $`=`$ $`4m_2\mathrm{cos}{\displaystyle \frac{\pi }{5}}\mathrm{cos}{\displaystyle \frac{7\pi }{30}}=(3.8911568233..)m_1,`$
$`m_8`$ $`=`$ $`4m_2\mathrm{cos}{\displaystyle \frac{\pi }{5}}\mathrm{cos}{\displaystyle \frac{2\pi }{15}}=(4.7833861168..)m_1,`$
where $`m_1(h)`$ is the lowest mass of the theory. which coincides with the inverse of the (exponential) correlation length. After Zamolodchikov’s paper several other interesting results were obtained, ranging from the explicit values of various critical amplitudes  to the values of the overlap amplitudes of correlators . All these predictions have been tested with numerical simulations both in the 1d Ising quantum chain , in the dilute $`A_3`$ IRF (Interaction Round a Face) model  (which is another realization of the scaling Ising model in a magnetic field) and directly in the 2d lattice Ising model  and in all cases a full agreement between S-matrix predictions and numerical results was found.

### 2.4 Delfino and LeClair-Mussardo proposals

One of the most interesting features of the approach discussed in  is that, by using the form-factor technology the author was able to give a very explicit and compact expression for the first few orders of the finite size corrections both for $`\sigma `$ and $`ϵ`$ for the Ising model in a magnetic field. This result will be of great importance for our analysis.

According to  a generic one point function $`\mathrm{\Phi }_R`$ evaluated on a cylinder of transverse size $`R`$ approaches exponentially its asymptotic value $`\mathrm{\Phi }_{R=\mathrm{}}`$ with the following law.

$$\frac{\mathrm{\Phi }_R}{\mathrm{\Phi }_{R=\mathrm{}}}=1+\frac{1}{\pi }\underset{i=1}{\overset{3}{}}A_i^\mathrm{\Phi }K_0(m_iR)+O(e^{2m_1R}),$$
(20)
where $`m_1,m_2`$ and $`m_3`$ are the first three masses (the only ones below the lowest pair creation threshold) of the Zamolodchikov’s solution discussed above and $`K_0`$ denotes the zeroth order modified Bessel function.

The major result of  was to show that the $`A_i^\mathrm{\Phi }`$ constants are indeed universal and can be evaluated exactly in the framework of the S-matrix description of the model. Their values in the case of the two operators in which we are interested here are given in tab. 1. As mentioned in the introduction these results turn out to agree at this order with those obtained by LeClair and Mussardo in .

## 3 Numerical test

In order to test the above prediction we performed a numerical study of the finite size corrections by using the results of the transfer matrix analysis discussed in  to which we refer for further details on the algorithm and on the raw data. We only recall here some general informations on the experiment. We studied the Ising model for the 19 different values of the external magnetic field listed in tab.2. For each choice of $`h`$ we studied lattices with values of $`R`$ ranging from 10 to 21. For each values of $`h`$ and each lattice size we evaluated the mean value of the magnetization and internal energy (see however the footnote below eq.(6)). We show in tab.3 a typical sample of our data. An important ingredient of eq.(20) are the values of the first three masses of the theory. In principle these masses could be evaluated directly from the S-matrix solution, using the suitable normalization constants to match with the lattice discretization and then imposing the standard $`h`$ dependence dictated by the RG analysis. However it is important to stress that at the level of precision of our analysis, these three masses show rather large corrections to scaling for the values of $`h`$ that we study. To avoid large systematic errors in the analysis of the finite size correction it is mandatory to take into account these corrections (they are discussed in great detail in ). An alternative route, which turns out to be much simpler, is to evaluate these masses, when possible, directly from the transfer matrix calculation. For the present analysis we used this second route. We report for completeness these values in tab.2 and refer to  for a detailed discussion of this table.

### 3.1 Analysis of the data

We performed a two steps analysis

1\] First of all, for each value of $`h`$ we fitted the values of $`<\sigma >`$ and $`<ϵ>`$ as a function of $`R`$, according to eq.(20): The results are reported in tab.4. Here are some comments on the fits:

* We always kept the asymptotic value of $`<\sigma >`$ and $`<ϵ>`$ as a free parameter.
* The error of the input data was always of the order of $`10^{12}`$. We started by fitting all the values of $`R`$ and then eliminated them one by one, starting form the smallest ones, until we reached a reduced $`\chi ^2`$ equal or smaller than 1. Then we accepted the result of the fit and stopped the analysis. These acceptable $`\chi ^2`$’s could be achieved for $`h0.075`$ keeping only the first term in eq.(20). For $`h=0.066..,0.055..,0.044..`$ we had to introduce also the second mass. and finally for $`h=0.033..,0.022..`$ all the three terms were needed444This trend is due to the fact that we studied, for each value of $`h`$, the same range of values of $`R`$ while $`m_1`$ decreases as $`h`$ goes to $`0`$. This means that the argument of the first Bessel function in eq.(20) (i.e. $`m_1R`$) decreases as $`h0`$ thus allowing to detect, within the errors of our data also higher order terms of the equation.. The values of the masses which could not be directly evaluated via transfer matrix (the empty spaces in tab.2) were evaluated using the S-matrix results and keeping into account the proper corrections to scaling.
* For all the values of $`h`$ the asymptotic values that we obtained for $`<\sigma >`$ and $`<ϵ>`$ from these fits are compatible within the errors with those quoted in . In general the present estimates are more precise. This gain in precision becomes larger and larger as $`h`$ decreases. For instance for $`h=0.075`$ we found

$$E=0.122801011162(5)M=0.884110944919(2)$$
(21)
to be compared with the values reported in

$$E=0.1228010112(1)M=0.88411094491(1)$$
(22)
with an improvement of more than one order of magnitude in precision.

Let us stress that both this improvement and the fact that all our results agree within the errors (which means typically an agreement within 10 significative digits as in eq.s(21) and (22) ) with those of ref  are rather impressive evidences of the correctness of eq.(20)555The results of  were obtained with an iterative method (see sect.5.1 of ref. and sect.3.1.2 of  for a discussion) based on the general assumptions on the finite size behaviour of one point functions summarized in eq.(1). In particular the discussion of sect.3.1.2 of  shows that the agreement with the present results is not a case but is due to the fact that the iterative method actually mimics the exact behaviour of eq.(20). It would be nice to see if this iterative algorithm keeps its predicting power even when the theory is not integrable and no exact prediction can be deduced from a S-matrix analysis.. In fact, as it can be seen by looking at tab.3, the finite size corrections which eq. (20) is able to describe are much larger (five order of magnitude!) than the uncertainties of the data: the corrections are of the order of $`10^7`$ while the precision of the input data is $`10^{12}`$.
* For all the values of $`h`$ we obtained stable and reliable values for $`A_1^\mathrm{\Phi }`$. These values are reported in tab.4. The same was not true for $`A_2^\mathrm{\Phi }`$. This is not strange. In all the similar analyses performed in  it was impossible to have reliable estimates for the amplitudes of the subleading exponents due to the large systematic deviations induced by the uncertainties in the leading corrections.

Looking at the data of tab.4 one can see that there are rather large corrections to scaling. This is true in particular for $`A_1^ϵ`$ (for which it is possible that the rather large magnitude of the correction is a consequence of our asymmetric definition for the internal energy, see footnote below eq. (6)). In any case it is clear that in order to test the predicted values for $`A_1^\sigma `$ and $`A_1^ϵ`$ it is mandatory to discuss these deviations. To this end we performed a second level of analysis

2\] The aim of this second level of analysis is to use the data reported in tab.4 to reach the scaling limit values for $`A_1^\sigma `$ and $`A_1^ϵ`$. This problem is exactly the same that we had to face in  to extract the continuum limit values of the critical amplitudes and thus we use here the same method developed and discussed in sect.6 of  to which we refer the interested reader for further details.

We used as fitting functions

$$A_1^\sigma (h_l)=A_1^\sigma (1+b_{1,\sigma }|h_l|^{\frac{16}{15}}+b_{2,\sigma }|h_l|^{\frac{22}{15}})$$
(23)
$$A_1^ϵ(h_l)=A_1^ϵ(1+b_{1,ϵ}|h_l|^{\frac{8}{15}}+b_{2,ϵ}|h_l|^{\frac{16}{15}}).$$
(24)
Differently from the cases studied in  keeping only two terms in the scaling function is enough due to the large errors in the input data.

The final results are

$$8.090<|A_1^\sigma |<8.125$$
(25)
$$17.0<|A_1^ϵ|<18.1$$
(26)
which perfectly agree with the predictions of  reported in tab.1

Then (similarly to what we did in ) we tried to extract the dominant correction to scaling by assuming as fixed inputs the values for $`A_1^\mathrm{\Phi }`$ reported in tab.1 and performing again the fits. The results are:

$$0.56<b_{1,\sigma }<0.66$$
(27)
$$0.96<b_{1,ϵ}<0.82$$
(28)
As it also happens for the energy itself (see ) the leading correction to scaling in the energy sector is rather large and negative in sign.

### 3.2 An attempt to obtain $`A_2^\mathrm{\Phi }`$

It is impossible to extract the $`A_2^\mathrm{\Phi }`$ by directly fitting. Moreover we cannot improve the situation by using the known continuum limit values of $`<\mathrm{\Phi }>`$ and $`A_1^\mathrm{\Phi }`$ due to the large correction to scaling whose uncertainty is of the same order of magnitude of the terms that we would like to observe. The only possible way out is to construct a combination of the input data which exactly eliminates $`<\mathrm{\Phi }>`$ and $`A_1^\mathrm{\Phi }`$. This can be easily done by combining the values of $`\mathrm{\Phi }`$ measured at three different values of the lattice size $`R_1,R_2,R_3`$. The combination is the following:

$$A_2^\mathrm{\Phi }=\frac{\pi }{\mathrm{\Phi }_{R=\mathrm{}}}\frac{[\mathrm{\Phi }(R_1)\mathrm{\Phi }(R_2)]\mathrm{\Delta }_1(1,3)[\mathrm{\Phi }(R_1)\mathrm{\Phi }(R_3)]\mathrm{\Delta }_1(1,2)}{\mathrm{\Delta }_2(1,2)\mathrm{\Delta }_l(1,3)\mathrm{\Delta }_2(1,3)\mathrm{\Delta }_l(1,2)}$$
(29)
where

$$\mathrm{\Delta }_l(i,j)[K_0(m_lR_iK_0(m_lR_j)].$$
(30)
The result obtained using as input data those for $`h=0.022034\mathrm{}`$ for the magnetization are reported in tab.5. The quality of the result is rather good and is compatible with the prediction for this constant reported in tab.1. Notice however that the quality of the results becomes worse and worse as $`h`$ increases (they are already almost completely unstable at $`h=0.055085\mathrm{}`$) moreover no result can be obtained for $`ϵ`$ even at $`h=0.022034\mathrm{}`$. This can be easily understood. It due to the fact that the absolute value of $`<ϵ>`$ is much smaller than that of $`<\sigma >`$ and as a consequence the relative errors (which are those which generate the observed instabilities) are larger. We have checked that in our range of values of $`h`$ higher contributions in the expansion of eq.(20) essentially give no effect.

## 4 Concluding remarks

Let us briefly summarize our main results.

* Our results nicely agree with the functional form for the finite size behaviour of the one-point functions proposed by Delfino and LeClair-Mussardo. We are also able to give a good estimate of the first correction to scaling terms, both for $`\sigma `$ and for $`ϵ`$ and a rough but reliable estimate of the $`A_2^\sigma `$ for the magnetization.
* Unfortunately we are not able to discriminate between the two proposals. This would require precisions which are definitely outside the range of our transfer matrix methods.
* The asymptotic values for the energy and the magnetization obtained implementing eq.(20) are more precise that those that quoted in  of more than one order of magnitude. As the critical point is approached this improvement in precision increases.
* The agreement between the present results and those of , besides being a further strong evidence of the correctness of the proposal of  and , also points out the power of the algorithm proposed in  to deal with the finite size correction. This is an important result since this algorithm is not constrained to integrable theories

It would be interesting to extend this analysis to other integrable models which could offer examples in which it could be simpler to discriminate between the two proposals of  and  and to compare them with other existing approaches (say, for instance, ). It would also be important to extend the analysis to models which could be more easily accessible for comparison with experimental results (see for instance the results of  on the Haldane-Gapped spin chains). At the same time it would be interesting to extend the present analysis to two-point functions, for which very precise numerical results exist for the 2d Ising model in a magnetic field . This would be a rather important test since it has been recently claimed  that the extension to two-point functions of the proposals  and  could fail in the case of interacting theories.

### Acknowledgements

We thank G. Delfino and G. Mussardo for several useful discussions on the subject.


text:
# In: D. E. Wolf, P. Grassberger (Eds.) “Friction, Arching, Contact Dynamics” World Scientific (Singapore, 1997) p. 265-273 Force Distribution and Comminution in Ball Mills

## 1 Introduction

Milling is one of the most important techniques in mechanical engineering and much effort has been done to optimize the comminution efficiency of milling techniques with respect to energy, space and time consumption $`^\mathrm{?}`$.

Mainly because of their simple construction and application ball milling is a wide spread milling technology, particularly in mining. Ball mills are (usually large) cylindrical devices of length and diameter up to several meters which revolve along their axis. The material to be comminuted moves inside the cylinder. There are two different methods: autogeneous comminution where the material is pure and heterogeneous comminution where the material is mixed with heavy spheres from steel of typical diameter of several centimeters to increase efficiency. Throughout the paper we will deal with autogeneous comminution. Unfortunately the efficiency of ball mills is not very high and, therefore, engineers did much scientific work to increase their efficiency. There exists much experimental knowledge on the operation mechanisms of ball mills and on the comminution of granular matter in these mills. An overview can be found in ref. $`^\mathrm{?}`$.

There are still effects, however, which have not been understood yet. The present paper deals with the spatial distribution of stress in the material. At first glance one could assume that the largest stresses will be observed close to the surface of the material where the granular particles have the largest values of relative collision velocity. Hence one could assume that large part of particle comminution occurs close to the material surface. To increase efficiency of the mill one would have to tune the rotation velocity so that the average collision velocity becomes maximum.

Experimental investigations, however, show opposite results: Rothkegel and Rolf measured directly the spatial distribution of intensive impacts $`^\mathrm{?}`$ in a ball mill. In an experiment using a small almost two-dimensional mill they used instrumented balls which flash whenever the acting force at a particle contact is larger than a threshold. The statistics of the spatial distribution of the flashes led to a surprising result: the spatial distribution of intensive particle contacts which in realistic ball mills might lead to particle comminution has its maximum deep inside the material close to the bottom of the cylinder. This result was new and there was no satisfying explanation for the measured stress distribution yet.

In the following section we present the results of a two dimensional molecular dynamics simulation of a ball mill. It will be shown that the spatial stress distribution is closely related to the properties of force chains in the material. The occurrence of such force chains, and even of entire networks of such chains has been reported before, see e.g. refs. $`^\mathrm{?}`$-$`^\mathrm{?}`$.

## 2 Molecular dynamics model

The rotation velocity of ball mills is high enough to keep the granular particles in the continuous flow regime (for explanation of the regimes see $`^\mathrm{?}`$) and on top of the material grains are thrown through the air caused by intensive forcing due to rapid revolution of the mill. Therefore static friction does not play an important role for the dynamics of the material and, hence, we can apply a sphere model in the MD simulations. We just want to remark that in other cases when static behavior of granular material becomes essential for the dynamics one has to apply more complicated grain models of non-spherical shape. For detailed discussions of the differences in the dynamics of spherical and non-spherical objects see references $`^\mathrm{?}`$-$`^\mathrm{?}`$.

For our simulation we apply a spherical MD model for granular materials by Cundall and Strack $`^\mathrm{?}`$ and Haff and Werner $`^\mathrm{?}`$ which was used by many authors in simulations of rapidly moving granular material.

Two particles $`i`$ and $`j`$ with radii $`R_i`$ and $`R_j`$, positions $`\stackrel{}{r}_i`$, $`\stackrel{}{r}_j`$ and velocities $`\dot{\stackrel{}{r}}_i`$ and $`\dot{\stackrel{}{r}}_j`$ feel a force

$$\stackrel{}{F}_{ij}=F_{ij}^N\stackrel{}{n}+F_{ij}^T\stackrel{}{t}$$
(1)
only when contacting, i.e. if the condition $`\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|<R_i+R_j`$ holds. The force consist of a component $`F_{ij}^N`$ acting in normal direction $`\stackrel{}{n}`$

$$F_{ij}^N=Y\left(R_i+R_j\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|\right)m_{ij}^{\text{eff}}\gamma _N(\dot{\stackrel{}{r}}_i\dot{\stackrel{}{r}}_j)\stackrel{}{n}$$
(2)
and a component $`F_{ij}^T`$ acting in tangential direction $`\stackrel{}{t}=\left(\genfrac{}{}{0pt}{}{01}{10}\right)\stackrel{}{n}`$

$`F_{ij}^T`$ $`=`$ $`\text{sign}\left(v_{ij}^{\text{rel}}\right)\mathrm{min}(m_{ij}^{\text{eff}}\gamma _T\left|v_{ij}^{\text{rel}}\right|,\mu \left|F_{ij}^N\right|)`$ (3)
with (4)
$`v_{ij}^{\text{rel}}`$ $`=`$ $`(\dot{\stackrel{}{r}}_i\dot{\stackrel{}{r}}_j)\stackrel{}{t}+R_i\omega _i+R_j\omega _j`$ (5)
$`m_{ij}^{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{m_im_j}{m_i+m_j}}.`$ (6)
$`m_i`$ and $`\omega _i`$ are the mass and the rotation velocity of the $`i`$-th particle. $`Y=810^6gs^2`$ is the Young Modulus, $`\gamma _N=800s^1`$ and $`\gamma _T=3000s^1`$ are the damping coefficients in normal and tangential direction, and $`\mu =0.5`$ is the Coulomb friction coefficient. The parameters have been chosen to give realistic results in comparison of the simulation results with the behavior of a typical granular material.

Eq. (5) describes the relative velocity of the particle surfaces at the contact point. This velocity consists of the tangential part of the relative velocity of the particles and of a term which originates from rotation of the grains. $`m_{ij}^{\text{rel}}`$ stands for the effective mass.

The Coulomb friction law is taken into account by equation (3) saying that two grains slide on each other if the tangential force is larger than $`\mu `$ times the normal force. Otherwise they roll mainly. For the detailed discussion of this and other models for molecular dynamics simulation of granular material see $`^\mathrm{?}`$.

Based on the force (1) Newton’s equations of motion have been integrated using a Gear predictor-corrector scheme of fifth order, e.g. $`^\mathrm{?}`$. This method was used in many simulations of granular material and has been proven to be very numerically stable.

## 3 Results of the simulation

The model described in the previous section was used to simulate a two-dimensional ball mill of diameter $`D=8cm`$ filled with $`N=800`$ spherical grains. The radii of the grains are distributed in the interval $`[0.05,\mathrm{\hspace{0.17em}0.11}]cm`$. To study the influence of the rotation velocity of the cylinder we performed simulations for three different velocities $`\mathrm{\Omega }_I=2Hz=2/(2\pi )\text{ revolutions per second}`$, $`\mathrm{\Omega }_{II}=10Hz=10/(2\pi )\text{revolutions per second}`$, $`\mathrm{\Omega }_{III}=19Hz=19/(2\pi )\text{revolutions per second}`$.

At $`\mathrm{\Omega }=\mathrm{\Omega }_I=2Hz`$ the flow at the surface is already continuous and the surface is characterised by an inclined line (see fig. 1, top). For $`\mathrm{\Omega }=\mathrm{\Omega }_{II}=10Hz`$ the material surface is still almost compact and consists of two sections a steep and a flat one (fig. 1, middle). This velocity regime is experimentally known to have the best comminution efficiency $`^\mathrm{?}`$. Large rotation velocity $`\mathrm{\Omega }=\mathrm{\Omega }_{III}=19Hz`$ leads to free flight of the grains at the surface (fig. 1, top).

Fig. 1 shows snapshots of simulations with $`N=800`$ grains for the three mentioned rotation velocities $`\mathrm{\Omega }_I`$, $`\mathrm{\Omega }_{II}`$, $`\mathrm{\Omega }_{III}`$. The color codes for the local pressure $`P_i`$ acting on the $`i`$-th particle

$$P_i=\underset{j}{}F_{ij}^N.$$
(7)
The index $`j`$ runs over all neighbors of the $`i`$th particle. Grains feeling a pressure $`P_i>1000gcms^2`$ are drawn bold lined. Obviously most of the hardly stressed particles are located near the walls inside the bulk of material, which agrees well with the experimental findings by Rothkegel and Rolf $`^\mathrm{?}`$.

The hypothesis we want to discuss in the following is, that the self organized formation of force chains is responsible for the spatial stress distribution and, hence, the relevant physical reason for comminution processes in ball mills. Without the existence of such force chains, ball mills would not work or at least would work much less efficient. In the following we will discuss the properties of these force chains.

Force chains are defined by a set of simple conditions: Particles are called members of the same force chain if:

1. they feel a pressure (eq. (7)) of more than $`1000gcms^2`$,
2. each member particle of the chain touches at least one other member, and
3. for all members $`k`$ of the chain having two neighbors $`i`$ and $`j`$, which belong to the force chain as well, the condition

$$\left|\frac{\stackrel{}{r}_k\stackrel{}{r}_i}{\left|\stackrel{}{r}_k\stackrel{}{r}_i\right|}\frac{\stackrel{}{r}_k\stackrel{}{r}_j}{\left|\stackrel{}{r}_k\stackrel{}{r}_j\right|}\right|<0.85,$$
(8)
holds, i.e. these three particles lie almost on a line.

These three conditions can be checked using a computer algorithm. In fig. 1 neighboring particles have been connected by a black line provided that either of them belong to the same force chain. Obviously most of the hard stressed (and therefore bold drawn) particles are members of force chains. We conclude that an essential part of the static and dynamic pressure in the ball mill propagates along force chains. Particles which belong to a force chain feel pressure which are up to 100 times larger than the local average pressure

$$P_i^{\text{(av)}}=\frac{1}{K}\underset{k=1}{\overset{K}{}}P_k,$$
(9)
where the index $`k`$ runs over all $`K`$ neighbors of the particle $`i`$ which are within a certain neighborhood $`\left|\stackrel{}{r}_i\stackrel{}{r}_k\right|4\overline{R}`$ and which do not belong to any force chain.

The occurrence of force chains in granular systems is not new. They have been observed before in other granular systems in experiments $`^\mathrm{?}`$-$`^\mathrm{?}`$ and computer simulations $`^{\mathrm{?},\mathrm{?}}`$. It can be shown that one of the (at least) three phases of an idealized granular system is characterized by the occurrence of force chains $`^\mathrm{?}`$.

Fig. 2 (left) shows the frequency of force chains of lengths $`L`$ where $`L`$ is the number of particles which join the chain. The frequency decreases almost exponentially with increasing length. For the highest rotation velocity $`\mathrm{\Omega }_{III}`$ (fig. 2, bottom) the frequency distribution breaks down for longer chains, which is mostly due to less compact bulk of material and finite size of the system.

The right hand side of fig. 2 shows the correlation of the length of a force chain and the pressure which acts on the particles which belong to this chain. Of particular interest is not the average but the maximum pressure inside a chain. The maximum pressure is criterion whether a grain will be comminuted in a mill. As soon as one of the particles of a chain has broken the load discharges and the force chain vanishes and reappears elsewhere. Therefore here we are mainly interested in the maximum pressure. For low velocities $`\mathrm{\Omega }_I`$ and $`\mathrm{\Omega }_{II}`$ (fig. 2, top and middle) the average value of the maximal force in a force chain of length $`L`$ grows monotonously with the length, i.e. a force chain acts similar as a beam in a framework: all the weight and momenta from particles above the chain are supported by the chain and propagated through the chain. Therefore one finds in most cases that the lower the position of a grain the higher the load.

In contrast, the corresponding plot (fig. 2, bottom) for high velocity $`\mathrm{\Omega }_{III}`$ shows no characteristic dependence of the maximal force on the length of the force chain, i.e. for this case the force chains do not act as described in the previous paragraph. Thus we conclude that the mechanism of stress propagation in force chains is much less efficient in the case of high rotation velocity.

In fig. 3 the average number of particles per force chain which feel a pressure $`P_i`$ larger than a given threshold $`P`$ is plotted versus the length of the force chain. This number is directly related to the breaking probability which can be described using Weibull statistics $`^\mathrm{?}`$. Again one can see that longer force chains lead to a higher probability of comminution. For the highest rotation velocity $`\mathrm{\Omega }_{III}`$ the distributions are significantly lower again.

The time averaged spatial distributions of the pressure $`P_i`$ for $`\mathrm{\Omega }_I=2Hz`$ are shown in fig. 4 (left part of figure). The color codes for the local pressure $`P_i`$, where blue means low pressure and red color stands for high pressure. The upper plots include only particles which do not belong to a force chain, whereas the lower figures include only particles which do belong to a force chain.

The right hand side of figure 4 shows the time averaged spatial distribution of the number of particles, which feel a pressure $`P_i>3000gcmsec^2`$. These results compare directly to the experiment by Rolf and Rothkegel $`^\mathrm{?}`$ (see above). The maxima of the spatial distribution are located at the bottom of the cylinder near the wall which coincides with the results of the experiment.

Figures 5 and 6 are equivalent to fig. 4 where the angular velocity of the mill is $`\mathrm{\Omega }_{II}=10Hz`$ and $`\mathrm{\Omega }_{III}=19Hz`$, respectively. The color scaling is the same for all left side figures and the color code for the right figures is given by the color bar.

From these figures we conclude the surprising result that direct impacts of particles with high relative velocities at or close to the surface do not lead to high stresses. Therefore they are almost neglectable for the comminution process. Regions of high pressure can be found near the walls deep inside the material. Particles which do not belong to a force chain feel everywhere in the material a low pressure, only particle which are a member of a force chain feel larger stresses, i.e. only particles which belong to a force chain have a high probability to break. Fig. 6 shows that the mechanism of force chains is of low efficiency for rather high rotation velocities. The pressure is much lower than in the other figures 4 and 5, which leads to a diminution in the breaking rate as observed experimentally $`^\mathrm{?}`$.

For $`\mathrm{\Omega }=\mathrm{\Omega }_{III}`$ (fig. 6) the rates resulting from impacts of particles at the surface overcomes the rates resulting from the force chains. Nevertheless the rates are clearly lower than for slower spinning mill.

## 4 Conclusion

Using two-dimensional molecular dynamics simulation of spherical granular particles we investigated the motion of granular material moving in a ball mill which spins with different angular velocity. The numerical results which agree with experimental data by Rolf and Rothkegel $`^\mathrm{?}`$ lead to an explanation for the experimentally observed spatial distribution of pressure in such a ball mill. A network of force chains which permanently changes its structure has been observed in the rotating cylinder. Detailed and separate analysis of the statistical properties of these chains and of the pressure acting on particles which either do belong to a chain, or which do not belong to a chain lead to the conclusion that force chains play an important role for comminution processes in granular materials. Force chains lead to concentration of the dynamic and static load of a large amounts of grains in the bulk of the material to few grains close to the bottom of the container. These particles finally feel very large pressure. The described mechanism of the action of force chains is able to interpretate experimental data while other explanations fail $`^\mathrm{?}`$.

## References