Assessing contact forces in granular materials from experimental measurements of kinematics

We propose a new approach to assess interparticle forces in granular materials, based on the combination of experimental measurements and numerical techniques. Experimental measurements mainly consist in the assessment of the full kinematics of grains in a quasi-statically loaded assembly; such measurements are performed by applying Digital Image Correlation (DIC) on high-resolution photographs taken during the tests. This experimental data is then combined with the application of numerical methods which account for grain equilibrium. Two different approaches are adopted for this purpose: a Static ElastoPlastic Computation (SEPC) and the Non-Smooth Contact Dynamics (NSCD). Both methods aim to find a set of contact forces that is mechanically admissible, based on kinematics measured from experiments. In this way, we aim to make an assessment of forces that is experimentally based, as it comes from experimental measurements, and realistic, in the sense of being mechanically admissible.


Introduction
The full understanding of the mechanical behaviour of granular materials needs a complete description of what occurs at the grain scale.In particular, it is necessary to focus on the complex mechanisms of force transmission, whose peculiar nature is well known, but not widely investigated yet.The transmission of forces between grains in contact plays a key role in the microscale mechanics of an assembly of particles; its study can prove essential also for the comprehension of phenomena occurring at the local scale, such as grain breakage.Nowadays, thanks to numerical methods, we can easily reproduce and investigate such microscale behaviour; on the other hand, it is still not easy to make an experimental assessment of contact forces.Photoelasticity [1] is commonly used for this purpose; anyway, the use of this technique is limited by some constraints.Another solution, proposed by Andrade and co-workers [2,3], is the Granular Element Method (GEM), which aims to analytically determine the contact forces acting on a volume on the basis of the stress field inside the volume and of (at least) one of the forces on its boundary.The stress field should be derived, through the definition of an appropriate constitutive law, from the strain field, which can be determined by means of Digital Image Correlation (DIC) measurements.However, this method sets constraints to the grain size and consequently to the number of grains, since each grain surface has to be large enough to be sufficiently discretised to allow strain measurements.Recently, a new e-mail: mathias.tolomeo@3sr-grenoble.frInstitute of Engineering Univ.Grenoble Alpes approach has been proposed by Hurley et al. [4], representing the first attempt to measure interparticle forces in a non-idealised 3D granular material.This technique combines 3D x-ray diffraction (to assess grain strain, and derive grain stress tensors from it) with a numerical optimisation algorithm which infers the interparticle forces on the basis of these measurements and of grain equilibrium.
The method we propose is a combination of experimental measurements and a numerical approach.The experimental measurements consist in the determination of the full microscale kinematics of a 2D granular assembly, i.e., grain motion (displacement and rotation).Two numerical techniques are used: a Static ElastoPlastic Computation (SEPC), which views the assembly as a system of nodes, connections, and supports, and solves the problem as a conventional structural mechanics one, by looking at successive states of equilibrium during a test; the Non-Smooth Contact Dynamics (NSCD) [5], which is used here as an indirect method to determine contact forces based on contact velocities.Despite the different approach, both these techniques aim to find a set of contact forces which is mechanically admissible and compatible with the experimental information (grain kinematics, macroscopic loading).This indirect use of NSCD was already proposed in [6]; however, some issues of non-uniqueness of the solutions were observed.In our application, it might prove useful to use some experimental information such as the status of a contact (i.e., whether it is sliding or not), in order to overcome such issues.

Experimental approach 2.1 Experimental device
The 1γ2ε apparatus, described in detail elsewhere [7,8], allows to quasi-statically strain a 2D assembly of rollers (Schneebeli rods), with independent control on the macroscopic deformation components (vertical and horizontal strain, shear strain).Different kinds of tests were performed (isotropic compression, biaxial compression − vertical and horizontal −, simple shear).However, in this study, the results being discussed are mainly from the biaxial vertical compression test.Such test was carried out by first applying an isotropic loading with a mean stress of 100 kPa, and then increasing the vertical stress while keeping the horizontal one constant.The 2D material used for these tests is an assembly of around 1850 6 cm long wooden rods (Fig. 1), with four different diameters (8, 12, 14 and 20 mm).

Kinematic field measurements
During the tests, 80 Mpixel digital pictures of the sample were shot with a Phase One IQ180 camera, with a frequency of 5s, corresponding to a strain increment in the vertical direction Δε y ≈ 0.01%.Each grain's visible surface was painted with a speckle pattern (Fig. 1) in order to track its displacement and rotation.The tracking procedure was carried out by means of a particular technique based on Digital Image Correlation, optimised for the case of discrete materials, for which the motion can be erratic; we refer to this technique as Particle Image Tracking (PIT) [9,10].
Figure 2: Sketch of a contact between two grains.The "fictitious" overlap is defined by Eq. ( 1).We assume the radii to be constant, while some deformation can occur in a small area around the contact point.With these assumptions, it is straightforward that x i j can be negative when the grains in contact are compressed.

