Numerical modelling of rough particle contacts subject to normal and tangential loading

Our understanding of the mechanics of contact behaviour for interacting particles has been developed mostly assuming that surfaces are smooth. However, real particles of interest in engineering science are generally rough. While recent studies have considered the influence of roughness on the normal force–displacement relationship, surface roughness was quantified using only a single scalar measure, disregarding the topology of the surface. There are some conflicting arguments concerning the effect of roughness on the tangential or shear force–displacement relationship. In this study, optical interferometry data are used to generate the surface topology for input into a 3D finite element model. This model is used to investigate the sensitivity of the normal force–displacement response to the surface topology by considering different surfaces with similar overall roughness values. The effect of surface roughness on the tangential force–displacement relationship and the influence of loading history are also explored. The results indicate that quantifying roughness using a single value, such as the root mean square height of roughness, Sq, is insufficient to predict the effect of roughness upon stiffness. It is also shown that in the absence of interlocking, rough particle surfaces exhibit a lower frictional resistance in comparison with equivalent smooth surfaces.


Introduction
Analytical expressions for the force-displacement relationships between contacting grains are required for accurate discrete element method (DEM) simulations. The Hertz-Mindlin model, which is used in many DEM studies, assumes that the contacting particles are smooth. While this model can capture key elements of contact mechanics, it is now clear that surface roughness has a measurable influence on the force-displacement relationship, particularly for normal loading and small deformation levels [2][3][4][5]. Several studies have emphasised the importance of roughness on the macroscopic soil stiffness (e.g. [6][7][8]). Otsubo and O'Sullivan [9] quantitatively related macroscopic stiffness to particle surface roughness in a combined experimental-DEM study.
The effect of surface roughness on small strain stiffness was investigated by Yimsiri and Soga [10] and Otsubo et al. [11] by means of analytical and DEM studies. Yimsiri and Soga [10] assumed that the tangential contact response is not influenced by surface roughness; this assumption was reformed in Otsubo et al. [11] by considering a reduction of stiffness in both normal and tangential force-displacement relationships as the surface roughness of contacting particles increases. In an experimental study using interparticle loading tests, Senetakis et al. [5] noted that tangential stiffness might not be significantly affected by surface roughness and Cavarretta et al. [3] observed a higher friction for rough contacts. Particle-scale experimental data on the effect of surface asperities is however very scarce given the difficulties in monitoring the position, deformation and force acting on the individual asperities during loading. In tribological research, which considers surfaces from a more general perspective, numerical prediction of normal contact stiffness in the presence of roughness and adhesion has been carried out using different numerical techniques, such as, the Boundary Element Method [12], DEM [13] and a multilevel multi-integration technique [14]. Regarding tangential stiffness, Medina et al. [15] showed that tangential stiffness is proportional to the normal load and independent of the asperity radius and the Young's modulus of the material. This is in agreement with experimental observations from Berthoud and Baumberger [16], which used a multilevel multi-integration technique to approximate the frictional energy dissipation.
The use of the finite element (FE) method to quantify the effect of roughness has been previously reported in tribology literature (e.g. [17][18][19][20]). The present paper uses a micro finite-element (μFE) model proposed by Nadimi and Fonseca [21,22]. This numerical approach uses a more accurate representation of the grain morphology, including contact topology and grain roughness; which is critical to advance our understanding of contact behaviour and to derive empirical parameters for more realistic contact models. In this paper, the numerical model and the effect of surface roughness on the normal and tangential force-displacement relationships are described prior to considering the effect of loading history and contact area on friction. This is followed by a comparison between analytical and numerical estimates of contact area.

Modelling surface roughness
Rough surfaces can be generated using a random algorithm (e.g. [23][24][25]). The present study focuses on rough surfaces that were artificially prepared in the laboratory using a milling process. These surfaces better represent realistic surfaces that exist in nature than a randomly generated rough surface.
The modelling of surface roughness comprised three main steps: (1) measuring roughness and developing the 3D topographical map, (2) converting the 3D topographical map into a numerical mesh, and (3) solving the partial differential equations using FE in order to calculate the stress and strain in the model.

