Breakage in quasi-static discrete element simulations of ice rubble

Breakage of particles plays a key role in force transmission in granular materials. Discrete element method (DEM) simulations are often used to model granular materials, but modeling particle breakage in them remains a challenge. Models for breakage of non-spherical particles are scarce and often the existing models are computationally heavy to be used in simulations with large numbers of particles. To address this, the present study develops a particle breakage model for quasi-static DEM simulations of non-spherical particles that fail due to shear. The model is novel, since it is based on experimental observations and high resolution modeling. Breakage models based on experimental evidence are rare as it is often virtually impossible to gain detailed data on the mechanisms related to breakage. The developed particle breakage model was integrated into a DEM code and direct shear box experiments on ice rubble, a granular material consisting of ice particles, were then simulated. Accounting for particle breakage in DEM simulations improved their accuracy: simulations were compared to experiments and the results were found to be in better agreement when particle breakage was taken into account. The effect of particle breakage on the shear strength of a granular material was found to be independent of particle size, decreasing fast with increasing particle strength. The combined effect of shear box length and breakage was also studied. The results showed that the strength of a granular material may be determined reasonably well with a shear box that has a box to particle length ratio greater than 60.


Introduction
Load transmission in a granular material under quasi-static loading occurs through force chains, chain-like sequences of its particles under high compression [1].Two important mechanisms that limit the magnitude of the load transmitted by a force chain are (1) buckling of the force chain [2] and (2) breakage of the particles within the force chain [3].A common way to model the behavior of granular materials is to employ simulations based on the discrete element method (DEM) [4].While force chain buckling is inherently described by DEM simulations, modeling breakage relies on estimates on stresses within individual particles.Despite earlier work on breakage, describing it reliably remains a challenging task.Experimental data on mechanics of breakage are rare, and often, depending on the material and application, impractical to obtain.
A particle breakage model has to include (1) a failure criterion for the particles and (2) a technique to generate particle fragments.For spherical particles, a common approach is to use contact forces to compute average principal stress for each particle and compare that with a relevant bulk strength measure, often compressive strength, of the material of the particles [5][6][7][8][9][10].The most straightforward approach for generating fragments in this case is to replace the breaking particle with progeny fragments whose dimensions are determined from a predefined fragment size distribution [11][12][13][14][15].This approach assumes a uniform stress for a particle and that the failure plane passes through the centroid of the particle.Thus, the approach may work well for spherical particles.The applicability of such an approach in simulations with even moderately elongated particles is questionable anyhow because such particles often have highly non-uniform stress distributions and the failure planes cannot be justifiably assumed to pass through the centroids of the particles.This is, for example, the case when the contact forces act mainly on one end of a particle only.Breakage models based on principal stress values have been, nevertheless, used in simulations with non-spherical particles [16][17][18][19][20][21].This is likely due to such models being fairly straightforward to implement.
To mitigate these issues, boundary elements have been used in DEM simulations to estimate stress distributions within the particles [22][23][24][25].In addition, techniques combining finite-and discrete element methods (FEM-DEM) have been used in simulating breakage [26][27][28][29][30].In FEM-DEM, particles are meshed into finite elements, which are then used to model the failure.Both techniques allow the stress state and the onset of breakage of a given particle to be more accurately represented and yield better estimates for crack paths to generate realistic fragments.Such approaches, however, result in high computational burden, rendering https://doi.org/10.1016/j.ijmecsci.2023.108595Received 20 March 2023; Received in revised form 30 June 2023; Accepted 2 July 2023 them unfeasible for simulations with large numbers of particles.In highly dynamic simulations, one the other hand, the cumulative impact energy of the particles is often compared to the fracture energy of the material to determine if breakage occurs [11,[31][32][33][34], but such an approach is not suitable for applications with quasi-static deformation.Hence, there exists a research gap for a simplified particle breakage model that can be used in large-scale quasi-static DEM simulations of non-spherical particles.
The present paper introduces a new, rather simple, model for breakage of particles in DEM simulations.The model is aimed at quasi-static simulations of granular materials consisting of particles of moderate aspect ratios.The particles are assumed to break due to shear failure governed by the strength of the particle material rather than instantaneous crack propagation caused by stress concentrations.The model is applied to simulations of ice rubble, a granular material consisting of a rubble pile of ice blocks, for which the engineering problems often exhibit the aforementioned features.Importantly, our particle breakage model is based on direct experimental observations and high-resolution bonded particle model simulations on breakage of particles in ice-to-ice contacts [35,36].We compare the results from simulations with and without the breakage model to experimental data of direct shear box experiments on ice rubble [37,38] and then study the effect of breakage on the shear strength of ice rubble using the results of direct shear box experiments.
Effect of breakage on mechanical properties of granular materials has been studied extensively for geo-materials, such as sands and rockfills [39,40].Similar studies on ice rubble are scarce.Hopkins and Hibler [41] used a DEM tool with a bending-failure-based particle breakage model to simulate the direct shear box experiments.The work showed that the dilatancy of the rubble in direct shear box experiments decreases when the particles are allowed to break.Kulyakhtin [42] used a continuum ice rubble model with a breakage parameter to simulate ice loads on a structure moving through ice rubble and concluded that the ice loads were governed by the accumulation of the rubble, but not breakage.However, continuum models are not capable of capturing force chains within ice rubble, while earlier work suggests that they may have a crucial role on ice rubble behavior [43,44].Therefore, there exists a research gap on the understanding of the effect of particle breakage on the mechanical properties of ice rubble, which motivated the present study.
In what follows, we first introduce our particle breakage model and describe the mechanics of our DEM simulations in Section 2. Section 3 demonstrates the applicability of the model by comparing our simulations to experiments.After this, we use the model to study the effect of breakage on the mechanical behavior of shearing ice rubble.We end the present paper with a discussion on our results in Section 4 before concluding it in Section 5.