Contact forces assessment 3.1 Direct method
A first, rough assessment of contact forces was performed with a method which we refer to as "direct" or "DEMlike".This method aims to find a direct relationship between the contact force and some microscale kinematic indicator derived from the measurements.For this purpose, we introduced a microscale kinematic parameter which is defined as a "fictitious" overlap, following the overlap definition in common DEM.This definition is based on two assumptions: the grain is treated as a rigid body, so its radius stays constant; a very small deformation occurs in a narrow area around the contact point, when two grains in contact undergo compression.With these two assumptions, the definition of "fictitious" overlap comes straightforward (see sketch in Fig. 2): Then, a simple linear contact law was assumed between the overlap and the normal contact force (in compression), also based on the results of a specific compression test which was performed on a single contact.
Having defined the normal contact forces in this way, it was possible to make a first assessment of their distributions in our assembly.Of course this method was not intended to give a precise evaluation of the exact value of the forces one-by-one, due to some unavoidable uncertainties in the determination of the overlap (coming mainly from the non perfect roundness of the surface of rods), but still some features of interest could be found: the statistics of normal forces in the isotropic compression test reproduce very well those obtained in similar DEM simulations in [11] (see Fig. 3); preferential contact orientations, with respect to those contacts subjected to the highest extent of compression, are aligned with the direction of maximum compression for each test respectively (biaxial vertical compression, horizontal compression, simple shear), thus showing a good correspondence between the microscopic and macroscopic behaviour of the assembly.of the probability distributions of normalised "fictitious" overlaps ξ = x/ x .The linear fit of weak and strong overlaps, in the semilogarithmic and log-log plot respectively, show that the trend can be modelled through a power-law decay for weak overlaps and an exponential one for the strong overlaps.
Except for the previously mentioned sources of inaccuracy in the measurements, there are mainly two elements that this method is actually lacking in order to give a correct estimation of the whole set of contact forces.One is the history of tangential displacements, essential for a correct assessment of tangential forces; since we only have measurements in discrete time steps, we are missing a lot of information in terms of the tangential displacements occurring between two consecutive steps, and consequently in terms of the corresponding forces.The other missing feature is the mechanical balance of the grains: as it is not taken into account here, the set of forces obtained is not, in general, statically admissible.Introducing this concept, by means of the two numerical tools presented in the following sections, can lead to a more realistic solution.
The output of this direct method, although trivial and rough, is yet not to be discarded: using it as a starting point (guess solution) for the numerical tools that we propose can help to reach a solution, and in particular to orient the solution towards the "real" one, in case of non-uniqueness issues (as in [6]).

Static ElastoPlastic Computation
In the modelling of granular materials, quasistatic methods represent an important alternative to the most common numerical methods, which account for inertial effects (such as Molecular Dynamics and Contact Dynamics).Their peculiarity is to view the system, under a variable loading, as a succession of equilibrium states [12].In this way, the assembly can be considered equivalent to a system of nodes, connections, and supports, and the problem can be studied as a conventional structural mechanics one [13].Such methods are based on building a rigidity matrix G which describes the granular structure of the assembly, i.e., it links grain-related quantities (forces and moments acting in correspondence of each degree of freedom and grouped in the vector F ext ) to contact-related ones (contact forces grouped in the vector f), from both a kinematic point of view and a static one, such as in the generalised equation of force equilibrium where H = G t .The method presented here, which will be referred to as Static ElastoPlastic Computation (SEPC), is meant to find a set of contact forces which is, at the same time, statically admissible (i.e., in equilibrium with the external loading) and plastically admissible (i.e., fulfilling the adopted friction law).This is done by means of an iterative procedure of double projection, of an initial set of forces (guess solution), respectively on the subspace of statically admissible solutions and on that of plastically admissible one.In particular, since we assume a Coulomb friction law, the boundary of the subspace of plastically admissible solutions corresponds to the Coulomb cone.This computation algorithm is fully detailed in [12].

