{'fraction_non_alphanumeric': 0.02956065810884108, 'fraction_numerical': 0.010594829144820106, 'mean_word_length': 4.619894330420646, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High compression of granular assemblies of brittle hollow tubular particles

This paper is devoted to the micro-mechanical origins of the high compressibility of brittle tubular particle assemblies. The material is extremely porous due to the presence of a large hole within the tube-shaped particle. The release of the inner void, protected by a fragile shell, gives the material a very strong ability to compress. The compressive response is investigated by means of the Discrete Element Method, DEM, using crushable-elements. To address the complexity of the model, a step-by-step break-down is developed. The paper comprises the comparison of the numerical results with both results obtained by the authors and existing experiments. With the insights provided by the DEM, we have sought to better understand the phenomena that originate at the grain scale, and that govern macroscopic behaviour. Grain breakage was proven to control the compressive behaviour, and thus, the importance of internal pores dominates the inter-particle voids. Then, a novel concept of compressibility analysis has been proposed using the separation of the double porosity and the quantification of the pore collapse through primary grain breakage. Finally, a general, geometrical development of a semi-analytical model has been proposed aiming the prediction of the evolution of double porosity vs axial strain.


Introduction
The research presented in the current paper originates from an engineering concept of the mono-block compressible archsegment (US Patent pending on "Voussoir Monobloc Compressible", VMC, developed jointly by CMC and Andra [1]). VMC is an innovative technique of fabricating the tunnel lining. This pre-casting technology combines a conventional reinforced concrete segment with a compressible layer on its outer side (extrados) [2,3,4,5]. In the engineering practice, there exist some cases in which the host rock around tunnel was partially replaced with a compressible material that prevents an overloading of already existing tunnels [6]. Such approach treats the tunnel and the compressible material as separate systems, whereas in Andra's double-layer lining the two systems cooperate closely. Another innovation of the VMC design concerns the compressible layer itself. Instead of a highly compressible material, like foam, a cemented assembly of highly porous and quasi-fragile particles was envisioned. From a mechanical point of view, the influence of the weak cement bonds between the particles is essentially dominated by the characteristics of the granular core. The key benefit of introducing the granular material is related to the charge transfer that it enables [7]: a local load can be minimised by spreading it on the extrados surface of the reinforced concrete segment. In addition, under sufficient load, the layer will experience some compressible deformations allowing for a stress dissipation ( Figure 1). As long as the outer layer exhibits its compressible capacities, it limits the load transferred on the inner layer. The quasi-fragile particles were manufactured from the excavated COx clay stone, first dried then baked at high temperature, and therefore have brittle properties. The compressible response (Figure 1) is predominantly achieved by closing the pores. The large internal porosity contributes significantly to a highly compressible mechanical response, e.g., [8]. This material owes its high porosity to the tubular geometry of the constituent particle, called shell (Figure 2a). The geometry was specially selected for this application, such that the shellscale porosity, i.e., volume of internal void to total volume of particle, reaches 50%.
First purpose of this work was to design a numerical model capable of simulating shell breakage. For this, methods based on peridynamics, e.g. [9], or the material points method, e.g. [10], would be suitable but expensive in computation time. Another approach like Discrete Element Method, DEM, as it was introduced in the 80s [11], is also particularly well suited for breackage granular materials. A couple of DEM techniques can include the capability of the particle to break: fragment substitution, FS, e.g., [12,13,14], and Bonded Particle, BP, e.g., [15,16,17]. In FS methods, the breakage criterion is defined as a function of the critical load that particle cannot withstand. Once broken, the original particle is replaced by a set of new particles. The replacement mode is usually structured under the condition of a self-similarity. In other words, the replacement particles keep the same shape but are smaller in size. Self-similarity suits a material like sand but is hardly applicable in the case of tubular shells. It is then more convenient to use BP approach that connects sub-particles by means of the cohesive links. The breakage is verified between sub-particles, thus it requires knowledge of a material strength. A rigid clump of elements can in detail reproduce highly complex shapes [18]. In the field of geomechanics, this technique was found considerably useful in the investigations of the particle shape, for example aiming the concavity effect [19,20]. This work also took advantage of clumping method when constructing a non-spherical polygonal-like sector. To generate the tube-shaped geometry (Figure 2b), the sectors were linked together using BP method. Section 2 elaborates on how the BP Model was adapted to the case of the cylindrical shells. The paper is then continued with shell-scale investigations, both numerical and experimental, in Section 3. This part is dedicated to the identification of discrete parameters and a calibration of the breakage model. Section 4 is a tribute to the true packing with the focus being redirected to the internal state variables. It mostly discusses shell orientation: the theoretical approach and the results of an experimental characterisation with X-ray acquisitions. Sensitivity analysis presented in Section 5 gives a better understanding of how the discrete parameters as well as the initial state affect the macroscopic behaviour to oedometric compression. At this point, the complexity of the model, the necessary compromises and the validity point will be also discussed. This section also begins the discussion over high compressibility exposed by the stress-strain curve (with a reference to the experiment). Yet, the analysis of compressibility is fully completed in Section 6. To finish, a geometrical and analytical development of a prediction model incorporating the breakage phenomenon is proposed. Finally, a brief summary is given at a very end, in Section 7.

Model of crushable shell
The experimental study of the behavior under oedometric compression of tube shaped clay-shell assemblies shows that before being completely crushed and reduced to powder (ε a = 0.9 in Figure 1), the particles break essentially radially, up to ε a 0.6, insets in Figure 1.
To simulate the breakage of the shell (Figure 2a), this model was based upon the Bonded Particle (BP) approacha DEM technique clustering smaller, rigid sub-particles into a "parent" particle. Therefore, in the current paper a numerical, tube shaped, breakable particle is referred to as the cluster, while term sector stands for its constituent ( Figure 2b). As presented in Figure 3, each sector has been generated as a polyhedron-like shape by clumping 3D solid figures (sphere, cylinder and rectangular prism) rigidly 1 together,and thus, preventing it from any deformations. Sectors of the same parent cluster are joined by means of the distinctive connections, hereinafter called links. For the sake of simplicity, the link is only located in-between spheres that are positioned at the sector "corners" (Figure 3). Relative movements evolve in all four links aligned on an interface between two sectors, and thus, those links develop the elastic forces as follows: where k • are the link stiffnesses, and δ • are the relative displacements between two spheres at there touching point. In other words, δ • can be considered as shortening or elongation of the local links in corresponding direction. The subscripts I and II refer to the fracture mode, but can also be seen as the normal and tangential direction of the link. To avoid pure mode-I fracture ( f II = 0 N), f I must not exceed a tensile yield threshold − f I (Figure 4a). The material also withstands shear of the pure mode-II type ( f I = 1 lack of relative movement between the solid figures 0 N), as long as the force f II remains in the ± f II range (Figure 4b). The yielding surface for any mixed loading (combined mode-I and mode-II) in the f I -f II space, is defined by the following function: where q enables to adjust the surface between a bi-linear shape (q = 0) and a semi-infinite rectangle shape when q tends towards infinity. Figure 4c shows an example of yield surface for q in the order of 2. All four links are broken as soon as in one of them the breakage criterion is met: Once the link is yielded, it is irreversibly broken. At this stage, frictional contact interactions become possible between the two sectors that have just been detached.
In case of frictional contact, the force-displacement relations used are quite classical, and can read as and where f n is a normal repulsive force proportional to the normal distance δ n in case of overlap, and the increment of tangential force ∆ f t is proportional to the increment of sliding distance (that is the sliding velocityδ t multiplied by the time step ∆t) as long as f t remains below µ f n . As a result, a contact involves only compressive f n as demonstrated in Figure 4a. The tangential force is constrained to remain within the range ±µ f n ; Figure 4b. This leads to a bi-linear yield function (Coulomb cone) ruled by a given friction coefficient µ; Figure 4c. In discrete element simulations of quasistatic regimes [22], energy dissipation is a major concern [23]. In our situation, the fracture of the brittle shells is characterised by high energy releases, so the efficiency of energy dissipation becomes even more critical. To handle it, two means are employed. The first is a viscous force added to the contact normal force, such that the damping is at 95% of its critical state corresponding to perfectly inelastic collisions. The second means of damping is not strictly physical but rather numerical. It is the non-local damping which has been set to 30% as prescribed by their designers [24]. The purpose of its use is to minimize the surges of kinetic energy induced by sudden link failures. These two damping parameters were kept identical for all simulations, and we obviously ascertained that they did not significantly affect the macroscopic mechanical responses. The numerical part of the work has been conducted using a parallelised tool named Rockable 2 which algorithm operates on complex, sphero-polyhedral shape features (here they are sectors). Regarding the force laws, one needs to specify 8 parameters: ¥ 5 for the links (ruling the breakage): k I , k II , f I , f II and q, ¥ 3 related to the frictional contacts (between detached sectors or clusters): k n , k t and µ.
In the following Section 3, their determination and/or tuning will be discussed. Notice that the adjustment of some parameters (k n , k t and µ) require investigations on assemblies so that clusters can interact. This aspect is to be discussed further in the section 5.

Assessment of shell mechanical properties
The model requires a calibration. This process of parameter tuning, which is essential to the reliability of the model, is not a straightforward task of bridging quantities at the micro and macro scales. We are focusing here on the breakage of a single shell. The strength of the constituent material of the shell must be considered in the evaluation of the threshold forces of Equation (2). The present modelling is such that the parameters f I and f II monitor the resistance of the shell 2 In house numerical code developed by V. Richefeu , in collaboration with the authors of this paper in combination with the chosen sectioning (i.e., the way the cluster is partitioned). In other word, the threshold forces are structural parameters and they are in no way constitutive parameters. Therefore, experiences of breaking a single shell have been relied on. A series of strain-controlled Uniaxial Radial Compressions, URC, has been performed to initiate a tensile failure. The loading conditions bear a strong resemblance to Brazilian compression test [25]. A ring-cross section of shell was diametrically compressed by the force F I in which direction the ring deformation ∆ d was measured ( Figure 5). A detailed description of the experimental campaign including the characterisation of tested shells can be found in [26]. Figure 5 shows the variety of the mechanical responses captured in the experimental campaign (green curves), but only the average response in terms of maximum force and displacement (F I 122.5 ± 19.4 N, ∆ d 1.83 × 10 4 ± 0.35 × 10 4 m) was numerically adjusted (red line in the plot). Since the force law in the links are basically linear, (see Figure 4a), the curve F I ↔ ∆ d also obeys the linear relationship with the constant slope K I -a shell radial-stiffness that depends on the ring geometry (size, relative thickness, and so forth). Note that also the experimental curves convert to linear trends after a contact stabilisation manifested by the non-linearity of the curve. This stabilisation is caused, in particular, by a surface chipping and/or a slight movement of shell that inhibit the increase of force. The numerical investigations have shown how to model the shell response to URC throughout two micro-parameters: the yield tensile threshold f I and mode-I stiffness k I . The average response presented in Fig [29] of F I can be easily converted to its numerical substitute shown in Figure 7. It characterises the probability of shell survival from the tensile failure with the Weibullian cumulative density function (cdf ) for scale and shape parameters, respectively x 0 and m [26,30]. Using this parameters to distribute f I on the links in a simulation, the heterogeneity of material, partly caused by geometrical imperfections of shell (with respect to idealised shape of the cluster, Figure 2b), can be introduced.
Finally, we could verify the shape partitioning in view of the breakage context. In case of URC, the cluster model with only circumferential slicing is sufficient to reproduce the experimental breakage into 4 fragments: in the vertical plane recognised as primary breakage, and in the horizontal plane as secondary breakage (see the insets of Figure 5). However, one must have noticed a distinction in the occurrence of secondary cracks -∆ d exceeds the extents observed in the experiments. This is mainly an effect of the model assumptions. More precisely, this modelling did not distinguish the initiation and the propagation of the crack as it arises ex-perimentally -cracks occurred in just one time step which led to a sudden, high release of energy stored in link (manifested by the drop of the force). After the force mounted up again but in the model we deal with half-rings so the cluster stiffness K I decreased due to change of geometry.
As seen before, the shells are subdivided into N circ sectors in its circumference. The sectors are the unbreakable elementary entities of the modelling. In the axial direction, sectors can also be divided into N axial parts, Figure 6. For oedometer test simulations at low stress levels, that is when the shells tend to fail due to cracks initiated in tension, this shape partition respects the experimentally observed fracture pattern [31]. But the higher is the stress level, the more deviate the longitudinal sectors from the actual shape of fragments 3 . In this work, shell clusters have been sliced by using N axial = 1 and N circ = 12, for most of the cases. In this way, the compression rate has been consciously limited in the model. Nevertheless, with the prospect of larger scale modelling, reducing the amount of particles in favor of lower computational cost was a necessary initiative. Local parameters must be adjusted when the slicing of the shell shape is changed such that: ¥ k I increases proportionally to N circ with 24 sectors as a reference ¥ f I changes with N axial referencing no axial division -4 links to break in the axial direction In-plane shear tests were performed using an in-house shear device inspired by the commonly used direct shear box such that shell breaks in the middle creating 2 parts; Figure 8. The shell is first caged in two-part shear box and then sheared by imposing a constant relative displacement rate. The tangential force was recorded during the test. Unfortunately, these attempts of a pure in-plane shear test were hardly satisfying as the shells mainly failed in axial compression for a where m is the Weibull's modulus, x 0 is the scale parameter and x = f I . Red points mark data deduced from the experimental campaign, while the curve is a theoretical trend, that was found from fitting equation force F II 2 000 N. As a consequence, a minimum tangential force of 2 000 N was assumed to be needed for mode-II failure. Fig. 8 The concept of in-plane shear failure. A shearing force F II triggers failure through two axial cracks, a breakage manner directly applied to the radial slicing shown in Figure 2b.

F II F II
In the current model, four links are set between subsequent sectors. The critical shear strength f II of each link is thus set to 250 N so that the triggers shell-failure in mode II, has illustrated in Figure 8 is F II = 2 000 N. For sake of simplicity, hereafter the link stiffness k II was set equal to k I .

Combined mode-I and mode-II failure (q)
The rupture of a link can result from pure mode-I tension or pure mode-II shear. In practice however, the loading state of a link is actually a combination of these two modes, and it is necessary that a failure criterion takes both loading modes into consideration. This is the reason why the yielding function (Equation (2)) was introduced (see Figure 4c). The yield envelop is cone-shaped when q = 1, which is similar to a linear Mohr-Coulomb surface with the major difference that the link response remains elastic and linear inside the surface. Increasing q, the yield surface tends to become a prismatic shape, and the combined influence of the two pure rupture tresholds tends to vanish (they become completely independent for q = ∞).
A number of URC and in-plane shear simulations, carried out on single shells, have shown a limited influence of q on failure provided that q ≥ 2. Finally, a parabolic yield function was used (q = 2), following [21].
3.4 Sector-to-sector frictional contact parameters (k n , k t , µ) Frictional contact can act between sectors of different clusters, or between sectors of the same cluster if their common adhesion plan has been broken. One can estimate the contact stiffnesses as follows. For Young's modulus E = 4 Ga and the Poisson coefficient ν = 0.29, the dimensionless stiffness parameter for 3D Hertz contact expresses [32] where σ 0 is a reference mean stress. For a typical shell size being the shell diameter d (Figure 2), the normal contact stiffness is k n = 5.34 × 10 6 N/m, in the case of a linear normal contact law, κ = k n /(dσ 0 ) [32]. Recall that we have always been keeping the ratio k n /k t = 1. Then, in case of the cluster with N circ = 12, it was possible to equal the link and frictional contact stiffnesses k n = k t = k I = k II at the same value of around 5.5 × 10 6 N/m. Concerning the contact friction coefficient, it was estimated using an experimental device of 3 shells described in [35] and [26]. The following value was used: µ = 0.24 ± 0.06.
There exists other non-numerical parameters that one might find influential for the shell assemblies, e.g., initial density, contact network, orientation of shell. Such investigations cannot be limited to one cluster, and therefore, Section 5 is dedicated to the parametric study at the sample scale, but first a realistic numerical sample needs to be generated. Thus, in the following Section 4, we present the experimental identifications of some parameters that serve us as the starting points.

Microscale characteristics of a real shell assembly
Let us start with the internal variables that on various complexity level refer to the fabric of granular material. Probably, the most basic measure is a density. This study relay on the number density n measuring the amount of shells in the volume unit. The true values that the modelling aims are n = 151 384 ± 1 319 m −3 for the sample with 2.09 m 3 volume and n = 155 129±3 952 m −3 for the volume of 6.43 m 3 . The given values concern the samples being a shell assembly, a situation we will address numerically (Section 5). Towards more comprehensive characterisation, a couple of cylindrical samples with diameter 12 cm and height of about the same size were extracted from the compressible layer of the tunnel segment, VMC. On that account they are the shell-mortar composites. Thanks to X-ray radiographs acquired in RX-Solutions scanner of Laboratoire 3SR, it was possible to reconstruct their volumes (Figure 9a). The acquisition capabilities of the RX-Solutions scanner are specified in Table 1, but the precise descriptions of the X-ray imagining technique and the scanner itself are not the object of the current paper, yet can be easily found, for example, in [33].

Shell detection in 3D images
Form a binary image, that is a 3D solid-void matrix, we can extract a couple of datasets. A spacial detection of the shells motivated us to perform the X-ray scanning in the first place. The full detection concerns two components: the position and the orientation. 3DShellFinder algorithm was specifically designed to suit the tubular geometry of the shell. A perfect tube is pre-specified such that it can be inscribed in the true body of the shell. This tube is then represented by a set of the search points as shown in the inset of Figure 9b. For each targeted shell, the algorithm probes different positions and orientations varying them by a considerably large amount of the small increments. An attempt to detect shell firstly consists in assigning the search points to the voxels in the image. Their grey level equals either 1 for voids or 0 for solid. Then, the algorithm proceed with the minimisation of an error function: where I is the 3D-image and N ZOI is the number of search points that are located within the zone of interest (ZOI). The function of Equation (7) is parametrised by the position (vector x x x shell ) and the orientation (quaternionq q q shell ) of the search points P i . For the search points perfectly aligned with the shell in the image E (x x x shell ,q q q shell ) = 0, but to increase the probability of the successful search the algorithm accepts a small error, e.g., in the following case E (x x x shell ,q q q shell ) ≤ 0.02. For the sample of Figure 9a, the search tube was circumferentially discretised into 18 points for two rings of radii 70 and 80 voxels. Axially, 30 points were densely distributed along the length of 120 voxels. Figure 9b shows 229 shells (green points) successfully located in the image. We have estimated 4 that the algorithm identified around 98 % of shells (positions and 3D orientations). The unsuccessful searches were caused by (i) the damaged of the boundary shells or (ii) a strong imperfection of the full shell within the sample core.

An approach to analyse shell orientations
More complex aspect of studying the sample fabric should include a combined treatment of the particle orientation and the contact network [34]. Experimentally, with the X-ray radiographs we were able to address only the first one. Considering the axisymmetric shell, the ring cross-section defines the weak plane, while a central axis indicates its strong direction. It is worth mentioning that in the global coordinate system the vertical axis corresponds to the loading direction. Hence, an investigation of the shell/cluster orientation can be limited to the inclination of the central shell axis with respect to the vertical. This inclination can be measured by an angle α, but herein it is quantified indirectly with its cosine, such that | cos α| = 1 for vertically resting shells and | cos α| = 0 for the horizontal orientation as in URC tests. Figure 10 presents the probability density function of the shell orientation with a strong but noncontinuous anisotropy proven by a well marked peak for | cos α| = 0. Despite a natural tendency of the shell to embed horizontally, there still exists a partial isotropy for the remaining orientations. As a consequence, in oedometer test the shells are mainly radially compressed, and one can expect mainly the tensile stress to cause the primary breakage.
Additionally, a void-to-solid ratio e could be easily extracted the from 3D image, where solid includes both the shells and the mortar. For the sample shown in Figure 9a the standard void ratio e equals to 1.283 for both shells and mortar. The calculation disregards the boundaries of sample that were spoiled due to a extraction process used to drill the specimen from the tunnel segment. A rough estimation of e for shells (without taking into account the mortar) only reaches 2.4. A difference between those values will play a crucial role in the evaluation of the shell-mortar modelling where the volume mortar will not be represented (not further discussed in current paper).

Sensitivity to the model inputs
As seen in the Section 3, the numerical model requires numerical, mechanical and geometrical parameters. In order to carry out this study, it was fundamental to first perform a sensitivity analysis and to try to really understand the effects of the parameters. The oedometer test consists in imposing an uni-axial compression while the boundaries are fixed in remaining directions ( Figure 1). The mechanical response is displayed as the relationship between the axial stress σ a and the natural axial strain: where H and H 0 are the current and initial high of the sample, respectively. The current analysis is essentially a bare judgement of the consequences of the parameter variations on the macroscopic mechanical response. The present section attempts to make clear which are the most influential parameters and how through their value control the mechanical response.
To generate a specimen of packed clusters, a deposit under normal gravity have been simulated. Initially, the clusters were placed, without contact, within a cylindrical container closed by a flat lid on top; the threshold forces of the bonds that hold the cluster segments were set so that the clusters can not break in this phase. Each cluster has been not only positioned in a separate node of a cylindrical grid but also randomly oriented. Afterwards, the suspension deposited due to the activation of gravity. Additionally, an off-grid movement has been imposed at each cluster. Such numerical trick is equivalent to the sample shaking which prevents formation of local cavities, especially within highly frictional assemblies. Once the assembly embedded at the bottom of the container, the energy dissipation process needed time to reach a satisfactory equilibrium state.
While falling, the clusters collide through frictional contacts (see Figure 4b, dashed line). The greater is the friction coefficient, the looser is the cluster packing at the end. With series of numerical deposits, an analytical relationship between the inter-granular friction coefficient µ shell and the density has been determined [26]. Therefore, it was possible to tune the initial density. The deposit was simulated with µ shell set to 0.08, with friction-less walls (µ wall = 0). The sample density is then in the order of the experimental one (see section 4).
In the following, the sensitivity to different parameters is discussed in the context of oedometric compressions of cluster samples, at constant loading rate. The loading rate is chosen low enough to ensure a quasi-static regime. To do so, we ensured that the inertial number I =ε a m/σ was less than 10 −3 [27,28]. Additionnaly, the time step was chosen to ensure the numerical stability of the integration scheme for the equations of motion. Figure 11 shows stress-strain trends obtained numerically for a dense sample for selected couple of friction coefficients, µ shell or µ wall . No clear influence of any friction coefficient can be noticed. It can however be noticed that for a given strain, the vertical stress increases with the friction coefficients. At the initial state, before loading, nearly 89% of the interactions are bonded and 11% are frictional contacts. During the compression, some links break and the later percentage increases. However, an oedometer compression, unlike classical shear loading, does not massively mobilize friction between the constituent rigid elements. This can explain the quite minor influence of the frictional parameters on the compression curves.

Parameters of non-bonded frictional contacts
The other parameters of non-bonded contacts are the stiffnesses. To avoid unreasonable overlaps, k n and k t were varied in the order of 10 6 to 10 7 N/m. In considered range, their values were not observed to significantly effect the macroscopic response.

Cohesive links
Analogous arguments must be made concerning the stiffness of the link, both in mode-I and mode-II direction. Particularly, k I was well estimated from the experiments, and it does not seem meaningful to study the sensitivity to this parameter. The key to the understanding of the mechanical response to oedometric compression is the strength of clusters, controlled through the yield thresholds (see Section 2), f I and f II . Whereas tensile yield limit f I has been kept constant over simulations the shear one f II was varied, f I / f II = 0.34, 1.00, 1.70. Figure 12 makes clear that the condition f II ≥ f I is sufficient to disable the influence of f II . For f I / f II = 1.70, a remarkable influence of the shear threshold exists. It results in a decrease of the plateau and a stiffening appearing at a lower axial strain. As discussed in Section 3, the shells exhibits strong variability of the tensile strength characterised by the weibullian distribution; Figure 7. In such a distribution [29], the probability density function is described as a combination of power law and exponential function:  It is ruled by the modulus or shape parameter m, and the scale parameter x 0 . A modulus m < 10 corresponds to a wide distribution, i.e., a strong variability. Here, instead of using a constant value of f I , we choose to take into account a slight variability whose statistical distribution follows a Weibull distribution whose parameters are determined experimentally: m = 7.2, x 0 = 92.4 N (see Figure 7). Since f II = f I /0.34 for each cluster the variability of cluster strength apply both on f I and f II . Figure 13 presents the macroscopic stress-strain response with and without weibullian variability and explores the increase of x 0 (that is an assembly of stronger and stronger shells having the same breakage variability). The modelling with no variability employs only the average strength ( f I = 85 N and f II = 250 N). It is shown that the variability has a minor influence on the mechanical behaviour. The sample structure already provides a diversity in local fabrics and thus on the individual cluster failures. Thus, we state that an homogeneous cluster strength is sufficient in our simulations ( f II and f I have constant values). On contrary, the level of this strength plays a key role. The higher is the tensile yield threshold f I (also represented by x 0 ), the higher axial stress withstands the sample.

Initial density of the packing
Initially, the numerical campaign focused on dense assemblies whose initial density was of the same order of magnitude as the experimental measurements. Here, it is pro- posed to account for the influence of the initial density of the samples on the oedometer test response; the density being controlled by the number of particles per unit volume, denoted n. Three densities were tested. The mechanical response curves are given in the Figure 14.
For very small strains, ε a < 0.01, the three curves appear superimposed, revealing equivalent initial stiffness. Once the stress reaches 0.01 MPa, the 3 curves differ, with the loos-est specimen (black curve on the Figure 14) showing a slow hardening phase while the other two denser specimens see their strengths take off towards a peak at 0.3 MPa for the red curve (intermediate n) and almost at 0.7 MPa for the densest specimen (blue curve). This result is finally quite classical. It is comparable to what can be observed on standard granular materials: the higher the density, the higher the average number of contacts per grain and the higher the strength of the material for low strains. For 0.05 < ε a < 0.5, the intermediate σ a -plateau has the highest slope for the lowest cluster density specimen, with all clusters being broken around ε a = 0.5. For the two other densest specimens, the σ a -plateaus have similar slopes and end for ε a = 0.4. Once this intermediate regime is over, the hardening phase begins ("linear" relationship in a semi-logarithmic scale). The initial density seems to play only a secondary role here, all the clusters having been broken.
On one side, Figure 14 makes clear that the lower is the initial density, the larger is the final strain range. On the other side, the stress level is significantly lower despite the identical cluster strength.
Against our expectations, the loosest specimen exhibits the oedometer behaviour closest to the experiments whereas it has a high density of shells. When the DEM specimens have an initial cluster density similar to that observed experimentally with shells, the oedometer behaviour differs from the experiments in that the initial stiffness of DEM specimen is higher.
A possible interpretation is that, for the same density, the DEM specimens have a higher coordination number than the experimental samples. This interpretation is motivated by the fact that the shells are of rather irregular shape, which is not the case of our digital clusters which then tend to organize themselves spatially and thus have a greater number of contacts.

Anisotropy of cluster orientations
In this section, we focus on the influence of another very influential parameter which is the cluster orientations in the initial packings. Indeed, as for clay shells, the maximum cluster strength depends on the loading direction; the maximum strength in mode-I being much weaker than in mode-II, as seen in the Section 3. Therefore, five different specimens were prepared using a specific numerical procedure (not explained in this paper) to control the orientation of the clusters. For specimens Oedo.H and Oedo.V, the clusters were all oriented horizontally and vertically, respectively, with a given tolerance. An illustration of these two specimens can be seen in the Figure 15. For specimens Oedo.h and Oedo.v, the clusters are mainly (but not only) oriented horizontally and vertically, respectively. Finally, for speci- men Oedo.r, the clusters are randomly oriented, like it is on the real samples, see Figure 9. Figure 9 shows the mechanical responses to uniaxial loading for these five specific specimens. As anticipated, vertically oriented particles tend to improve compressive strength; it is perhaps more appropriate to say that the greater the number of horizontal particles, the weaker the compressive strength. The specimen Oedo.r (red curve) shows an intermediate strength.

Comparison of the grading curves throughout the shell-breaking process
In the previous sections, we presented the influence of various numerical and geometrical parameters on the numerical macroscopic behaviour of the clusters specimens. A broad experimental campaign has been performed by collaborators of Navier Laboratory in Paris [31], to study the mechanical behaviour from samples made of around 2 000 clayshells of identical size. In each uniaxial compression test, the evolution of stress and strain is characterized by particle rearrangements but also, and more importantly, by particle breaking -both in experiments and in DEM simulations. Therefore, we propose to perform a DEM-experiment confrontation at the particle scale to compare the evolution of grain breakage of initially identical cluster size throughout the uniaxial compression test. In the Figure 16b Once broken, cluster pieces hereafter called fragments can be made of 1 to 11 segments, the intact clusters being made of 12 segments. Assuming that all fragments pass the sieve perfectly vertically, each fragment is associated to the dimension of its annular cross section. Therefore, six types of fragment can be distinguished, and then sizes: 1 to 5 sectors and sectors with 6-11 sectors that have a sieve size identical to that of the intact cluster. Unlike clay-shell fragments, the sizes of cluster fragments are perfectly known. We can therefore assume that a sieve with a mesh size exactly equal to that of the fragment would allow the fragment to pass through, while a sieve with a slightly smaller mesh size would retain the fragment. Thus, for a sieve size corresponding to a given fragment size, we chose to consider that these fragments could pass through and be retained by that sieve. Therefore, the particle size curve obtained from the analysis of the fragments produced during the DEM simulation is not a curve but an area (hatched part in the This result is crucial in a sense that it gives confidence in the relevance of the experimentally identified values of f I and f II used for the DEM model.

Setting the model parameters
In the previous sections, all the numerical, mechanical and geometrical parameters (internal structure parameters) that can influence the macroscopic behaviour have been presented and tested. It was thus shown that, although the brittle grain model proposed here has a great amount of parameters, few of them seem to have a significative influence on the kinetics of grain breakage during sample deformation subjected to oedometric loading; in particular, the number of segments used to define a cluster is not a sensitive parameter for non-extreme compression, provided that a minimal number is used. The only really sensitive parameters are, at the grain scale, the mode-I fracture-triggering parameter f I (recall that f II does depend on f I ), and the Coulomb friction angle between the grains -all other numerical parameters showed a relatively small influence. The other relevant parameters are, at the sample scale, the initial density of the sample (number of particles per cubic meter) and the initial orientation of the shell-particles.
Therefore, using a cylindrical specimen of about 2 000 clusters, we propose to make a comparison with a real sample of about 2 000 clay-shells, both subjected to oedometric compression. For clusters as for the clay-shells, the grains are randomly oriented during the packing preparation. the set of numerical parameter used for the DEM oedometric compression is summarized in Table 2.  Table 2 Summary of the parameters used for the oedometric compression (see Figure 17).  Table 2.
The experiment performed at Euro-Géomat-Consulting EGC [36]. The axial strain ε a is computed using Eq. (8). Figure 17 compares the outcome of the DEM modelling with the experimental σ a : ε a curve. The experimental oedometer test performed by Euro-Géomat-Consulting EGC [36] includes series of the unloading-reloading (UR) cycles -UR cycles were also performed on the DEM model. With the experiments, it is noticeable that the clay-shells assembly subjected to an oedometer-type loading path can reach a compression level close to a vertical strain ε a = 1. At this stage of deformation, the clay-shells are totally crushed and reduced to powder. If the compression test had been continued, the stress would have increased sharply, as the sample would hardly increase its vertical strain. With a DEM simulation, the clusters cannot be reduced to powder, the smallest element of the model being the sector. Consequently, the sample stiffness will tend to the sector-to-sector stiffness for smaller axial strain. Thus, while the two experimental and numerical σ a £ε a curves are initially very similar, they diverge for ε a 0.45. At ε a 0.45, the numerical curve presents an inflexion point. An exponential law [37,38] seems well suited to describe the trend of the classic compression curve, i.e., for ε a ≤ 0.45. In the case of tests with a large strain range, those laws loose their validity after the inflexion point. On the experimental curve the inflexion point is less clearly marked, but it appears at ε a 0.6 such that oedometric modulus rises (σ a rises up faster). Our discrete element model limits ε a to about 0.7. This limitation results from the large size of the sector chosen (sectors are as long as the tube shape), because the tube shape was partitioned only circumferentially: N circ = 12 and N axial = 1. Increasing N axial , as shown in the Figure 6, would not only have widened the range of the maximum ε a but also strain-postponed the inflexion point. On the contrary, we checked that increasing N circ brings no benefit as long as N axial remains unchanged. By reducing the number of sectors with N circ = 12 and N axial = 1, we intended to limit the computational time (a DEM simulation of a compression test on 2 000 clusters from ε a = 0 to ε a = 0.7 runs during 6 to 7 weeks on 6 cores).
The good agreement observed between the DEM and experiments is mainly due to the choice of the parameters f I and f II , µ and the initial density of the sample. The contact elasticity parameters k I and k II play only a very limited role in the observed behavior, as long as the specimen is monotonically loaded. To justify the correct choice of elasticity parameters in the model, we have, as in the experiment, carried out UR cycles, which can be considered as quasi-elastic and therefore very dependent on the elasticity parameters of the contact between the particles. However, as it can be seen in Figure 17, the UR moduli obtained via the DEM are equivalent to those observed experimentally -the contact elasticity parameters have therefore been correctly chosen. This validates these parameters but also the procedure used to determine them (see Section 3). 6 An analytical model for compression-strain arising from shell comminution

Grain breakage analysis
In Section 3, we introduced a distinction between two breakage regimes: the cluster breakage is said to be primary when an intact cluster breaks and consequently make its internal void accessible. The breakage is called secondary when parts of already broken cluster can continue to break without "releasing" too many voids. In primary breakage, when voids become physically accessible, they play a key role in the compressibility of the whole granular specimen. Note that we refer to compressibility in terms of void reduction, exclusively (cluster sectors being unbreakable and assumed perfectly rigid).
To express the level of primary breakage, a simple ratio can be used: where N / 0 is the initial number of clusters in the specimen. The number of primarily broken clusters N broken increases from the beginning of the oedometric compression and cannot exceed N / 0 . Hence, b ranges by definition between 0 and 1; it is a proportion of broken clusters. In Figure 18, the evolution of b (red curve) compared to the stress-strain response (black curve) clearly demonstrates that cluster breakage influences the mechanical behaviour.
The evolution of b with respect to ε a shows 3 phases: in the I phase, the sample compresses mainly due to cluster rearrangements with little cluster breakage (b < 0.05). At some point, the clusters can no longer rearrange, higher inter-cluster forces develop, and the particles are more prone to break up. This corresponds to the beginning of the second phase, II . Then, the breakage occurs rapidly with an approximately constant slope of 2.093 with the increase of ε a . At the same time, the macroscopic stress increases slightly and linearly as a function of the axial strain. The phase II begins when the mechanical behavior curve becomes nonlinear, even exponential-like, for b ≥ 80%. In this phase, the higher the breakage rate b, the faster σ a increases. When all the clusters are broken (b = 1), the internal voids of each cluster are released.

An amended concept of the void ratio
The void ratio e is commonly used to describe the volume changes of compressible soils. It is a dimensionless parameter measuring the voids as a fraction of the solid phase: where the solid volume is a sum of cluster volume V i s (being itself the sum of the volumes of sectors composing the cluster) over all the clusters N / 0 :  Figure 19a shows the classical division of total volume V tot into two sub-volumes: the voids V v (grey) and the solid phase V s (red). Within an assembly of intact clusters, the hollow space inside the tubes is inaccessible by other clusters. This geometric exclusion implies that the space V i v trapped inside the intact clusters i can be assigned to the solid phase as presented in Figure 19b. When the cluster breaks, its internal void area becomes accessible and is considered as a free space (Figure 19c). Hence, we suggest to amend the definition of e into a suitable void ratio e that accounts for the evolution of accessible voids into inaccessible voids: where V accessible and V inaccessible are the whole volumes that can or cannot be filled by matter, respectively. To calculate e the breakage ratio b is introduced in the calculation of the internal inaccessible voids treated as the solid: The compression curves with respect to axial strain. Evolution of the standard void ratio e (Eq. (11)) and of the modified void ratio e (Eq. (13)) in the course of uniaxial compression at constant velocity rate. This simulation is the same as the one shown in Figure 17, with the parameters of Table 2. Figure 20 shows the evolution of both e and e as a function of the axial strain ε a , computed with Eqs. (11) and (13). As already mentioned, ε a was defined in the framework of the natural strain, that is calculated using a logarithmic function (Eq. (8)). Therefore, e always decreases non-linearly as a function of ε a . On the contrary, e evolves in a nonmonotonous manner. In the phase I , e remains around its initial value e 0 because of the clusters rearrangements and the onset of breakage. A rather stable increase of e is observed in the phase II of extensive breakage. This is a consequence of many voids becoming prone to the progressing overall irreversible compression (decrease of V ). The curve exposes the breakage inhibition with a maximum value of e (for ε a 0.42 in Figure 20). It corresponds to a maximum of e after which the sample becomes denser. This maximum indicates the axial strain at which the highly compressible regime ends. Once all the clusters are crushed, (b = 1, or equivalently V = 0 m 3 ) the two curves must merge because Eqs. (13) and (11) become identical.
The control of the maximum of e enables to tune the compressible capacities of the cluster assembly. This opens the door to the optimisation of the material strength and the cluster geometry. To this end, a prediction of the e £ε a relationship is proposed. By dividing all terms of Eq. (13) by V s , the modified void ratio can be expressed as a function of the standard void ratio e: Using Eqs. (12) and (14), the ratio V /V s becomes where the constant E 0 = V i v /V i s is the cluster scale void ratio. As a consequence, combining Eqs. (15) and (16) lead to where e depends on the void ratio e, the initial void ratio E 0 of the intact clusters, and the evolution of the rate of breakage b during the compression.
To express e as a function of ε a , both b and e need to be a related to ε a . Firstly, b £ ε a relationship was assumed to be linear with slope of 2 (which corresponds to the slope 2.093 shown in Figure 18 with the logarithmic strain definition in uniaxial compression defined as where subscript / 0 denotes the initial state. We can then write that exp(ε a ) = e 0 + 1 e + 1 (20) and finally, with the initial void ratio e 0 serving as the input.
In summary, the model distinguishes two regimes depending on the breakage level: Eqs. (17), (18) and (21) constitute the prediction of e for b < 1. Once all the clusters are broken, b = 1, Eq. (18) is no longer valid and the classical definition of the void ratio is used. Hence, the model can be written as: 6.3 Testing the model Figure 21 shows the outcome of the prediction model for different thicknesses of shells (t = 2.4; 3.6; 4.8 mm, which corresponds to E 0 1.062; 0.498; 0.234, respectively). These predictions (solid lines) are also compared with the equivalent simulations (dashed lines). The model was built upon a very basic concept and a simple assumption on breakage kinematics, Eq. (18). Although the model is unable to replay all the details in the evolution of the e , it captures quite well the main features, especially when the void ratio switches from e to e. In essence, this model has structured a solid foundation for a more complex prediction of the constitutive e £σ a relationship to be derived in the future.

Summary
This paper gives some insight into the micro-mechanical origins of the high compressibility of the assembly that comprises extremely porous grains. Hereinbefore, we broke down the complexity of the model into several steps and build it up from micro to macro scale, as follows.
In the Section 2, by bonding non-spherical particles (called sectors), we have designed a discrete element (DE) model of a tubular grain (i.e., a cylindrical shell) capable of simulating brittle crushing. An improvement of the rigid clustering technique that allows efficient modeling of particles of any complex shape has been introduced. To our knowledge, clustering sphero-polyhedral elements into breakable grain is a rather innovative idea in the field of crushable geo-materials. The necessary numerical implementation were included in the in-house code Rockable.
The tensile failure has been achieved at the shell-scale using uniaxial radial compression (URC) in the Section 3. The variability of the shell strength was described through a Weibullian distribution. The experimental campaign contributed to numerical study not only with the determination of the tensile strength and its variability, but also with a characterisation of the breakage manner, helpful to construct the geometry of the cluster. The model successfully captured the shell mechanical behaviour in URC under the assumption of the linear elasticity. Still, the first consequences of the model simplifications have been identified due to the lack of (i) non-linearity of the contacts laws (delaying the arise of intra-cluster forces) and (ii) separation between the crack initiation and propagation. A parametric study has been undertaken to clarify the influence the elastic and yielding parameters. The number of sectors per cluster have been reported to affect the values of discrete parameters. Finally, a robust set of parameters and a cluster structure tradeoff were selected. Section 4 comments on the "structure" of actual specimens. Acknowledging the influence of initial state on the mechanical behaviour, we characterised the bulk density and the intrinsic anisotropy of shell orientations. While the former was a simple and straightforward measure, the latter required the use of more advanced techniques. First, we have proposed a geometry-based approach to efficiently detect the tubes in 3D binary image resulting from X-ray scanning. It was then followed by the statistical analysis of the orientations. Despite a natural tendency of the shell to embed horizontally, the remaining orientations are uniformly distributed within the sample. As a consequence, under oedometeric compression, the tensile failure of shell will be dominant.
In the Section 5, a systematic investigation of the responsiveness of the macroscopic model to changes in the DEM parameters has been conducted. Beyond that, the sensitivity analysis has been extended to the influence of the initial state and of the cluster geometry. To moderate the mechanical response numerically, two parameters have been found crucial. First, it has been shown that the axial stress, that is the sample bearing capacity, depends on the tensile strength of constituents which is controlled by the contact parameter called the tensile yield force (see Section 2). Second, the initial density shapes (either smoothen or sharpen) the trend of the constitutive relationship. It has been clearly demonstrated that the significant size of the sector in longitude size of longitude sector restricted the range of the strain. On the other hand, the evolution of the grading size distribu-tion highlighted that the micro-mechanics of packing has not been harmed with this simplification. Decrease of the sector size extends the functionality of the model by postponing the end-deformation of the primary breaks, but becomes computationally expensive. The effort undertaken has made it possible to achieve an important objective: to reproduce the experimental result provided by EGC and Laboratoire Navier, permitting to satisfy the technological specification from fabricating the shell, to the prefabrication of segment lining.
The Section 6 evidenced that the high compressibility originated from the extensive internal pore collapse triggered by the primary breakage. The standard void ratio has been redefined in the framework of accessible and inaccessible space, instead of solid and void. The overall void ratio decreased progressively respecting the stain-controlled compression. In contrast, the redefinition of porosity taking into account the accessibility of voids has led to a strong change in this evolution. The accessible void ratio increases despite the fact that the macroscopic volume shrink. In other words, the packing experienced a "porosification", and once the extensive breakage has been finished, only densification took place (i.e., a limited rearrangement of rigid sectors allowing high transmission of forces through the elastic stiff contacts). The end of this porosification was associated with the point at which the high compressibility could not be exploit any longer. A general, geometrical development of a prediction model, including some analytically determined relationships, has been conducted to foresee how far in the macroscopic compression the internal porosification take place. In the comparison with modelling of the thicker shells, the estimations were rough but robust. This novel concept for compressibility analysis is then left as the basis for next upgrades and developments of the authors.