Roughness measurement
The particles considered here were glass ballotini with a diameter of approximately 1.2 mm. The choice of spherical particles is justified by the need to exclude the effect of particle shape on the investigation of particle roughness. The manufacturer-supplied ballotini is considered to be only nominally smooth. Artificially rough ballotini particles were produced by milling the manufacturer-supplied ballotini following the technique described in Cavarretta et al. [26]. The surface roughness for each particle was controlled by varying the milling time [27]. Table 1 summarises the 8 particle surfaces considered here: one manufacturer-supplied bead (named MS01); three medium-rough beads (named MR01, MR02 and MR03), which were milled for 5 h; and, four very rough (VR) beads (named VR01, VR02, VR03, VR04), which were milled for 25 h. More details of the process can be found in Otsubo [27]. A computer generated, perfectly smooth surface named PS01 was included for comparison of the numerical model with existing theories.
The roughness measurements were made with a Fogale Nanotech optical interferometer [28]. As illustrated in Fig. 1, data are acquired by transmitting white light from a source located above the sample. The light is split into two half beams. The first half-beam reflects off the surface being measured; this is compared with the second half-beam that is reflected from a reference mirror. A charge-coupled device (CCD) camera captures the wave interference corresponding to the difference of length between the paths of the two beams. These data are used to generate a 3D surface map. The vertical resolution of the apparatus is in the order of 10 nm, depending on the reflectivity of the surface. Here, the motif analysis programme available in Fogale 3D Viewer [28] was used to remove the surface curvature with a shape motif size of 25% of the window length of the measurement. The X, Y and Z coordinates of the points on the surface were extracted for the analysis. The X and Y values are determined on a regular square Cartesian grid, with a grid interval of 0.184 μm in both directions. Further details on the application of this technology are given in Cavarretta [29] and Otsubo [27]. This measurement approach was also applied by Altuhafi and Coop [30], Yang et al. [31] and Yao et al. [32] to measure the roughness of sand grains. The roughness maps measured for each particle were obtained for 106 × 106 μm 2 region of interest, as shown in Fig. 2. Roughness was quantified here using the root mean square height of the surface, S q , and root mean square gradient of the surface, S dq , given by: where m = number of discrete data points; and Z i = elevation relative to the reference surface. Table 1 includes the S q and S dq values for each bead. Yang et al. [31] proposed an alternative, fractal-based approach to quantify roughness using interferometry data. In the present paper, the more widely accepted metric of S q is used to better contextualize the research with a broad range of previous studies [19].
Images of two representative surfaces are presented in Fig. 3. Figure 3a, c correspond to the rough particle VR04 with S q = 767 nm, while Figs. 3b, d, represent the nominally-smooth (manufacturer-supplied) particle MS01 with S q = 96 nm. Figure 3a, b are contour plots where the grey level indicates the relative elevation of each point rendered. Figure 3c, d are cross-sectional profiles through the data presented on Fig. 3a and 3b, respectively. Each line profile is taken at the centre of the contour image, i.e. at a y value of 53 μm. Fig. 1 Schematic diagram of the process for acquiring the interferometry images, as developed by Fogale [28] and Yao et al. [32] CCD Camera

FE mesh generation process
As explained above, the data were extracted on a regular grid in the X-Y plane with a grid interval of 0.184 μm. To generate the μFE mesh the data were converted into a 3D volumetric matrix, where each cubic cell in the matrix has dimensions of 0.184 × 0.184 × 0.184 μm 3 . The volumetric matrix was built as stacks of cubes starting from a reference plane located beneath the lowest point of the surface.
Each stack was centred on a given X, Y grid point and constructed sequentially by adding cubes indexed as "1" until the midpoint of the next cube added was higher than the Z coordinate recorded for that grid point. That cube and all subsequent cubes were indexed as "0". For each stack, the elevation of the top non-zero cell is denoted Z max . The upper limit of the stack was set as 9.2 µm (50 cells). This algorithm generates a binary volume with solid and void elements. This volumetric matrix was used to generate the FEM  Delaunay refined meshing mesh as detailed in Fig. 4. A refined Delaunay triangulation algorithm was used and more details can be found in Nadimi and Fonseca [21]. To improve computational performance, a small volume was used instead of the full sphere. This simplification however restricts the amount of normal force that can be applied since the deformation zone must be within boundaries of the small volume considered. A cross-section through this volume of interest (VoI) is shown in Fig. 5. To ensure that the nodes associated with the asperities were sufficiently distant from the domain boundaries, a 44 μm thick layer of cells was placed below the assembly of cubes. The appropriateness of this model configuration is confirmed by the parametric study detailed in the following Section. The nodes and elements were generated in MATLAB (Mathworks, 2016) and imported into Abaqus software package [33] using an input file containing the nodal coordinates and all elements forming the mesh. Each mesh contains approximately 500 thousand nodes and 3 million tetrahedral elements.

