Three-dimensional numerical assessment of axial and torsional interface shear behaviour

Interfaces between geo-materials and soils play a critical role in a wide spectrum of geotechnical structures and soil/site characterization techniques in geotechnical engineering. Consequently, understanding the mechanics of interface shear behaviour at different scales can benefit both soil characterization and the design of geotechnical systems. This paper presents a series of numerical simulations that utilize the 3D discrete element modelling (DEM) technique and compares the results with those obtained from laboratory counterpart tests under axial and torsional axisymmetric interface shear. The difference observed in macroand meso-scale responses under these loading conditions, such as shear strength, volumetric change, and shear zone characteristics are evaluated. In addition, responses at microscale including particle displacement trajectory, particles rotation, and local void ratio evolution are assessed allowing for links to the results obtained at larger spatial scales. These 3D numerical model studies expand the micromechanical processes under different shearing conditions previously presented by the authors from 2D to 3D space.


Current practice in interface testing
In the past three decades, studies on the mechanisms of interface shear in a variety of geotechnical problems have drawn significant attention. In the current practice of geotechnical engineering, interfaces between soil and construction geomaterials (e.g., concrete, steel, geotextiles, geomembranes, etc.) were designed to fulfil mechanical requirements in the applications of tunnels, piles, pavements, foundations, earth retaining structures, and so on. Similarly, many in-situ tests and laboratory experiments were also driven by or designed for quantifying shear strength at interfaces. For example, the modified direct shear device with large shear box has been a standard for examining the interaction between soil and geotextiles during shearing [1]. The cone penetration test (CPT) was designed to use the measured sleeve friction-the acting force at the equipment-soil interface-as a robust indicator for site characterization in-situ. However, the sleeve friction, fs, is not always reliable since the measurement is dramatically affected by the lack of control on the sleeve roughness as well as the soil disturbance induced by the cone tip insertion.

Recent research on interface shear
In response to the need in CPT practice to accurately evaluate interface shear strength, multi-sleeve friction attachments (MFA/MFPA) were developed at the Georgia Institute of Technology [2,3]. The devices are configured to equip sleeve attachments with controlled surface roughness at four distinct locations behind a conventional 15 cm 2 CPT module ( Figure 1a). As shown in Figure 1b, the surface roughness of these exchangeable sleeves was controlled by the dimensions and layout of the staggered asperities [4]. These four additional measurements provide the multi-sleeve devices the capability of simultaneously collecting sleeve friction measurements against surfaces of varying roughness from one CPT sounding. The most recent generation of the multi-sensor device-the multi-sleeve piezo-friction-torsion attachment (MPFTA)-provides the ability to shear the soil in torsional direction in addition to the existing axial shear mode [5]. Moreover, pore pressure measurements are enabled through the sensors before and behind each friction sleeve. These readings can be used to achieve a more robust evaluation of the mechanical and hydraulic properties of the tested soil layers. The aforementioned references provide detailed information regarding the design of the multisensor tools.
Numerical simulations are valuable alternative approaches to not only qualitatively but also quantitatively evaluate shear behaviour at interfaces due to their low cost and high reproducibility. Discrete element method (DEM) is well suitable to model boundary value problems involving large displacements and localizations, for the fact that soils are modelled with discrete entities that are subject to Newton's second law, which agrees well with the natural state of soil particles. Huang & Ma (1994), Jiang et al. (2006), and Butlanska, et al. (2014) pioneered research in DEM simulations of CPT in both two-dimensional (2D) and threedimensional (3D) space [6][7][8]. To extend these studies, 3D DEM models were developed to unveil the links between micro-and macro-scales responses at particulate-continuum interfaces.

Laboratory test setup
A series of axial and torsional axisymmetric shear tests between Ottawa 20-30 sand and a textured sleeve with a maximum surface roughness Rmax = 1mm were performed. The specimens were prepared to relative density Dr = 65 % and sheared against the textured sleeve at a rate of 5 mm/min under a constant 50 kPa confining pressure. The schematic of the testing equipment used in this study, the axisymmetric shear device, is shown in Figure 2a. The body of the device is composed of a cylindrical steel chamber and a two-piece smooth aluminium rod connected by a customized exchangeable friction sleeve. The steel chamber provides uniform lateral confining pressure to the soil sample via external air pressure applied on the latex membrane around the soil sample. Fixed aluminium plates close the chamber from top and bottom with rubber seals such that a BC4 boundary condition is established. At the centre axis of the chamber, the assembled rod-friction sleeve module is in position. The sand is air-pluviated around the module with predetermined density to simulate a "perfect insertion" scenario that enjoys the absence of the disturbed zone close to the penetrometer resulting from probe insertion. The axial and torsional loads and shear displacements are recorded for further analysis [9].