Validation of the method
A first validation of the method was performed on synthetic data produced by Molecular Dynamics simulations of biaxial vertical compression tests, on an assembly which resembles a 1γ2ε sample in terms of number of particles and size distribution.The purpose of this validation was to test the ability of the method to determine a set of contact forces, given an initial geometrical configuration and an external loading, which fulfils the two conditions of static admissibility and plastic admissibility.
Having applied the method on some thousands of intermediate steps of the simulations, the first outcome is that in most cases it is possible to find a solution (Fig. 4), which, for strong normal forces, matches very well the MD forces (with a maximum relative error of around 10% 1 ).However, some issues have been encountered, leading to an impossibility of finding a statically and plastically admissible solution in some of the intermediate steps analysed.In some situations, the problem was simply overcome by initialising the set of forces to a guess solution which is close to the expected one, e.g., using the set of forces from the previous step, instead of starting with all zero forces; this shows that the success of the method is affected by the choice of the initial guess solution, as introduced in Sec.3.1.
Yet this was not enough for some other situations, for which other issues were encountered.One typical issue, especially when dealing with MD simulations, is the occurrence of highly dynamic effects: in such situations, no static equilibrium could be reached, but the set of forces obtained was surprisingly consistent with the MD one (still with a maximum relative error for strong normal forces of about 10%), as in Fig. 4. Another typical issue of quasistatic methods, that should be taken into account when applying SEPC, is the stability of an equilibrium state, which can be studied through the condition that the second-order work is positive [12,13]. 1 The maximum relative error is computed as max The map on the left corresponds to a configuration for which the solution obtained is statically admissible; the map on the right, on the other hand, corresponds to a solution that is not statically admissible.In both cases, force chains clearly appear.Despite being not equilibrated, the set of forces in the configuration on the right is still clearly consistent with the external loading applied.

Non-Smooth Contact Dynamics
The application of Non-Smooth Contact Dynamics as an indirect method to assess contact forces from experimental kinematic data was introduced, at a preliminary stage, by [6].The idea is to apply NSCD to find a set of contact forces which is mechanically admissible, i.e., compatible with both the grain kinematics and the external loading applied.The kinematic information required by this method mainly consists in the relative velocities at the contact between each couple of grains; such information can be easily extracted from the grain kinematics assessment performed by means of the PIT technique.This method might result complementary to the one in Sec.3.2, as it can prove particularly useful for those situations in which the system undergoes dynamic events, and it is not possible to assume a static equilibrium.
The main issue concerning this method, as it was observed in [6], is the non-uniqueness of the solution, which sometimes can lead to a set of forces which might be mechanically compatible with the overall loading, but totally inconsistent with the kinematics, due to possible inaccurate inputs.To overcome such issues, some extra information can be exploited, e.g., the status of a contact (i.e., if it is a sliding contact or not).Another way to orient the method towards a "good" solution is to initialise the set of forces to a guess solution which should be close to the expected one, in the same way as in Sec.3.2.1.To this purpose, a possible initial set of forces can be represented by the one derived from the "direct" method in Sec.3.1.

Summary and outlook
In this study, we have presented the preliminary results of a combined experimental-numerical approach for contact forces assessment in granular materials.The experimental data comes from kinematic measurements performed by means of the PIT technique on a 2D rods assembly which behaves equivalently to a granular material.Such information is then combined with two different numerical methods, a quasi-static approach (SEPC) and the Non-Smooth Contact Dynamics, which, by injecting grain equilibrium, aim to find a set of forces which is compatible with the external loading and the grain kinematics.Some testing of the first method has already been performed on synthetic data from MD simulations, as a validation, proving the reliability of this method, and its robustness with respect to some inaccuracies in the inputs which may represent an issue when dealing with experimental measurements.In a next future, the applicability of both methods to experimental data will be addressed.

Figure 1 :
Figure 1: Image of an isotropically compressed sample in the 1γ2ε device, with zoom on a small amount of grains to show the speckle pattern on their visible face.The frame is 599.4 mm × 444.7 mm.

Figure 3 :
Figure 3: Semilogarithmic plot (a) and log-log plot (b)of the probability distributions of normalised "fictitious" overlaps ξ = x/ x .The linear fit of weak and strong overlaps, in the semilogarithmic and log-log plot respectively, show that the trend can be modelled through a power-law decay for weak overlaps and an exponential one for the strong overlaps.

Figure 4 :
Figure4: Maps of normal contact forces obtained with SEPC for an intermediate step (vertical strain ε y ≈ 5% for both − left and right − configurations) of a biaxial vertical compression MD simulation, where the thickness of each line is proportional to the normal force magnitude.The map on the left corresponds to a configuration for which the solution obtained is statically admissible; the map on the right, on the other hand, corresponds to a solution that is not statically admissible.In both cases, force chains clearly appear.Despite being not equilibrated, the set of forces in the configuration on the right is still clearly consistent with the external loading applied.