Numerical model description
The boundary value problem considered is illustrated in Fig. 5. The numerical formulation was defined in the framework of a combined finite-discrete element model using a dynamic explicit formulation [21] that employs a central difference time integration scheme as implemented in Abaqus. Previous studies that used the finite element method for contact mechanics problems include Vu-Quoc and Zhang [34], Li et al. [35] and Rathbone et al. [36]. The contact algorithm implemented in Abaqus generates contact forces to resist node-to-element penetration and uses a finite-sliding formulation that allows for arbitrary separation, sliding, and rotation of the surfaces in contact [33]. The finite-sliding formulation assumes that the incremental relative tangential motion between surfaces does not significantly exceed the dimensions of the master surface (in this case the surface of the rough particle). However no limit is imposed on the overall relative motion between surfaces. In this study, a coefficient of friction of 0.2 was set for shearing. The physical and mechanical parameters of the glass ballotini were assigned to the model, as listed in Table 2. This includes the Elastic modulus (E), Poisson's ratio ( ) and density. Normal loading was first applied to the domain by moving the rigid platen illustrated on Fig. 5 downwards until a target normal force was achieved. Shearing was simulated by inducing horizontal movement on the rigid platen. In order to investigate the effect of loading history on the tangential stiffness, an elastic perfectly plastic constitutive model was assigned to the particles. This has a minor effect on the normal stiffness of glass beads as shown previously by Nadimi and Fonseca [22].
To confirm that the thickness of the analysis domain, i.e. 44 + Z max μm, was sufficient, the stress distribution throughout the model was investigated. Figure 6a, b show two YZ cross-sectional views of the stress distribution in the domain for the model of VR04, at x = 30 µm and x = 70 µm respectively. The location of these cross sections relative to the specific features on the surface topology can be found using the x-axis coordinates in Fig. 3. Stress concentration can be seen around the asperities with a maximum value of 130 MPa. At the bottom of the domain, the stress value is zero, which means that the stress distribution was not affected by the boundary conditions. This confirms that the choice of 44 + Z max μm for the domain thickness was reasonable. For the semi-smooth and smooth samples considered, the stress concentrations are smaller in extent; a maximum stress concentration value of 68 MPa was observed in MS01. Figure 7 illustrates load-deformation response during the application of normal loading. It can be seen that the load-deformation response for PS01 perfectly matches the Hertzian theory [37,38], which can be taken as a validation of the model used here. It is clear that there is a systematic variation in the response with increasing surface roughness. The manufacturer-supplied particle MS01 shows a response that is broadly similar to the Hertzian theory prediction (albeit slightly softer). For the mechanically roughened  particles (cases MR and VR), an initial very soft response was induced as evidenced by the low initial stiffness. This behaviour is in agreement with the observations reported by Cavarretta et al. [3]. Comparing the MR and VR datasets, it is clear that an increase in surface roughness results in a softer contact behaviour. It seems also that the variability in the response increases with increasing roughness; there is closer agreement amongst the MR data than amongst the VR data. The scatter in the VR data does not link directly to the S q values; VR01 has the lowest S q value and should have, therefore, the stiffest response. This is, however, not the case. The variability in the data at large S q values indicates that the S q measure alone cannot be used to accurately predict the load-displacement response. These observations support the use of a roughness measure that can account for variability, such as that proposed by Yang et al. [31]. Greenwood and Tripp [39] showed that rough surface contact responses can be described by the Hertzian theory for contact response if the normal load exceeds a threshold normal force, N GT :