Numerical approach
To closely simulate the laboratory experiments, a 3D DEM model (PFC suite 5.0) was built, which is shown in Figure 2b. The lateral chamber and the planar caps at both end (not shown) were modelled as rigid walls with low friction coefficient, through which 50 kPa servocontrolled confining pressure was applied to the specimen during simulation. The particles of the simulated specimen cannot be too small because the computation effort for a specimen generated with small particles would be unrealistically large. On the other hand, shearing a sample consisting of too large particles can result in contact deficiency on the frictional sleeve during simulation. Thus, a trade-off between computation efficiency and numerical stability was made by enlarging the size of the Ottawa 20-30 sand by four times. The particle size distributions of the real and simulated specimens are compared in Figure 3. In the simulations, the sleeve was not scaled in an effort to avoid having a very large number of particles that would unrealistically increase computation time. In addition, since the shear stress on the sleeve is independent to the sleeve length, we shortened the sleeve and the specimen by 60 % to further reduce simulation time. Nevertheless, the assumption of spherical particles is too ideal since particle shape is ignored in conventional DEM models.
In answer to the shortcoming, a linear contact model with rolling resistance was implemented. The rolling resistance coefficient is a parameter that quantifies resisting moments in the counter directions of relative particle rotation. The details of the rolling resistance linear contact model can be found in the PFC 5.0 suite documentation or Ai et al. (2011) [10,11]. In the simulations, all the parameters were calibrated with the results of triaxial tests on Ottawa 20-30 sand. A summary of theses parameters is presented in Table 1. Normal-to-shear stiffness ratio (κ * ) 2

Interface shear stress
The results of axial and torsional axisymmetric shear tests with smooth and the textured sleeve mentioned above are presented in Figure 4. In addition, the curves for the 3D simulations with the textured sleeve are also compared with the laboratory results. In the figure, the interface strength is evaluated by the stress ratio calculated using the stress on the sleeve normalized with the lateral confining pressure applied to the specimen. Note that the shear stress on the sleeve was calculated by the summation of all the contact forces acting on the textured sleeve projected in the direction of central axis divided by the side area of the sleeve surface. For both shearing directions, much higher interface strength was observed on the surface of the textured sleeve than that of the smooth sleeve. The development of interface shear stress in both axial and torsional shear simulations agrees well with the results of the laboratory tests: The torsional shear loading mobilized higher peak and residual strength at the sleeve-particle interface than the loading in axial direction. However, in the torsional shear simulation, deep plunges of interface strength were observed when the sleeve is sheared for 10 mm and 22 mm. The difference between these two shear displacements (12 mm) is roughly equivalent to the distance between two consecutive columns of textures on the sleeve surface. The enlarged particle size was considered to be the cause of the strength degradation: when large particles roll over the texture elements, they fall into the "ditches" between diamond-shape textures resulting in the breakage of the strong force chains on the sleeve texturing elements.

Volume change
To verify our assumption, as shown in Figure 5, global volumetric changes of the specimens when sheared against the textured sleeve of Rmax = 1 mm were monitored during simulation. While shearing, the volumes of both samples increase for most of time. However, in the torsional shear test, an abrupt drop in specimen volume was observed when the sleeve displacement is in the range between 9 mm and 15 mm (highlighted in the figure). This phenomenon is in accordance with the recession of interface strength shown earlier in Figure 4. It suggests that as soon as the particles close to the sleeve rolled over the textures and fell into the "ditches", the constant confining pressure serves as the resistance to volume expansion-it drove peripheral particles towards the surface of the textured sleeve. The gaps created by the "rollover" effect near the sleeve were filled by the particles that were further away from the sleeve. As a result, the specimen shrunk. In contrast, the volumetric contraction of the specimen in the axial shear simulation was not as obvious. This is attributed to the configuration of the texturing patterns on the sleeve (Figure 1b): The untextured passthrough area provides flow paths around the asperities for the contacting particles and prevents them from clogging the texture when shearing soil specimens against the sleeve in the axial direction which is absent in the torsional direction.

