Image based simulation of one-dimensional compression tests on carbonate sand

High factors of safety and conservative methods are commonly used on foundation design on shelly carbonate soils. A better understanding of the behavior of this material is, thus, critical for more sustainable approaches for the design of a number of offshore structures and submarine pipelines. In particular, understanding the physical phenomena taking place at the microscale has the potential to spur the development of robust computational methods. In this study, a one-dimensional compression test was performed inside an X-ray scanner to obtain 3D images of the evolving internal structure of a shelly carbonate sand. A preliminary inspection of the images through five loading increments has shown that the grains rearrange under loading and in some cases cracks develop at the contacts. In order to replicate of the experiments in the numerical domain, the 3D image of the soil prior to loading was imported into a micro Finite Element (µFE) framework. This image-based modelling tool enables measurements of the contact force and stress map inside the grains while making use of the real microstructure of the soil. The potential of the µFE model to contribute insights into yield initiation within the grain is demonstrated here. This is of particular interest to better understand the breakage of shelly grains underpinning their highly compressive behavior.


Introduction
Carbonate soils cover over 40% of the world's seabed, where offshore structures, pipelines, artificial islands and other marine structures are founded. For the most part, carbonate soils are of biogenic origin comprising skeleton bodies and shells of small organisms, the shelly carbonate sands [1][2][3]. These soils are a complex and poorly understood material, in particular the high compressibility of the material is not well predicted using current modelling techniques. As a consequence, shelly sands have been placed into a niche classification of ''problematic soils'' in most design guides [4,5]. While failures are now relatively rare, conservative methods and high factors of safety are commonly used. Understanding the behavior of shelly carbonate sand is critical for the design of foundations for offshore structures. In particular, understanding the physical phenomena taking place at the microscale has the potential to spur the development of robust computational methods.
Shelly carbonate sands differ considerably from more commonly studied silica sands. An important distinction is the composition; carbonate soils are rich in calcium carbonate which has a much lower hardness value when compared to quartz. But composition plays only a part in explaining the less conventional behavior of shelly sands. A significant characteristic of shelly sands is the shape of the constituent grains, comprising a mixture of plate-like and hollow bulky shaped grains both characterized by high angularity values. Despite the variety of grain morphologies found in carbonate sands, typical shapes can often be associated to a finite number of species in a given sediment [6].
The significance of grain shape dependent behavior is widely reported in the literature. Santamarina and Cho in [7] have shown that while smooth grains are prone to contact slippage, angular shapes tend to intensify particle deformation at contacts. The shape contribution to the macro response derives from the arrangement of the grains and the grains interaction. Earlier studies have used thin sections to quantify the shape of shelly grains [8]. More recent work using X-ray micro computed tomography (lCT) images and advanced image processing tools [9] has enabled the quantification of the morphology of shelly carbonate sands in three-dimensions (3D) and paved the way to discrete modelling of the material.
This study makes use of a recent developed micro Finite Element (lFE) model [10,11]. The model takes as input the geometrical grain-scale data obtained from lCT to model the individual grains and their interactions in the framework of combined discretefinite-element method [12,13]. The soil fabric is therefore virtualized by meshing the constituent grains and allowing them to interact and deform according to appropriate constitutive model and frictional contact conditions. The contact response results from the deformation of contacting bodies, which accounts for the specificities of each contact surface. The rational is that using a true representation of the internal structure of the material and deformable grains will enable a more realistic representation of the physics of granular behavior.
A small sample of a shelly sand imaged during oedometer compression is used in this study to (1) obtain the microstructure to be translated into a numerical sample and (2) validate the numerical simulation by comparison against the experimental macroscopic behavior. The present paper comprises a brief description of the material and the experiments, followed by the numerical framework used, including the formulations for body motion, body deformation and contact algorithm. Particular emphasis is given to the numerical meshing of the shelly grains. The analysis and discussion of the results from the simulation are presented together with future directions of research.