Load-deformation response
where E* is the effective contact stiffness given by and R 1 and R 1 are the radii of the contacting spheres. The particle compression tests documented by Cavarretta et al. [3] conformed to this theory. Greenwood and Tripp [39] proposed that the N GT value varies proportionally with S q 1.5 , and in this study the semi-smooth case approaches the N GT value (≃ 7.5 N), while rougher cases do not reach their N GT values (> 35 N). It is important to clarify that the load in the simulations considered here was lower than the N GT values.
The effect of surface roughness typically reduces with increasing contact normal force due to deformation and yield of the asperities; and so, the contact load-deformation behaviour will vary in a simple compression test. Otsubo et al. [40] proposed a mathematical expression to model rough contact responses considering three regions of behaviour: (1) asperity-dominated (Eq. 4), (2) transitional (Eq. 5), and (3) equivalent Hertzian contact (Eq. 6) where the region shifts with the overlap (δ) of the two surfaces relative to Sq, given by: Otsubo's model was previously verified by considering the experimental data by Greenwood et al. [41] and analytical 24.47 S q ≤ Fig. 7 Comparison between theory and the results from the numerical simulations in terms of normal force-displacement for the smooth and rough glass beads described in Table 1 Fig . 8 Comparison of the results from the numerical simulation and the data from Otsubo et al. [40] model, in terms of normal force-displacement, for a the medium rough glass beads (MR set) and b the very rough glass bead (VR set) expressions by Yimsiri and Soga [10]. The normal force-displacement relationships for the analytical model are compared with the simulation results in Fig. 8a, b for MR series and VR types, respectively. In the model, the diameter of the glass beads was assumed to be 1.2 mm. The range of normal forces considered in Fig. 8 is approximately N GT /20 and N GT /50 for the MR and VR types, respectively. Thus, the contact responses are in the transitional region of response that exists between asperity-dominated contact behaviour and Hertzian contact behaviour as proposed Otsubo et al. [40]. When comparing the simulation data with Otsubo's model, as shown in Fig. 8, the initial stiffness values at small displacements are in good agreement for both MR and VR types. This suggests that the simulation results capture the effect of asperity shape on the force-displacement response when compared with the averaged response predicted by the model. However, Otsubo's model gives a stiffer response as the displacement increases, when compared with the simulation results, especially for the VR particles. This is most apparent for sample VR01 (simulation) which has a significantly softer response similar to VR04, despite having a lower S q value of 514 nm compared to S q = 767 nm for VR04. Referring to Fig. 9, the surface of the sample VR01 exhibits distinct spikes and it takes a larger displacement to reach a point where the response conforms with Hertzian contact mechanics. In contrast, Otsubo's model assumes a Gaussian distribution of the surface asperities following Greenwood and Tripp [39]. In other words, while the initial contact in the FE model is taken at the interaction between the tips of two asperities that may not be the case in the experiments used to develop Otsubo's model, where some interlocking may occur. What is evident from the simulations is that it is not only the average surface roughness that is important; each surface has a force-displacement response that depends on the specific topology of that surface. For the surfaces considered here, the effect of this variability is finite relative to the overall force-displacement relationship.

Contact area
Hertzian theory predicts that the contact between two smooth spheres with effective radius R* is a circle with radius a 0 , where a 0 is given by: As indicated above, the Greenwood and Williamson [42] and Greenwood and Tripp [39] models approximate the rough surface by considering spherical asperities of identical radius with a Gaussian distribution of heights. The idealisation of Greenwood and Williamson theory, which leads to nearly linear variation of the real contact area with F N , was debated by Persson et al. [43] and Campaná and Müser [44]. Subsequently, Pastewka and Robbins [1] proposed an analytical relationship to relate the contact area ( A PR ) and the normal force: where κ is a constant that typically takes the value of 2.0 [45]. Recently, Müser [45] compared the approximation of Pastewka and Robbins solution with numerical data and observed that the theory predicts the real contact area with less than 10% error. Müser has also improved the original formula by cancelling the mean-field approximation, for a better scaling at large loads, for which Hertzian theory dominates.
In the numerical simulations performed in this study, the actual contact area can be quantified by specifying 'CAREA' in the Abaqus history outputs. Figure 10a shows the numerical results of MS01 in terms of both contact force and contact area versus normal displacement. Figure 10b-e provide schematics of the contact area at different normal displacement values, δ n = 0.6, 1.2, 1.8, and 2.4 µm, respectively. In the schematics, the contact area is formed by the elements with contact pressure higher than zero, as illustrated using the grey colour. Figure 11a shows the numerical results of the VR04 simulation in terms of both contact force and contact area versus normal displacement. Figure 11b-e present representations of the contact area at normal displacement values of, δ n = 0.6, 1.2, 1.8, and 2.4 µm, respectively. The differences in contact area caused by the physical asperities can be clearly seen by comparing Figs. 9 and 10.
The relationship between contact area and normal force is plotted for different values of S dq , as shown in Fig. 12. The Hertzian response and numerical predictions of the rough particle response are illustrated by a grey dash line and a black bold dash line, respectively. The good agreement observed between the model and the theoretical equation, suggests that future development of contact laws for DEM models should be based on the reduced contact area as proposed by Pastewka and Robbins [1]. Figure 13 compares the results of the simulations for VR04 under tangential loading for different normal load values: F N = 2, 4, 6, and 8 N. Mindlin [46] and Mindlin and Deresiewiez [47] investigated the elastic deformation of two contacting spheres under tangential loading. Based on their results, the tangential force-displacement can be described as follows:

Tangential load-displacement response
where µ is the friction coefficient, t max is the maximum t angential def lection before sliding, t max = 0.5 n (2 − )∕(1 − ) , and the condition for sliding to occur defined as | | t | | ≥ t max . The response predicted by Mindlin and Deresiewicz (M&D) theory adopting the same mechanical properties is shown in Fig. 13 for comparison (denoted as theory).
It can be seen that in all cases, the tangential load required for sliding of the rough surface is about 80% of the load predicted by the theory. This means that, although the   Figure 14 compares the results of simulations under tangential loading with M&D theory for all surfaces considered. The ratio of tangential force to normal force is considered. It can be seen that PS01 perfectly matches the theory. It is clear that an increase in surface roughness results in a softer contact behaviour. However, the MR surfaces and two of the VR surfaces have similar responses; this contrasts with the lower stiffnesses observed for the VR surfaces when compared with the MR surfaces under normal loading (Fig. 6). As before the response depends on the surface topologies of the individual particle. For these data, sample VR02 is an outlier, having a response that is significantly softer than the other VR and MR surfaces. These data provide further evidence that quantifying a surface only in terms of Sq value cannot effectively predict the load-deformation behaviour.

The effect of loading history
Contacts have a short life, constantly, new contacts are being created and others destroyed (e.g. [48]). Load reversals are also possible. Therefore, it is important to investigate the effect of loading history on contact behaviour. For this purpose, an elastic numerical simulation is not adequate because it is necessary to account for the deformation of asperities under loading. It was assumed here that plastic behaviour initiates at 10 MPa stress using an isotropic hardening model and the material was allowed to harden up to 110 MPa at 0.05 strain (hardening modulus, E t = 2 GPa), after which it behaves perfectly plastic (obtained based on curve-fitting and it is consistent with the plastic assumption in Nadimi and Fonseca [21]). Three simulations were carried out to investigate the effect of loading history on tangential force-displacement. The rough surface was loaded up to F N = 2, 4 and 8 N. Subsequently, the normal load was reduced to half of its initial value, that is F N * = 1, 2 and 4 N respectively (superscript * denotes loading history due to reduction in normal loading). Finally, the model surface was sheared under reduced normal loading. Figure 15 shows the result of the simulations under tangential loading. To facilitate the interpretation of results, M&D theory, roughelastic response, and rough-elastoplastic response with unloading history, under a given normal load, is presented (Fig. 16). It can be seen that the contact with loading history requires higher tangential force to slide when compared with a 'virgin' contact. This can be attributed to the larger contact area in the rough-elastoplastic model in comparison with the rough-elastic model. These observations indicate that understanding of contact behaviour can be improved by means of a theoretical method that relates contact area and surface roughness.

Conclusions
The contact behaviour of rough particles was investigated here using a combined discrete finite-element approach coupled with optical interferometry. The following conclusions can be made: 1. In the normal direction, contact stiffness clearly reduces with increasing surface roughness. Furthermore, the observed results indicate that the specific topology of the surface influences the contact response, and so, using a single value of S q is not sufficient. This suggests the need Fig. 14 Influence of surface roughness on tangential load-deformation behaviour presented by comparing the response of perfectly smooth case (theory and PS01), medium rough glass beads (MR set) and very rough glass beads (VR set)

Fig. 15
Comparison of the M&D theory and tangential force-displacement behaviour for a very rough glass bead (VR04) that experienced a significant reduction in the normal load (F N ). The superscript * denotes loading history due to reduction in normal loading Fig. 16 Comparison between the M&D theory and a very rough glass bead (VR04), in terms of tangential force-displacement, including elastic and elastoplastic behaviour with unloading history under normal force F N * = 2 N to use more sophisticated measures of surface topology, such as, the fractal-based approach proposed by Yang et al. [31]. 2. The contact model developed by Otsubo et al. [40] predicts a stiffer response when compared with the results from the FE simulations presented here, and the observed discrepancy increases with increasing surface roughness. Direct comparison of the two approaches is not simple since the FE analyses assume the initial contact to be defined by the extremities of the asperities, while this may not be the case in the experiments on which Otsubo's model is based. The expression proposed by Pastewka and Robbins [1] can capture the variation in contact area with normal force and may provide an improved approach to evaluating surface roughness effects on the normal force-displacement relationship. 3. In the tangential direction, it was observed that rough surfaces have a less stiff response than smooth surfaces. At a given normal force, the tangential force at which sliding occurs is lower for rough surfaces than for perfectly smooth surfaces. However, a clear trend of stiffness reducing with increasing roughness is not obvious; most of the MR surface responses are very similar to the VR sample responses. 4. The effect of loading history was also investigated by introducing plasticity into the finite-element model. The simulations show a slight increase in the tangential load at which sliding occurs due to flattening of asperities.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.