Methods
This section introduces the breakage model.The model is integrated into a three-dimensional DEM code, described in detail by Polojärvi [45].The code is parallelized and fairly standard on a general level: simulations are explicit, particles interact through pairwise contact forces, a central difference scheme is utilized for time-stepping, and polyhedrons are used to describe arbitrarily shaped particles.We first describe the calculation of contact forces, then the breakage model, and finally the modeled direct shear box experiments.

Contact forces
The DEM simulation tool used here utilizes a soft contact model, that is, the contact force for a pair of interacting particles is determined based on a small overlap between them.The point of application of the contact force is at the centroid of the overlap volume.In brief, the contact force between a pair of particles,  =   +   , consists of a normal and a tangential component,   and   , respectively.An elasticviscous-plastic contact force model is used to calculate   [46].The elastic and viscous components of   are solved by using the gradient of overlap volume and its rate of change, respectively [47,48].The plastic component of   , describing local yielding or crushing of the material of particles in contact, is based on area of contact.Tangential compliance and friction between the particles result in   [46].Importantly, the contact model can be parameterized by using the material properties of the particles as described by Polojärvi [45].

Breakage model
The breakage model is based on checking if a predefined failure criterion is met on potential failure planes within each particle in the simulation.This is performed by utilizing so-called base planes, defined for each particle during pre-processing (Fig. 1a).The base planes are used to reduce the originally three-dimensional problem into two dimensions.Actual failure planes, when formed, are perpendicular to the base planes.The base planes are defined in a local coordinate system of each particle, which moves with the particle.The model is foremost developed for particles with polyhedral shapes, but works for particles of arbitrary shapes, including spherical ones.Fig. 1a shows an example of a prismatic polyhedral particle for which the base plane is chosen to align with the face with the largest surface area.This would be a natural choice for an ice floe floating in water, often under inplane loading by other floating floes [21].Another example would be to assume a particle having anisotropic strength and known to be likely to fail due to shear on the plane having the lowest shear strength.In this case, a natural choice for a base plane would be a plane perpendicular to the plane of lowest shear strength, since then the failure plane would align with it.For a cuboid-shaped particle, on the other hand, the first choice could be three orthogonal base planes.Further, depending on the desired accuracy, a near-spherical particle of isotropic strength may need a number of base planes.Below we study an application, where choice of base planes is straightforward, but applications where the choice of base planes is more uncertain are left for future work.
At each time step of a simulation, contact forces applied on each particle and its geometry are projected onto its base plane, the latter forming a two-dimensional polygon (Fig. 1b).The quasi-static force equilibrium of this polygon is then analyzed to solve nominal force and stress components on potential failure planes (Fig. 1b).The contact forces are projected to the base plane simply by where  *  is the projection of contact force   to the base plane and   is the unit vector perpendicular to the base plane.The force component perpendicular to the base plane causes a bending moment on the particle, but it is neglected as the bending failure of particles is not the dominant failure mode under compressive contact forces.The points of application of all contact forces are similarly transformed onto the base plane.
The breakage model assumes that each failure plane goes through a point of application of one of the contact forces.Fig. 1b illustrates how this feature was implemented.For each contact force, the polygon boundary is discretized into a number of mesh points for testing planes of different orientations for potential failure.The figure shows a mesh created for checking if  * 1 leads to breakage.The figure also shows a line belonging to a potential failure plane, spanned between the point of application of  * 1 and mesh point  (as mentioned above, the actual failure plane is always assumed to be perpendicular to the base plane).
Each of the planes is then tested for failure by analyzing the quasistatic equilibrium of the two-dimensional polygon on the base plane.Nominal shear and normal force components,  *  and  *  , respectively, acting on a potential failure plane can be calculated as where  is the number of contact forces acting on the fragment and  *  and  *  , respectively, are the tangential and normal unit vectors of a potential failure plane.Nominal shear and normal stresses are then, respectively, where  * is the area of the potential failure plane, calculated by multiplying the length of the line spanned for describing the orientation of a given plane on the two-dimensional polygon by the average thickness of the particle.The sign of compressive normal stress is here assumed to be positive.Nominal stresses  and  are then used to check for particle failure by using them in conjunction with a failure criterion,  , suitable for the material studied.For many materials, the Mohr-Coulomb failure criterion is used.The material parameters in this criterion are internal cohesion and friction,  and , respectively.For admissible nominal stresses and breakage occurs if  ≥ 0. When the failure criterion is met, a failure plane is defined along the plane with  ≥ 0, dividing the parent particle into two fragments.Thus, both the shape and the size of the two fragments are determined by the location of the failure plane within the parent particle.This ensures realistic fragment shapes as well as mass conservation in the simulations.The breakage model is then applied to new fragments as the simulation progresses.In the simulations below, we allowed the particles to fragment two more times each after the initial failure.This limit was set to avoid very small particle fragments.
A pseudocode for the breakage model is presented in Algorithm.