Material and experiments
The material investigated was a shelly carbonate sand from the Persian Gulf with a median grain size of approximately 2.2 mm, the particle size distribution is presented in Fig. 1. This sand comprises a mixture of plate-like and hollow bulk shapes with high angularity features [9]. The enclosed voids within the grains were quantified using the intra-granular void ratio parameter (ratio between the volume of the internal void and the volume of solid) and a value of 0.095 was obtained. This corresponds to an intra-granular porosity of 8.67%. The tensile strength of these shelly grains was investigated by means of single compression tests [14] and it has been observed a wide range of normal force:displacement responses, which was attributed to the large variety of grain shape associated with biogenic origin of the constituent grains. The one-dimensional compression tests were performed under dry conditions inside an X-ray scanner. The system used was a Nikon XTH 225 ST located at the Research Centre at Harwell (RCaH), Oxfordshire (UK). The experimental set-up consisted of a load cell, a vertical piston a micrometer and a small oedometer, as shown in Fig. 2. The set-up was designed for this particular scanner and was mounted on the rotating table of the scanner. The size of the specimen was 13.5 mm in diameter with an aspect ratio (height/diameter) between 1.0 and 1.1. The force was exerted by a micrometer with 450 N axial loading capacity [15]. The sample container was made of Perspex with 2 mm thickness for which a value of less than 3 lm deflection under the maximum applied force was measured. The lateral friction has been minimized by considering a 1 mm gap between the container and the X-ray window. The exerted force was monitored by a low profile 'pancake type' load cell with a 500 N capacity.
The X-ray tomography data was acquired at an accelerating voltage of 90 keV. A total of 3142 projections were collected per scan, with an exposure of 500 ms per projection. The 3D images acquired had a resolution of 9.5 lm (length of voxel edge). In an X-ray scanning, the objects within the sample attenuate different levels of X-ray beam energy, depending on the material composition and density. Denser materials attenuate more than less dense materials and this difference in attenuation is represented by the intensity values of the voxels (3D pixels). The contrast of intensity level allows for differentiation of the features within the image. The scanned images have the dimensions of 1536 9 1536 9 1600 voxels. The large variety of sizes and shapes and the contrasting colors defining different grains can be seen in the top and side views cut through the experimental sample, shown in Fig. 3.
The experimental test was performed at five loading stages up to a maximum 10% axial strain. In-house imaging processing codes were employed to segment the images in order to identify the individual grains. This includes firstly to binaries the images using Otsu's thresholding and subsequently apply an iterative watershed algorithm to overcome the challenges posed by the large diversity and complexity of the shapes associated with the biogenic nature of shelly sands. The codes are fully described by Kong and Fonseca [9].

The numerical framework
In the framework of combined finite-discrete-element method, grain deformability can be described by a finite-element method formulation, whereas the

Body motion
The explicit dynamic analysis is established on the implementation of explicit integration and the use of diagonal mass matrices. The acceleration (€ u) is calculated based on Newton's second law. The equations of motion are integrated throughout time (t) using an explicit central difference integration rule: where _ u is velocity, u is displacement and i refers to increment number.
An advantage of using explicit time integration is the possibility of utilizing the diagonal lumped mass matrix. Computational efficiency can be improved by using the inversion of the mass matrix, for which the computation for the accelerations at the beginning of the increment can be reduced to a simple operation: where M is the diagonal lumped mass matrix, F is the applied load vector, and F I is the internal force vector. The explicit procedure requires no iterations and no tangent stiffness matrix.
In the explicit scheme, the time step must be small enough to ensure the stability of the integration. Abaqus finite-element software [16] used here automatically adjusts the time increment during the analysis based on a global estimation method. The advantages of using a global time increment estimation is the constant update of the maximum frequency of the algorithm leading to a better and more stable simulation. The trial stable time increment is calculated for each element in the mesh as follows: where x element max is the maximum eigenvalue of the element. A conservative estimation of the stable time increment is given by the minimum value taken over all elements.

