Stress-Strain diagrams for non-convex particles

While most granular materials in nature and technology consist of non-convex particles, the majority of discrete element (DEM) codes are still only able to cope with convex particles, due to the complexity of the computational geometry and the occurrence of multiple contacts. We have reengineered a code for convex polygonal particles to model non-convex particles as rigidly connected clusters. Constricting non-convex particles along the symmetry axes by 30% leads to an increase of the materials strength of up to 50%.


Introduction
Most discrete element methods use convex particles: Tangential forces are taken up only by (dry, Coulomb-) friction forces.Technically, in many applications like railway ballast, non-convex grains are preferred, as the higher strength is well known.For computer simulations, the controlled modeling of non-convex shapes is tedious.For clusters of round particles [1,2] or ellipsoids [3], the continuous variation of the non-convexity and the influence of interlocking is difficult to quantify.Rigid body simulations [4] correspond to the limit of vanishing strain, which for strength tests (like the bi-axial compression we want to investigate here) are as unsuitable as softly connected particles [5].We have reengineered a simulation of "soft" convex polygons in rigid non-convex clusters, overlap algorithms and force laws for non-convex are unfeasible.

Force computation
(Justifications and more details about for the interaction law beyond the short summary here can be found in Ref. [6].)For convex polygons, the elastic force is proportional to the overlap area A and the Young's modulus (in units of [N/m] for two dimensions).The "characteristic length" is e-mail: hg@mce.uec.ac.jp for the contact vectors r 1 , r 2 (see Fig. 1) between the centers of mass and the contact point.This definition fixes the propagation speed of (linear) deformations perpendicular to the sides of rectangles in space-filling packings to the sound velocity c = Y/ρ.Then the elastic force is the normal dissipation in analogy with the linear oscillator with the reduced mass m of the contacting particles with a damping constant γ.In contrast to the linear oscillator, in the DEM the dissipative force may not overcompensate and reverse the elastic force, which introduces unphysical attraction and noise.Therefore, a cutoff is needed so that the normal force cannot become attractive.Finally, a static (Coulomb-) friction force with friction constant μ is necessary to maintain the basic properties of physical grains (particles sticking on a slope tan α< μ, forming stable heaps also on flat surfaces).As friction law, we use the Cundall-Strack-model [7] for Coulomb friction where the tangential friction force between two particles is incremented proportional to their relative sliding distance, and a cutoff μN is introduced proportional to the normal force  above).For polygons, Coulomb friction is much more effective than for round particles which can escape sliding friction by rolling.The parameters are given in Tab. 1.For the theoretical contact time of a square particle with the smallest mass m min and Young's modulus Y, our simulation is stable with integration time-steps below 1/10T.Y in DEM-simulations should not be understood as that of the bulk Y bulk , but as the effective elasticity of the actual contact in the presence of micro-roughness as in Fig. 2 d): The actual contacts are much "softer" due the thinned out material density at the interstices, so the use of Y < Y bulk is appropriate.
For elastically connected monomers, one is faced with a choice between two evils: Either one has to use very soft connections to preserve the time-step for the monomers [5].or uses more rigid connections with much smaller integration time-steps.Instead, we implement rigid clusters where the the equations of motion (including rotation) act on the center of mass of the (rigid) cluster, while the interaction laws (forces, torques) between monomers of different clusters are the ones above for convex polygons, as layed out in Ref. [6].Additionally, only a data structure is needed which assigns monomers to clusters and vice versa.The definition of the characteristic length (see Fig. 3) is dealt with in Sec. 4. Different from clusters of round particles [1], polygonal monomers of the same particle (see Fig. 4) don't overlap, so there is no local variation of the Young's modulus due to multiple overlaps, with unforseeable consequences for the stability.

Simulation Geometry
For biaxial compression (see Fig. 5) we have implemented an approximately quadratic cell with 1400 particles.In experimental three-dimensional triaxial compression, a rubber membrane under water-pressure forms the side-walls and the stress is measured on the floor.Instead, we use straight walls with a force equilibration and measure the Composite particles with various degrees of nonconvexity, obtained by constricting the corresponding convex particle with the same elongation by 0%, 15 %, 20 % and 30 %, abbreviated in the following as c-0%, c-15% . . . .stress σ 1 on the ceiling.(This may lead so deviations from Ref. [8] where the measurement is on the bottom.The pressure on the floor is not equal to the lid pressure plus the weight of all particles combined because arching with the frictional contacts on the left and right wall deflects a part of the weight.)The bottom is fixed, while the lid moves down at constant velocity (0.01[m/s]) and the side walls are held at constant pressure σ 3 = 10000[N/m].(For artifacts resulting from more rectangular geometries or one side wall fixed, see [8]).The granular particles move according to Newton's equation of motion (2nd order Gear Predictor corrector as time integrator) while the positions x of the side-walls are determined from the external pressure p and the sum of the forces on the side-walls F as with M = 0.002 × m total to mimic the adaption of a rubber membrane to the whole system.(In Ref. [8], we had used M = 0.02 × m total , the differences are insignificant.)Trying position the side-walls at the equilibrium position with a zero-order approximation of Newton's equations of motion [9] fails as for a constantly deforming cell, there is no equilibrium position, but an (on average) constant drift.