Direct shear box experiments
Below we compare our DEM simulations to pseudo two-dimensional direct shear box experiments by Pustogvar et al. [37] (Figs.2a and  b).The ice rubble specimen modeled had the dimensions of 600 mm × 400 mm × 40 mm (length × width × thickness).Two rubble types were used, one consisting of ice blocks of 30 mm × 20 mm × 40 mm and another consisting of ice blocks of 60 mm × 40 mm × 40 mm, referred to as small and large particles, respectively (the experiments were pseudo two-dimensional as the thicknesses of the particles and the shear box were equal).In the simulations, the box was filled with ice rubble by first using gravity deposition and then by applying a vibrating motion to the box-the initial porosity of the simulated rubble was Define two new particles in the DEM simulation for the fragments 19: end if same as in the experiments.Confining pressure on the rubble,  , was applied by a varying force on the shear box cover, which was allowed to rotate.Other walls were not allowed to move.In the experiments,  had the magnitudes of 5.75 and 11.03 kPa, while simulations involved values of  varying in the range 5.75 … 22.06 kPa.The experiments were quasi-static, since the rate of displacement, δ, of the upper half of the box was only 0.02 ms −1 .Table 1 summarizes the main parameters used in the simulations.
Shear force, (), as a function of displacement  was recorded (Fig. 2a).() could then be used to derive the instantaneous shear stress,   (), within the rubble specimen from where () is the area of the shear plane.For ice rubble and ice engineering applications, it is also useful to determine the mean shear Polojärvi et al. [38] simulated the first quarter of the direct shear box experiments (0 <  < 0.14 m) studied here.Pustogvar et al. [37], however, reported that particle breakage affected the results of the experiments during the later stage of the experiments; thus, the present study simulated the whole experiment to capture the particle breakage and its effect on rubble behavior.

Results and analysis
This section starts with a comparison of our DEM simulations with experimental data.Then it continues by discussing the effect of breakage on the direct shear box experiment results.Finally, the combined effect of box length and breakage on the Mohr-Coulomb material model parameters is investigated.The present study is based on 395 simulations in total with 79 different parameter configurations and 5 simulations per each configuration with different initial rubble packings.Supplementary material of the present paper includes a MATLAB implementation of the breakage model, which the reader can use for model verification.3a and b are from simulations with small and large particles, respectively.The confining pressure,  , was 5.75 kPa and in the simulations with the breakage model, the internal cohesion, , and internal friction, , had the values 250 kPa and 0.75, respectively, chosen after Prasanna et al. [36].The figure also shows the  −  curves from the experiments by Pustogvar et al. [37] for the corresponding cases.Experimental records only reached  ≈ 0.4 m due to technical difficulties with the setup [37].

Comparison to experimental data
As Fig. 3 illustrates, the general features of the  −  records from the simulations and experiments are similar.Both show several distinct load peaks, here referred to as peak load events.The breakage model appears to improve the simulation results, as the simulations without it yield several load peaks of considerably higher magnitude than those in the experiments.It is worth emphasizing here that the simulated and experimental  − records do not show one-to-one agreement, since the initial rubble configurations were not the same.Nevertheless, records from the simulations with the breakage model are more representative of experimental curves with regard to the peak load magnitudes.The  −  records, however, from the simulations without the breakage model initially follow those from the simulations without it.The two records deviate from each other as the shear force, , builds up towards the first peak load; clearly breakage affects the mechanism that limits peak loads, as reflected further by the magnitudes of them being in general smaller in simulations with the breakage model than in those without it.Fig. 4 shows that, overall, the simulations allowing breakage yield  and   data more closely resembling that from the experiments than the simulations without it. and   in the simulations and experiments differed on average by 11% and 14%, respectively, with the breakage model.Simulations without the breakage model led to errors of 23% and 41% respectively.Thus, the values for  and   are more close  to experimental ones in the simulations with the breakage model than in the simulations without it.This indicates that the particle breakage model improves the accuracy of DEM simulations of shearing granular materials.Furthermore, the  and   values from the simulations with the breakage model are lower than those from the simulations without it.It is also interesting to notice that the values of   from the simulations show smaller scatter when the particles are allowed to break.