Body deformation
Deformability results from the external strains that the particulate material is experiencing. If there is no strain, the particle flows like a rigid body. The motion of the element nodes is governed by internal forces acting on them. The nodal forces include the contribution of the force transmitted through contact interaction, deformation of a discrete element and the external loads, as described in the following formulation: where x is the nodal displacement vector, F int is the internal resisting force vector, F ext is the applied external load vector and F c is the contact force vector. The computational contact algorithm requires specifying the geometry and the constitutive relations for the contact interactions. A general contact algorithm from Abaqus/Explicit is used. This relies on a robust, yet computationally efficient search tools for tracking the motion of the grains. These sophisticated algorithms enable fully automated contact detection between all grains at each time steps. Grain-to-grain interaction during grain rearrangement is optimised by allowing for the entire surface of the deformable body to become a potential contact. This presents a clear improvement when compared with previous studies on single grain simulations, for which the contact area is specified to be limited to a small region. The general contact algorithm generates contact forces to resist grain penetration and uses a finite-sliding formulation that allows for arbitrary separation, sliding, and rotation of the surfaces in contact. The finite-sliding formulation assumes that the incremental relative tangential motion between surfaces does not significantly exceed the dimensions of the master surface facets, but there is no limit to the overall relative motion between surfaces (more details can be found in section 5.1.2 of Abaqus manual). A surface-to-surface contact discretisation is used. In this discretisation each contact constraint is formulated based on an integral over the region surrounding the contact node (Fig. 4a). In contrast with the commonly used node-to-surface approach (Fig. 4b), for the surface-to-surface discretization (Fig. 4c), the continuity constraints are not enforced only at the discrete nodal points but along the entire boundary. Using a surface-to-surface formulation, the shape of both the slave and master surfaces in the region of the contact constraint is considered, and this is critical for grains with irregular shape such as real sand grains. For these reasons, surface-to-surface approaches provide, generally, more accurate stress measurement when compared to node-to surface approaches.
The contacts are modelled with properties of hard contact in the normal direction and Coulomb friction in the tangential direction. In hard contact, when two particles are in contact, i.e. the overclosure value equals zero, all the pressure is transmitted from master to slave surface. This is shown in the following expressions: where p is pressure and h is overlap. The contact constraint is enforced with a Lagrange multiplier representing the contact pressure. The virtual work contribution is:  Fig. 4 a Schematics of the contact domain between two contacting bodies X 1 and X 2 ; b node-to-surface approach, c surface-to-surface discretization Matlab [17] script to generate the image-based mesh. [10]. Triangular iso-surfaces were extracted from the 3D image with pre-set values for density, which controls the size and number of triangles representing the surface of each grain. Figure 5 shows two grains represented using different mesh sizes. The mesh size and quality is controlled by four parameters, which includes the minimum angle of the triangles, the maximum volume of an element, the maximum surface element size and the maximum distance between the center of surface bounding circle and the center of the element bounding sphere. Figure 5a and b show a coarse mesh with parameter values of 40°, 30°, 30°and 10°, respectively. It can be observed that this mesh size cannot capture well the irregular features of the shelly grain. The mesh was progressively refined in the other two examples presented and a minimum angle of 40°was used in all cases due to efficiency issues. The four parameters used in the two other cases were: 40°, 20°, 30°and 5°(intermediate mesh, Fig. 5c, d), 40°, 30°, 30°and 1°(fine mesh, Fig. 5e, f). It was found that the fine mesh is able to capture well the grain outline while avoiding the more expensive computation cost of using a more refined mesh and was then chosen for the simulation. It can be seen that for computational efficiency larger triangles are used to represent flat areas in the grain surface whereas more angular features are described with smaller triangles. The key advantage of this meshing technique is its ability to preserve the original boundary of the grain with no restrictions for complex topologies or intragranular voids. Figure 6 shows the meshed full sample, the element mesh size that is, on average, five times the voxel size of the original image.