Mesoscale response-Shear zone detection
To understand the mechanisms behind the difference in macroscale responses between axial and torsional shear tests, investigations at smaller scales are needed. Measurements of local void ratio have been used to characterize the shear zone, for the fact that volume expansion always comes with interface shear failure. Hence, the volume of specimen that undergoes obvious increase in void ratio defines the shear zone. Figure 6 illustrates the moving average technique adopted in the simulation to extract local void ratio in the specimen. The shape of the averaging zone resembles a hollow cylinder with a wall thickness of 5 mm. This is a logical volume to consider since the purpose is to obtain a relationship between local void ratio and the distance to the sleeve. The thickness of the moving average zones were determined by a trial and error process to balance the stability and precision of the measurements. A sequence of measurements as a function of distance to the sleeve surface was obtained by gradually expanding the radius of the cylindrical zone while keeping its thickness constant.  ratio measurements near the textured sleeve are prominently larger than the initial state e = 0.6. Represented by the vertical dashed lines, the boundaries of the shear zones are determined by the extents of the volume dilation zones. The results show that the shear zones in the laboratory tests (5 mm) are much smaller than the ones in the simulations (12 mm), which is also due to the effects of particle size scaling. To seek for a rule of thumb for estimating shear zone thickness, Hebeler et al. (2016) summarized the relationship between the shear zone thickness (normalized with D50) and mean particle size (displayed in Figure 8) based on laboratory soil shear tests, direct and axisymmetric interface shear tests and a review of existing literature [12]. A well-defined bi-linear relationship was observed in the plot: When particle size is smaller than 0.7 mm, the thickness of the normalized shear zone decreases monotonically. On the other hand, samples composed of larger particles tend to have similar shear zones in terms of the normalized thickness. The results of the laboratory experiments and numerical simulations performed by the authors are denoted by a red arrow and a star, respectively, on the plot and they were shown to be consistent with the bi-linear relationship.  Fig. 8. The bi-linear relationship between the shear zone thickness and the mean particle size (Modified from [12]). Figure 7 also shows void ratio reduction immediately outside the expansion zones in the laboratory experiments. The specimen in the torsional shear test behaved more contractive outside the shear zone. However, these findings were not observed in the numerical simulations, which is again attributed to the large particle size adopted in the simulations. The particle expansion makes the numerical models insensitive to local volume compression possibly because the contraction is very small compared to the particle size.

Particle trajectories
It is postulated that the different mesoscale responses shown in Figure 7 are due to the distinct patterns of force transfer at the particle-sleeve interface. By decomposing the force on the sleeve into friction and passive resistance components, Martinez (2015) proposed mechanisms of particle-sleeve interactions under axial or torsional shear loadings [13]. Figure 9 shows the mechanisms: In the axial shear test, load is mostly transferred in the vertical direction, resulting in particles in the shear zone being displaced along the shearing direction and causing only modest changes in void ratio outside the shear zone. In contrast, during torsional shear, the displacements of the particles close to the sleeve can be decomposed into radial direction and the direction of shearing. Fig. 9. The movements of particles in the shear zones of (a) axisymmetric axial and (b) torsional shear tests ( [13]).
These patterns of particle movements were validated in the DEM simulations by tracking the centroids of particles in the shear zones. Shown in Figure 10a, for the axial shear simulation, particles mainly migrated downward along the shear direction while these particles were restrained in a narrow band laterally. On the other hand, in conformity to the proposed mechanism, under torsional shear loading, the particles in the shear zone were not only displaced along the tangential direction of the rotation, but also in a centrifugal direction from the axis of the sleeve (Figure 10b). These particle movements show different interactions at interfaces and are considered to be the cause of the difference observed in macro-and meso-scale behaviour mentioned earlier.

Particle rotation
Another informative measurement at the microscale is the particle rotation, because it correlates with the rotation "frustration"-rolling incompatibility of contacting bodies-which is considered the main source of volume change under shear loading. The particle rotation was evaluated in terms of Euler angle which is used to describe the orientation of a rigid body with respect to global 3D coordinates. Figure 11 visualizes the variations of particle rotation in the planar cutting planes across the centre of either a column or a row of sleeve asperities depending on the shearing direction. Similar developments were observed in both simulations: At the beginning of the shearing phase, particles in a thin zone close to the sleeve rotated less than 6°. However, the torsional shear triggered more particles to rotate than the axial shear at the beginning of the simulations. When peak strength was established on the sleeves during shearing, the zones of intensive particle rolling expands further toward the boundaries of the chambers. Especially in the axial shear simulation, the rolling-intensive zone extends even ahead of the sleeve in the shearing direction. Finally, at the residual stage, the zone of intensive rotation (red area) was expanded to about 12 mm away from the textured sleeve in the lateral direction which is in line with the size of the shear zone detected in Figure 7. In the study of 2D DEM simulations, the correspondence between the zone of intensive particle rotation and shear zone was also reported. ( [14]).

Conclusions
This paper has shown that 3D DEM simulations are useful tools for revealing the links between the observed macroscale behaviour and smaller scale responses such as local void ratio, particle displacements, and particle rotations. The following conclusions can be made based on the numerical simulations: The proposed DEM model is able to replicate the key features of the macroscale response-shear stressdisplacement curve in both axial and torsional shear directions, although some fluctuations in the torsional shear simulation were observed because of the larger particle size adopted in the model.
In mesoscale, the distributions of local void ratio in the specimens after shear provide clear cut-off boundaries of the shear zones. However, the contraction zones immediately outside the shear zones observed in the laboratory tests were not observed which is also likely caused by the particle scaling.
Particle displacement trajectories helped to identify the difference of interfacial interactions when different shear motions were applied. It also agrees well with the previously proposed mechanisms of load transfer at sleeve-particle interfaces.
Particle rotation highly correlates to shear zone formation and volume dilation. The zone of intensive particle rotation near the sleeve has almost identical size of the shear zone.
Finally, simulations on specimens that consist of smaller particles are needed to achieve more accurate results, although more computation effort will be required to perform these simulations.