Breakage and peak load events
Breakage clearly limits the peak shear load values.In the simulations without the breakage model, load peaks are associated with the formation and subsequent buckling of force chains [38].In the simulations with the breakage model, particles can fail before the load transmitted by a force chain reaches that needed for its buckling.This difference is illustrated by Fig. 5 presenting two sequences of images from the first peak load event in the simulations with the shear forcedisplacement,  − , record presented in Fig. 3a.In Fig. 5, the left and right columns, respectively, show the ice rubble in the simulation without and with the breakage model.The figures also illustrate the first principal stress of the so-called particle stress tensor [49,50].This is used to illustrate the main compressive load and its direction within the particles.For the figures, the principal stress is normalized by its maximum value during the peak load event in the simulation without the breakage model, and it is only shown for particles for which its normalized value exceeds the threshold of 0.1.At the start of the peak load event, a force chain, that is, a sequence of particles under high compression, starts to form in both simulations as indicated by the series of highly stressed particles (Fig. 5a).As the load increases, some of the particles within the force chain break in the simulation with the breakage model (Fig. 5b).After breakage, the force chain disappears, while in the other simulation the particles in the force chain are compressed further (Fig. 5c).Particle breakage is, as expected, one of the force chain failure modes (the force chain in the simulation without the breakage model eventually buckled).As described above, the  −  record of Fig. 3a reflects this behavior: the magnitude of the first peak load in the simulation with the breakage model is significantly lower than that in the simulation without it.After the first peak load event, the  −  records deviate between the two simulations, as the arrangement of the particles differs between them, mostly due to new fragments in the simulation allowing breakage.Also, in the simulations with the breakage model, some force chains collapsed due to buckling, and thus, breakage was not the only failure mode in them.