Simulations
The numerical mesh was imported into the lFE model that was implemented in Abaqus with an explicit algorithm that uses a dynamic framework. The sample was confined laterally using a rigid cylinder and the displacement was applied at the top of the sample. Elastic behavior was assigned to the model. Average values of 25 GPa for Young's modulus, 0.3 for Poisson's ratio and 0.3 for the coefficient of friction were used. These parameters were calibrated using experimental data from single grain compression tests and curve fitting the force-displacement numerical results (Fig. 7). The experimental single grain tests were carried out using the Instron 5969 straincontrolled machine (Instron, Norwood, Massachusetts). The instrumentation accuracy was measured to be \ 1 lm for displacement and \ 0.1 N for load. The numerical grain was compressed between two rigid plates to reproduce the experimental conditions. Another input parameter required for the numerical simulations was the density of the material, for which a value of 2.5 tonne/m 3 was used. The analysis of the full sample took approximately 3 days (running on DELL Precision T7610 with 12 dual cores 3.2 GHz CPUs). Figure 8 shows the stress distribution, computed using the von Mises criterion, in the full assembly at five different strain values. It can be seen that at the initial stage (e a = 0%), since the sample is not loaded yet, all grains have a stress field constant and equal to zero (Fig. 8a). At e a = 5%, the stress concentration starts at the grain contacts and propagates through the grains before being transmitted to other neighboring grains, again by way of their contacts (Fig. 8b). At e a = 7.5% and e a = 10% (Fig. 8c, d), the assembly becomes more heavily loaded and this is shown by the larger internal stress values exhibited at the grains.
The stress at the contacts can be better visualized in the vertical cut through the sample presented in Fig. 8e. These data enable the identification of the stresstransmitting particles and the investigation of the micro-mechanisms that lead to the breakage initiation within the grain 6 Results and discussion The stress:strain curves obtained from both experiments and numerical simulation are shown in Fig. 9.
Overall, we can say that there is a good agreement between the two curves. The numerical samples shows a slightly softer response in the initial part of the test (up to 3% strain), which seems to suggest stronger grain rearrangement of the numerical grains when compared with the experiment. Although it is important to note that given the large size of the grains, even small rearrangements can have a significant effect on the overall settlement of the sample, the difference in the responses may also suggest that the coefficient of friction needs to be adjusted. In fact, because shelly carbonate sands are formed of different organisms/ species the numerical parameters to be used will have to come from tests on representative species and assign to the numerical grains according to their grain type instead of using a global value.
To investigate the reliability of the numerical model in predicting the stress distribution within the individual grains, the data from the single grain compression simulations were used. In this case, it is known that the stress will start increasing at the two contacts with the top and bottom platens and propagate from there. This can be confirmed by observation of the grains shown in Fig. 10.

Conclusions
This paper presents a preliminary investigation into the numerical simulation of shelly carbonate sands that takes into account the detailed representation of grain shape and intragranular voids. The capability of the model to measure the stress concentration at the contacts and its propagation throughout the grain skeleton for materials with complex grain morphology is demonstrated. The macro response of the simulation shows generally a good agreement and gives confidence on the robustness of the methodology to simulate complex morphologies numerically. The results presented here appear to suggest that for the case of sands of biogenic origin with a diversity of species forming the granular system, more accuracy can be obtained by assigning the model parameters according to the grain type instead of using a global value.
Future work will include a more in-depth physical and mechanical characterization of the grains to refine the numerical representation of the material and carry out a micro-scale comparison between the experiments and the simulation. Finally, the enhancement of the lFE model to include simulation of breakage, taking advantage of the yield initiation within the grains, will constitute a step change for the numerical modelling of shelly carbonate sands.