Convex particles
For a base radius of r = 0.0055 [m], we set up regular dodecagons with elongations = 1.0, 1.4, 1.8.and a size distribution for the half axes with different, equally distributed random numbers rand ... .The resulting size distributions (see Fig. 6 for nonelongated particles, the ones for elongated ones were similar) were sufficiently wide to prevent ordering in latticelike structures, which could decrease the material strength considerably due to slipping along the "crystal axes" of the particle packing.Throughout the paper, the same sequence of random numbers is used, so initially, the particle areas are the same, only the elongations and the degree  of convexity are different.For increasing elongation, typical stress-strain curves and the behavior of the volume (inverse proportional to the packing density) are shown in Fig. 7.For all systems, the initial positions and relative particle sizes were generated with the same random number seeds.In this and the following graphics, the value for every thousand time-step is shown (the difference to a plot with every time-step was insignificant).The curves are aligned with respect to the linear regime.The peak of the stress-strain diagram for the dodecagons will usually not appear for circles.In the upper part of Fig. 7, one clearly sees the typical patterns of stress-strain curves for dense granular materials, the linear regime with constant gradient, the peak strength and the failure regime.In the lower part, we have shown the volume (area) occupied by the cell: The peak strength in the upper part corresponds to the minimal volume.One issue for the use of composite particles is the definition of the characteristic length: It is either possible to define it with respect to the center of mass of the contacting monomers or with respect to the center of mass of the whole particle, see Fig. 3.We will discuss both possibilities and the relation of the results in this paper.The  modellization is stable if penetrating contacts are treated as outlined in Ref. [6].When the approach it generalized for polyhedra in three dimensions, accordingly the most demanding task is the appropriate treatment of penetrating contacts.Our simulation advantages over that for clusters of round particles: Clusters for the particles and simple "monomer" polyhedra for smooth walls allows an implementation with the minimum amount of particles, without the remaining doubt whether the walls are smooth enough.When we go from simple to composite convex particles, with the same outline, there is the question whether the properties of the aggregate will change due to the internal "seams" in the particles as shown in Fig. 8.As one can see in Fig. 9, for composite particles, the corresponding stress increases slightly, and the stability (less fluctuations, delayed onset of the "breakdown") is slightly improved.(On a side-note: This means that a surface uneveness in spherical beads, e.g.scratches on glass beads or flashes on plastic beads, may increase the mechanical stability in comparison to simulations of spheres with ideal surfaces.)The definition of the inverse characteristic length with respect to the monomers of half of the size of the original particles (blue line) corresponds to an effective doubling of the Young's modulus (green line).

Non-convex particles
With the non-convex shapes from Fig. 4 and respective elongations of 1, 1.4 and 1.8, we obtain the stress-strain diagrams in Figs. 10, 11 and 12.The material strength (gradient of the linear regime, position of the maximal obtainable stress) increase with particle elongation and nonconvexity, as do the fluctuations in the data.As in the case of convex particles, the same seeds are used for the initialization, so that the relative volume (area) of the particles is the same throughout all simulations.For an increase of the non-convexity for the same elongation, the material strength increases up to 50The effect comes from the tangential forces which for non-convex particles are picked up and deflected also by the interlocking surface geometries via the elastic normal forces.This mechanism is more stable against load, slip and vibrations than the Coulomb friction between smooth surfaces alone.Nevertheless, the Young's modulus must be sufficiently large so 0 0.02 0.04 0.06 0.08 0.1 0.12 4.5 that interlocking asperities are not "pulled through" each other.For the external pressure uses, a Young's modulus of 1/4th would have annihilated effect.So for non-convex particles, a prudent choice of the Young's modulus is more important than for the convex case.

Summary and conclusions
We have successfully adapted a simulation for convex polygonal particles to non-convex clusters of such particles.For polydisperse dodekagons inscribed in ellipses, a variation of the elongation from 1 to 1.8 increases the material strength by a factor of ≈1.7.Constricting the particle shape along the symmetry axes by up to 30% increases the material strength by a factor of ≈1.5, so a total variation of the strength of up to a factor of ≈2.5 is possible.

Figure 1 .
Figure 1.Overlap computation with centers of mass C 1 , C 2 (black ×), force point F (center of the overlap polygon, gray ×) and contact vectors r 1 , r 2 .The dashed line indicates the tangential contact direction, given by the intersection points of the outlines.

Figure 2 .
Figure 2. Rigid contact in a), idealized deformed contact in b), DEM overlap in c) and actual contact with micro asperities in d).

Figure 3 .
Figure 3. Possible definitions of contact vectors and characteristic lengths, from the center of mass of the monomers (left) and the centers of mass from the clusters (right).

Figure 7 .
Figure 7. Stress-strain curves (measured at the lid) (above) and variation of the volume (below) for increasing particle elongation of a system with smooth size dispersion.

Figure 8 .
Figure 8. Non-composite convex particles experience tangential resistance to sliding (left) only via Coulomb friction, while composite particles (right) may interlock with a normal component.

Figure 9 .Figure 10 .
Figure 9. Stress strain diagram for non-composite dodecagons with the base Young's modulus Y (red) as well as 2Y (green) and for composite polygons of the same shape with two different definitions of the inverse characteristic length.

Figure 11 .
Figure 11.Stress strain diagram (above) and volume behavior (below) for non-convex particles with elongation 1.4.

Figure 12 .
Figure 12.Stress strain diagram (above) and volume behavior (below) for non-convex particles with elongation 1.8.

Table 1 .
Parameters for the computation of the interaction for 1m-long rod like particles.