Particle strength and specimen strength
As seen above, particle strength affects the shear force, , measured in a direct shear box experiment which in turn affects the measured shear strength,   , of the rubble specimen.Fig. 6 illustrates how this effect changes as a function of internal cohesion, , of the particles.The figure shows the magnitudes of mean shear resistance and maximum shear strength of the rubble,   and    , respectively, plotted against .
Results are shown for small and large particles and for both confining pressures,  = 5.75 and 11.03 kPa, used in the experiments by Pustogvar et al. [37].Again, five simulations were run for each case, and thus, the figure shows the mean values and the standard deviations for the data.The values presented were calculated for the shear box displacement interval  = 0 … 0.3 m to mitigate the effect of shear box walls [38].The figure further presents fits for the mean values defined by using the MATLAB curve fitting toolbox with an asymptotic function,  =  ′ − ( ′ −  ′ ) exp −∕ ′ , where  ′ ,  ′ and  ′ are the fitting parameters.
In each case, the fit was required to go through the mean value of the simulations without the breakage model, indicated by the data at  = ∞ kPa in the figure.Data in Fig. 6 show that both   and    increase with particle strength and size, and with confinement.This behavior is in line with earlier studies on granular materials, where shear strength increases with increasing particle size [51,52] and confinement [53].Moreover,   and    increase towards the values of   and    of simulations without the breakage model as  increases; thus, the effect of particle breakage diminishes with increasing .A similar effect of  is seen for simulations with both values of  .The mean relative difference between the results from the simulations for the two levels of  was 33% in   and 31% in    .For   , a measure for shear resistance, the effect of breakage virtually vanishes at around  = 400 kPa.With  > 400 kPa, the values of   are about 90% of those from the simulations without the breakage model.For    , a measure of instantaneous shear strength, breakage affects the results still at  = 500 kPa and the effect of  on     is more pronounced with a higher  .Fig. 7 shows normalized values of mean shear resistance and maximum shear strength of the rubble,   and    against internal cohesion, , of particles.These data were normalized with the results from the corresponding simulations without the breakage model.The figure combines the results from simulations with confining pressure  = 5.75, 11.03, 16.45, and 22.06 kPa, as importantly, the normalized data in all simulations behaved very similarly.The figure also presents fits of similar form to those in Fig. 6 on the normalized data.Fig. 7 shows that the data for normalized   and    for small and large particles overlap.This suggests that the effect of breakage on shear resistance,   , and instantaneous maximum shear strength,    , of the rubble is independent of particle size.Normalized   and    converge toward unity as  increases, and the effect of  on the data is predicted to vanish at about  = 1.0 MPa.
Further, we also investigated the effect of internal friction, , [Eq.( 4)] on the measured shear strength of ice rubble.Simulation results did not show significant difference with regard of shear force, , or the block failure events for the tested  values of 0.25. . .0.75.In addition, the effect of material properties of particles (Young's modulus and friction coefficient from Table 1) on  were studied briefly by varying them within one order of magnitude.The values of mean and maximum of shear force, S and   , showed an increasing trend with increasing Young's modulus and friction coefficient.However, they remained in general lower when the breakage was allowed compared to those in the simulations without the breakage.

Breakage and experimental set-up
Direct shear box experiment results are known to be affected by the box length to particle length ratio, ∕, of the setup used [38].Is this effect similar for simulations with and without the breakage model?We studied this by simulating direct shear box experiments with boxes having lengths in the range of  = 0.6 …3.0 m (Fig. 2a), that is, box length to particle length ratios of ∕ = 20 … 100, where  is the largest side length of small particles; only small particles were used.Simulations with and without the breakage model were performed, and for the simulations with the breakage model, internal cohesion  = 250 kPa [36] was used.Four different confining pressure  values of 5.75, 11.03, 16.45, and 22.06 kPa were tested to see if specimen shear strength,   , values follow the Mohr-Coulomb material model for the given shear box geometries [Eq.( 6)].For each parameterization, five simulations with different initial particle arrangements were run.
Fig. 8 presents the data for the mean shear resistance and maximum shear strength of the rubble   and    (Section 2.3), respectively, against  for shear box lengths  = 0.6 and 3.0 m.Data for simulations with and without the breakage model are shown for both lengths.Both   and    show a linear increase with  as illustrated by the trend lines fitted to the data.This indicates that the Mohr-Coulomb material model applies for the results yielded by a box with a given .Introducing breakage affects   and    as internal cohesion of ice rubble,   (intercept of the trend line), and friction angle,   (gradient of the trend line), are lower for simulations with the breakage model than for those without it.The first of these results would be expected based on the above considerations, as the shear force, , was lower in simulations with the breakage model than in those without it (Fig. 4).
Even more importantly, Fig. 8 shows that an increase in  from 0.6 to 3.0 m leads to a decrease in the values of   and   in all simulations.This should not be the case if the measured   and   were material properties, which motivated us to look for ∕ ratios at which they become independent of .Figs.9a and b show the values of   and   , respectively, against ∕ derived by using the results from simulations with different  values.The data in the figure were calculated by using both   and    from simulations with and without the breakage model.Additionally, Fig. 9 shows a fit of form  =  ′ − ( ′ −  ′ ) exp − ′  for   and a linear fit for   .
As seen from Fig. 9a, internal cohesion of ice rubble,   , reaches an asymptotic value at box length to particle length ratio, ∕ ≈ 60, that is, at  ≈ 1.8 m.For   derived using the instantaneous maximum shear strength,    , the asymptotic value is about 7 kPa.This is double the asymptotic value of   derived using   , which describes the shear resistance.Breakage improves the results with small ∕, or in other words, when experiments are performed for materials with particles that break, a shorter box appears to suffice for estimating   .Interestingly, when ∕ > 60 breakage has no effect on the values of   .
The values of ice rubble friction angle,   , however, did not appear to converge as Fig. 9b shows, but instead decreased with increasing   ∕.One likely reason for   not converging is non-uniform distribution of confining pressure,  , within the rubble: Polojärvi et al. [38] showed that  is often transmitted through few nearly vertical force chains and, thus, an increase in  does not necessarily results in increased frictional resistance on the whole shear plane but only on few negligibly short parts of it.The simulations with the breakage model yield smaller   values than the ones without it, but the linear fits on the data suggest the rate of decrease is about the same for both types of simulations.The gradient of the trend line derived using mean shear resistance,   , is −0.16 and −0.14 in the simulations with the breakage model and without it, respectively, while the same using the maximum shear strength,    , is −0.23 and −0.16, respectively.

Discussion
The approach to modeling breakage taken here, based on nominal force and stress components, and the Mohr-Coulomb failure criterion, was demonstrated to predict breakage reliably by our experiments [35,36].In these experiments, the material tested was ice, but we believe this rather simple approach can be used for modeling the breakage of particles of other materials, given the breakage is due to shear failure governed by the strength of the particle material, rather than instantaneous crack propagation due to stress concentrations.An example of such material is rocks [54].In addition, the model can be modified to work with other materials by choosing an appropriate failure criterion.It is also worth mentioning here that there have also been analytical models on ice edge failure based on shear and the Mohr-Coulomb failure criterion [55,56], but not for ice features other than an ice sheet.
DEM simulations of ice usually only account for local contact failure of particles by setting a plastic limit for the contact force.This limit is practically always based on the compressive strength of ice [43,44,46,57,58].Our work in Prasanna et al. [35,36] shows that this approach is often inaccurate for ice-to-ice contacts, which in turn motivated us to implement the breakage model introduced here.The results from the simulations with and without the breakage model confirm that the previously used approach is not sufficient, but rather breakage has to be accounted for when modeling ice rubble.This is because load transmitted by force chains is often limited by breakage of particles within it (Section 3.2).Moreover, fragments generated from breakage fill the voids within the rubble during particle rearrangement, which in turn reduces the porosity of the rubble.Hopkins and Hibler [41] were the first ones to include any type of a breakage model for particles within ice rubble.They used a model assuming bending failure.Our work in Prasanna et al. [35,36] does not support this assumption, and indicates instead that the model should be based on shear failure of ice particles when the aspect ratio of particles is low.For particles with high aspect ratios, however, a breakage model based on bending failure may work well.
Our results on breakage and peak load events indicate that often the force chain collapse due to particle breakage occurs at lower load levels compared to force chain buckling.Thus, the load-carrying capacity of ice rubble is lower when particles are allowed to break, which in turn decreases the shear resistance,   , and the maximum instantaneous shear strength,    , of ice rubble.The normalized   and     results of Fig. 7 showed about 20% reduction of shear strength of ice rubble for internal cohesion of ice,  = 250 kPa, chosen following the experiments of Prasanna et al. [35].This means that previous ice load estimates obtained from DEM simulations of ice-structure interaction processes without a particle breakage model [43,44,57] may have had some disparities, and that their results have to be used cautiously for interpreting limit loads on structures.Moreover, the normalized   and    values (Fig. 7) further showed that the effect of breakage on the shear strength of ice rubble is independent of particle size.This implies that using smaller particles in DEM simulations in order to mitigate the effects of force chains is not effective, but that a particle breakage model has to be used.This is due to the distinct role of particle breakage as a load limiting mechanism.It is also worth emphasizing here that normalized   and    of Fig. 7 results are not contradictory to the established understanding of the breakage of granular materials, where particle breakage increases with the increasing particle size [59][60][61].In an experiment, larger particles are prone to contain more flaws, which in turn increases the probability of their failure and decreases their strength.Normalized   and     results indicate that there is no apparent effect of particle size on breakage since normalized   and    of small and large particles overlap; however, as  decreases for a given particle size, normalized   and    also decrease indicating increasing particle breakage.Thus, an ice rubble specimen consisting of large particles with lower particle strength would exhibit more breakage than a rubble specimen consisting of smaller particles but higher particle strength.This is in line with the hypothesis that the effect of particle size on breakage is truly related to the strength of the individual particles.Moreover, this implies that the developed breakage model can be used to model particle breakage in different scales simply by defining internal cohesion, , as a function of the particle size.
Our simulations further revealed that internal cohesion of ice rubble,   , and friction angle,   , derived from the results of a direct shear box experiment depend on the length, , of the shear box used.This has been observed to occur in previous studies as well [62,63].Breakage affects   when the box length to particle length ratio ∕ < 60, but with ∕ > 60 the difference between the values for   from simulations with and without the breakage model was less than 10%.The tendency of breakage to affect the results with small boxes is partly explained by the ∕ ratio changing during an experiment due to breakage: when  is relatively small, the ∕ ratio changes significantly due to new ice fragments yielded by breakage.Further, with small ∕ ratios, breakage is more likely to limit the load transmitted by force chains due to their limited length.Surprisingly,   showed a continuous decrease with increasing  instead of converging.It is, however, clear that   > 0 for any material and, thus, it is evident that the experimental set-up does not produce physically sound results for   of the rubble.As was described in Section 3.4 and based on Polojärvi et al. [38], we believe this is due to non-uniform distribution of confining pressure,  , within the simulated rubble specimens.A potential remedy for this would be to consider an experimental set-up allowing more even distribution of  on the ice rubble [64].

Conclusions
The present paper developed a simplified model of particle breakage for inclusion in DEM simulations of granular materials.The breakage model is based on experimental observations [35] and high resolution numerical analysis [36] of ice blocks of meter scale breaking under quasi-static compressive ice-to-ice contact forces.It was here applied to simulations of direct shear box experiments on ice rubble [37].While the application here was on ice rubble, we believe our approach can be applied to other materials as well.Based on the foregoing results, the following conclusions can be drawn: • The particle breakage model developed in the present paper is able to successfully capture the effect of particle breakage in DEM simulations.• Breakage limits the load transmitted by force chains, as particles within a force chain may break before the chain buckles (Fig. 5).Modeling this in DEM improves the accuracy of the simulations and leads to 20% lower peak load values in direct shear box simulations (Section 3.1).• The effect of breakage on the shear strength of a granular assembly is independent of the particle size for a given particle strength, and decreases with an increase in the internal cohesion of the particles (Fig. 7).For ice rubble, the effect vanishes at about 1.0 MPa.• The internal cohesion of ice rubble can be estimated reliably with a shear box having a box length to particle length ratio ∕ > 60 (Fig. 9).Friction angle values did not converge for the tested ∕ ratios.
Including particle breakage in DEM simulations improves agreement between simulations and experimental observations.The first steps of future work should explore the effect of base plane choice on simulation results.

CRediT authorship contribution statement
Malith Prasanna: Conceptualization, Developed the particle breakage model and integrated it into the existing DEM tool, Prepared and performed the simulations, Processed the data from them, Analyzed the results, Writing -original draft.Arttu Polojärvi: Conceptualization, Discussed the results, Reviewed and edited the manuscript.

Fig. 1 .
Fig. 1.(a) Illustration of the base plane and the contact forces resolved onto the base plane.In the figure,  are the contact forces and  * are the contact forces resolved onto the base plane.(b) Discretization of the particle boundary projection on the base plane.Here,  is assumed to be the failure plane and,   and   are the nominal shear and normal forces on the failure plane respectively.

Algorithm 6 : 9 :
The particle breakage algorithm 1: Choose a base plane for the breakage model 2: Project particle boundary onto the base plane 3: Project contact forces and their points of application onto the base plane:  * contact force resolved onto the base plane [Eq.(1)] 4: : No. of contact forces acting on the particle 5: for  = 1 to  do Discretize the projected particle boundary of the element 7: : No. of mesh points on the projected particle boundary 8: for  = 1 to  do Define a plausible failure plane between the point of application of  *  and the mesh point   10: Find the nominal shear and normal forces,  *  and  *  , acting on the failure plane [Eq.(2)] 11: Find the nominal shear and normal stresses,  and , acting on the failure plane [Eq.(3)] 12: Apply the Mohr-Coulomb failure criterion and calculate the critical value for the failure criterion,  [Eq.(4)] 13: end for 14: end for 15: if max ( ) ≥ 0 then 16: Failure will occur along the plane with max ( ) 17: Split the particle along the failure plane 18:

Fig. 2 .
Fig. 2. (a) Set-up of the direct shear box experiments simulated.In the figure,  and  are the length and the height of the shear box, respectively.Further,  is the shear force,  is the confinement and  is the area of the shear plane.(b) Photograph of an experiment with small particles from Polojärvi et al. [38].

Fig. 3
Fig.3presents a comparison of typical shear force-displacement,  − , records from simulations of direct shear box experiments with and without the breakage model.Records in Figs.3a and bare from simulations with small and large particles, respectively.The confining pressure,  , was 5.75 kPa and in the simulations with the breakage model, the internal cohesion, , and internal friction, , had the values 250 kPa and 0.75, respectively, chosen after Prasanna et al.[36].The figure also shows the  −  curves from the experiments by Pustogvar et al.[37] for the corresponding cases.Experimental records only reached  ≈ 0.4 m due to technical difficulties with the setup[37].As Fig.3illustrates, the general features of the  −  records from the simulations and experiments are similar.Both show several distinct load peaks, here referred to as peak load events.The breakage model appears to improve the simulation results, as the simulations without it yield several load peaks of considerably higher magnitude than those in the experiments.It is worth emphasizing here that the simulated and experimental  − records do not show one-to-one agreement, since the initial rubble configurations were not the same.Nevertheless, records from the simulations with the breakage model are more representative of experimental curves with regard to the peak load magnitudes.The  −  records, however, from the simulations without the breakage model initially follow those from the simulations without it.The two records deviate from each other as the shear force, , builds up towards the first peak load; clearly breakage affects the mechanism that limits peak loads, as reflected further by the magnitudes of them being in general smaller in simulations with the breakage model than in those without it.Fig.4comparesthe data from the simulations and experiments further by showing the mean and maximum of the shear force,  and   , respectively.Data for confining pressure  = 5.75 and 11.03 kPa are included.Since in both the experiments and simulations the details of the  − records for a given parameterization depended on the initial arrangement of the particles, the figure shows data collected from five repeated simulations for each case.Meanwhile, the experimental data are from three repeated experiments.For each simulation and experiment,  was defined by simply taking the average of , and   as the maximum value of .The results are for the displacement interval  = 0 … 0.3 m, since the effect of shear box wall on  is more pronounced for larger .Fig.4showsthat, overall, the simulations allowing breakage yield  and   data more closely resembling that from the experiments than the simulations without it. and   in the simulations and experiments differed on average by 11% and 14%, respectively, with the breakage model.Simulations without the breakage model led to errors of 23% and 41% respectively.Thus, the values for  and   are more close

Fig. 4
Fig.3presents a comparison of typical shear force-displacement,  − , records from simulations of direct shear box experiments with and without the breakage model.Records in Figs.3a and bare from simulations with small and large particles, respectively.The confining pressure,  , was 5.75 kPa and in the simulations with the breakage model, the internal cohesion, , and internal friction, , had the values 250 kPa and 0.75, respectively, chosen after Prasanna et al.[36].The figure also shows the  −  curves from the experiments by Pustogvar et al.[37] for the corresponding cases.Experimental records only reached  ≈ 0.4 m due to technical difficulties with the setup[37].As Fig.3illustrates, the general features of the  −  records from the simulations and experiments are similar.Both show several distinct load peaks, here referred to as peak load events.The breakage model appears to improve the simulation results, as the simulations without it yield several load peaks of considerably higher magnitude than those in the experiments.It is worth emphasizing here that the simulated and experimental  − records do not show one-to-one agreement, since the initial rubble configurations were not the same.Nevertheless, records from the simulations with the breakage model are more representative of experimental curves with regard to the peak load magnitudes.The  −  records, however, from the simulations without the breakage model initially follow those from the simulations without it.The two records deviate from each other as the shear force, , builds up towards the first peak load; clearly breakage affects the mechanism that limits peak loads, as reflected further by the magnitudes of them being in general smaller in simulations with the breakage model than in those without it.Fig.4comparesthe data from the simulations and experiments further by showing the mean and maximum of the shear force,  and   , respectively.Data for confining pressure  = 5.75 and 11.03 kPa are included.Since in both the experiments and simulations the details of the  − records for a given parameterization depended on the initial arrangement of the particles, the figure shows data collected from five repeated simulations for each case.Meanwhile, the experimental data are from three repeated experiments.For each simulation and experiment,  was defined by simply taking the average of , and   as the maximum value of .The results are for the displacement interval  = 0 … 0.3 m, since the effect of shear box wall on  is more pronounced for larger .Fig.4showsthat, overall, the simulations allowing breakage yield  and   data more closely resembling that from the experiments than the simulations without it. and   in the simulations and experiments differed on average by 11% and 14%, respectively, with the breakage model.Simulations without the breakage model led to errors of 23% and 41% respectively.Thus, the values for  and   are more close

Fig. 3 .
Fig.3.Typical shear force-displacement,  − , plots from shear box simulations and experiments[37] with (a) small and (b) large particles.The confining pressure,  , in the simulations and the experiments was 5.75 kPa.

Fig. 4 .
Fig. 4. Mean shear force, , in simulations and experiments in the case of (a) small and (b) large particles with confining pressures of  = 5.75 and 11.03 kPa.Maximum shear force,   , in the case of (c) small and (d) large particles.The figure shows data points from five repeated simulations for each case. and   are for the shear box displacement interval 0.0 … 0.3 m.

Fig. 5 .
Fig. 5. Snapshots from the simulations, without the breakage model (left column) and with the breakage model (right column) with the shear force-displacement,  − , records in Fig. 3a.Snapshots are from the first peak load event.First (a), in both simulations, a force chain forms near the shear plane.Then, (b) particles within the chain break when breakage is allowed (locations of particle breakage are circled with dashed line).After this, (c) particles within the force chain are further compressed in the simulation without the breakage model, while the force chain has disappeared and new fragments have formed in the simulation with the breakage model.

Fig. 6 .
Fig. 6.Effect of particle strength on the shear strength of ice rubble, with two different confining pressures,  : mean shear resistance,   , against the strength of particles in the case of (a) small particles and (b) large particles.Maximum shear strength,    , against the strength of particles in the case of (c) small particles and (d) large particles.Strength was varied by changing the value of internal cohesion,  [Eq.(4)].Data from the simulations without the breakage model ( = ∞ kPa) are also presented.The figure shows data points from five repeated simulations for each case, their mean and standard deviation, and a curve fitted to the mean values (solid line).

Fig. 7 .
Fig. 7. (a) Normalized mean ice rubble shear strength,   and (b) normalized maximum ice rubble shear strength,    , values against the internal cohesion, , of the particles.Here, values were normalized by the mean   and mean    of the corresponding simulations without the breakage model.Solid lines in the figure are curves fitted to the mean values.

Fig. 8 .
Fig. 8. Mean shear resistance,   , against the confinement,  , as observed for shear boxes with lengths, (a)  = 0.6 and (b)  = 3.0 m, that is ∕ = 20 and 100, respectively.Maximum shear strength,    , for shear boxes with lengths, (c)  = 0.6 and (d)  = 3.0 m.The figure shows data points from five repeated simulations for each case, their mean and standard deviation, and Mohr-Coulomb material model [linear fit of Eq. (6)] fitted to the mean values.

Fig. 9 .
Fig. 9. (a) Internal cohesion,   , and (b) friction angle,   , of ice rubble deduced from the numerical direct shear box experiments against the ratio of box length to particle length, ∕.Solid and dashed lines in the figure are curves fitted to the results from mean shear resistance and maximum shear strength,   and    , respectively.

Table 1
[35,36]mulation parameters, chosen based on experiments in Pustogvar et al.[37]where applicable.The breakage model parameters are from Prasanna et al.[35,36].Commonly, direct shear box experiments are used to determine parameters for the Mohr-Coulomb material model,   and   = tan   , where   is the internal cohesion of rubble and   is the angle of internal friction, often reported for granular materials.These are determined by conducting experiments with various values of  and fitting  −   results to -Width × height  × ℎ 30 × 20, 60 × 40 mm × mmresistance and the maximum shear strength of the rubble, τ and    , respectively.Mean shear resistance is here defined by taking the average of   () for a given interval of .