Numerical simulation of the effect of grain fragmentation on the evolution of microstructure quantities

In this paper the effect of grain fragmentation of a cohesionless granular material on the change of microstructure quantities is investigated using a micropolar continuum model. To this end the change of the grain size distribution is related in a simplified manner to a reduction of the mean grain diameter, which enters the constitutive equation as the internal length. The additional densification of the grain skeleton is modelled by a reduction of the incremental stiffness, which is related to the so-called solid hardness defined in the sense of a continuum description. The reduction of the mean grain diameter and the solid hardness caused by grain fragmentation is described by corresponding evolution equations, which depend on an increasing mean pressure and an increasing deviatoric stress. The focus of the numerical investigations is on the evolution of microstructure quantities like the void ratio, mean grain diameter and couple stresses under monotonic compression. It is shown that for an initially inhomogeneous distribution of the void ratio polar quantities like microrotations and couple stresses can develop even outside of shear bands. With increasing compression the polar quantities can increase, while the range of fluctuation of the void ratio decreases.


Introduction
Grain fragmentation within a stressed granular body changes the grain size distribution and is usually also accompanied by a reduction of the void space caused by a reorientation of particles in the grain skeleton. The mechanisms of grain fragmentation are complex and various concepts for analysing experimental investigations and mathematical modelling are proposed in the literature, e.g., [2, 16, 20, 23, 24, 30, 43-45, 51, 53]. Experimental investigations reveal that the evolution of grain fragmentation strongly depends on the morphology of the granular material, the grain strength and the loading path, e.g., [3,18,32,39,[47][48][49][50].
The focus of the present paper is on modelling the influence of grain fragmentation using a hypoplastic material description. The concept of hypoplasticity of the Kolymbas type was originally developed within a non-polar continuum [37,38]. It differs fundamentally from the concept of elasto-plasticity as no decomposition of the deformation into elastic and plastic parts is needed to model non-linear and inelastic material properties. In the simplest case the constitutive relation for the stress rate depends on the current E. Bauer (&) Á S. Safikhani Á L. Li Institute of Applied Mechanics, Graz University of Technology, 8010 Graz, Austria e-mail: erich.bauer@tugraz.at stress state and strain rate. More enhanced hypoplastic models include additional state variables such as for instance the pressure dependent relative density, e.g., [5,26], and an inter-granular strain [46]. As the hypoplastic constitutive equation is incrementally nonlinear, different incremental stiffnesses for loading and unloading are modelled. Extensions of the concept of non-polar hypoplasticity to a micropolar continuum are, for instance, proposed by Tejchman and Bauer [52] and Huang and Bauer [36]. In the micropolar model a constitutive equation for the couple stress rate is needed in addition to the constitutive equation for the stress rate. When shear strain localization takes place, the benefit of a micropolar model is twofold: the results obtained from finite element simulations are independent of the element size, e.g., [17,52], and the internal length of the constitutive equation can be related to the mean grain size, which also triggers the thickness of the shear bands predicted, e.g., [36]. In the model proposed in the present paper the effect of grain fragmentation on the change in the grain size distribution is taken into account in a simplified manner by a reduction of the corresponding mean grain diameter. The additional densification of the stressed granular body is modelled by a reduction of the incremental stiffness. Herein the so-called solid hardness is a key parameter to reflect the influence of grain fragmentation on the additional densification. The solid hardness is defined in the sense of a continuum description and it is a parameter in the compression law by Bauer [4,5]. This compression law allows the simulation of the effect of grain fragmentation on the mechanical response under isotropic loading for an unlimited pressure range. The concept of the solid hardness was embedded into non-polar constitutive models, e.g., [5,7,8,26,31,46,54], into micropolar models, e.g., [13,21,27,28,35,36,52], and extended to moisture sensitive and weathered rockfill materials, e.g., [9,12,41].
In order to reflect the fact that grain fragmentation not only depends on the pressure level a more general concept was proposed to model the evolution of the solid hardness depending on the loading path [14,15]. The focus of the present investigation is on the effect of an increase in the mean pressure and the deviatoric stress. To this end specific evolution equations for the reduction of the solid hardness and the mean grain diameter are considered and embedded into a micropolar hypoplastic model by Huang et al. [35].
The present paper is organized as follows: Sect. 2 gives the definition of the solid hardness in the sense of a continuum description. To allow the calibration of the compression law by Bauer [4,5] to experimental data for a low pressure range a fictive solid hardness is introduced in Sect. 3. In Sect. 4 the evolution equation for the reduction of the fictive solid hardness is extended to model the effects of both increasing mean pressure and deviatoric stress. A similar evolution equation is introduced to describe the reduction of the mean grain diameter. Section 5 outlines the micropolar hypoplastic model which is used for the numerical simulations in Sect. 6. The focus of the numerical simulations is on the evolution of microstructure quantities like the void ratio, mean grain size, microrotation and couple stresses under monotonic plane strain compression. For an initial random distribution of the initial void ratio the numerical results are discussed for oedometric compression and biaxial compression in Sects. 6.1 and 6.2, respectively.
Throughout the paper indices on vector and tensor components refer to an orthonormal Cartesian basis. The summation convention over repeated indices is employed. A superimposed dot indicates the material time derivative, and compressive stress and strain and their rates are negative as in the sign convention of rational continuum mechanics.
2 Definition of the solid hardness in the sense of a continuum description For a granular material in the initially loosest state, the reduction of the corresponding maximum void ratio e i under increasing isotropic pressure p ¼ À r ii =3 can be modelled by the following compression law proposed by Bauer [4,5]: Equation (1) describes the pressure dependent upper bound of possible void ratios, i.e., there exists no grain skeleton with viable grain contacts beyond e i ðpÞ. The compression law includes three parameters. Parameter e i0 denotes the maximum void ratio for p % 0, parameter h s is called ''solid hardness'' and defined for the pressure 3 p, where the compression curve in a semi-logarithmic representation shows the point of inflection, and parameter n is related to the inclination of the compression curve at the point of inflection (Fig. 1). For p ! 1 the value of the void ratio asymptotically tends towards e ! 0. Thus, the compression law (1) is applicable to a wide range of pressures and can model the history of grain rearrangement and grain fragmentation within the whole pressure range in a phenomenological manner. The introduction of the term ''solid hardness'' for the parameter h s was motivated by the results obtained from high pressure isotropic or oedometric compression experiments [5] and from numerical simulations with the discrete element method [40] using an assembly of breakable particles. The data obtained show that the amount of grain crushing becomes extremal for pressures close to 3 p ¼ h s . In this context it is important to note that the term solid hardness is related to the compression behaviour of an aggregate of grains and should not be confused with the hardness of a single grain.

Concept of the fictive solid hardness
Laboratory experiments carried out for instance with quartz sand show that the point of inflection of the compression curve illustrated in Fig. 1 can only be reached with laboratory devices appropriate for high pressures of p % 80 MPa. If experimental data from compression tests are only available for a lower pressure range, the calibration of Eq. (1) will lead to an overestimation of the value of the solid hardness [31]. This is justified by the fact that under lower isotropic pressures the solid grains do not reveal their intrinsic ultimate strength, i.e., under low pressures the grains almost behave like rigid bodies and the densification caused by the reorientation of the grain skeleton is small. Therefore, the calibration using data only available for a rather low pressure range leads to a fictive solid hardness, which is higher than the corresponding solid hardness h s defined in Fig. 1. This fictive solid hardness obtained can be used to simulate the compression behavior for the corresponding restricted pressure range.
For the extended concept of the solid hardness introduced in the next section it is convenient to model a smooth reduction of the fictive solid hardness with an increase of the pressure range. To this end the constant solid hardness h s in Eq. (1) is replaced by the pressure dependent quantity h sp . The following evolution equation for the fictive solid hardness was proposed [15]: Herein _ p denotes the change of the mean pressure with time, i.e., _ p ¼ dp=dt, and b p and n p are material parameters. Equation (2) describes a rate independent behaviour, and with respect to the reference state h sp ðp % 0Þ ¼ h s0 the integration yields for the decrease of the fictive solid hardness h sp with an increase of the mean pressure p: 4 Decrease of the solid hardness and the mean grain diameter with increasing mean pressure and deviatoric stress In this section a model is proposed to describe the effect of grain fragmentation on the reduction of the mean grain size and the additional densification caused by a reduction of the solid hardness. Experiments show that the amount of grain damage does not only depend on the pressure level but it is also influenced by the loading path, e.g., [3,29,33]. Thus, for a refined modelling of the reduction of the solid hardness and the mean grain diameter the corresponding evolution equations should take into account the essential mechanisms of grain fragmentation depending on the loading path as discussed for instance in [13][14][15]. The focus of the present paper is on the constitutive modelling of the effect of grain fragmentation caused by an increase of the mean pressure and of the deviatoric stress. To distinguish these two effects the following evolution equations for the solid hardness h st and the mean grain diameter d 50 are proposed: The structure of Eqs. (4) and (5) is similar. The first term on the right-hand side of these equations takes into account the influence of an increase of the mean pressures as outlined in the previous section while the second term is related to an increase of the deviatoric stress. As the Cauchy stress tensor in a micropolar continuum is usually non-symmetric, in Eqs. (4) and (5) only the symmetric part of the stress deviator is considered to avoid non-real eigenvalues, i.e., the components of the symmetric part of the stress deviator are: where d ij denotes the Kronecker delta. As grain fragmentation is irreversible, a reduction of h st and d 50 takes place only for _ h st \0 and _ d 50 \0. In both functions the first term only leads to a contribution for an increasing mean pressure, i.e., _ p [ 0, while the second term takes into account an increase of the deviatoric stress. This distinction also allows an easier calibration of the material parameters involved. In particular the parameters b p and b d in Eq. (4) are related to the additional compaction depending on the mean stress and diviatoric stress, while parameters b dp and b dd in Eq. (5) can be calibrated based on the change of the grain size distribution caused by isotropic and deviatoric stress paths.

Micropolar hypoplastic model
In a micropolar continuum, e.g., [19,22,25], a point possesses displacement degrees of freedom u i (i ¼ 1; 2; 3) and rotational degrees of freedom called microrotations x c i (i ¼ 1; 2; 3). The rate of deformation and the rate of curvature are defined as [36]: respectively. Herein kij denotes the permutation symbol. In Eq. (6) the velocity gradient can be decomposed into the symmetric part and the skew-symmetric part Hence the rate of deformation _ e c ij can alternatively be represented as: where W ij denotes the macro-spin and W c ij ¼ À kij _ x c k denotes the micro-spin or so-called Cosserat spin. Equation (8) indicates that in the case where the macro-spin is equal to the micro-spin the rate of deformation tensor reduces to the stretching tensor of the standard non-polar continuum. The kinematic quantities are associated with the stress tensor r and the couple stress tensor l defined in the current configuration. For a micropolar continuum the local equilibrium equations read: Herein b and c represent the specific body force and the specific body couple, respectively. It follows from Eq. (10) that the stress tensor is non-symmetric with the exception of states with ol ij =ox j ¼ 0 and c i ¼ 0.
The evolution of the void ratio e is based on the assumption that the density of the solid grains remains constant. Then a relation between the rate of the void ratio _ e and the volume strain rate _ e v ¼ _ e c kk ¼ _ e kk can be derived from the balance equation of mass, which leads to: As the volume of the solid material is assumed to be constant Eq. (11) also holds when grain fragmentation takes place. For modelling the evolution of the stress r ij and the couple stress l ij the concept of hypoplasticity is used, where the objective stress rate r ij and couple stress rate l ij are functions of the current void ratio e, the mean grain diameter d 50 , the solid hardness h st , the rate of deformation _ e c ij , the normalized quantitiesr ij of the nonsymmetric stress tensor, the normalized componentsl ij of the couple stress tensor and the normalized rate of curvature _ j ij . In this paper a particular micropolar hypoplastic version by Huang et al. [35] is used, where the constant solid hardness and mean grain diameter in the original constitutive relations is replaced by the state dependent quantities defined in Eqs. (4) and (5), respectively. Then the constitutive relations for the components of the objective stress rate and couple stress rate read: For the relation between the objective stress rate, the objective couple stress rate and its time derivatives the formula given by Zaremba and Jaumann is used, i.e., In Eqs. (12, 13) the following normalized quantities are used: According to the general concept of hypoplasticity the constitutive Eqs. (12) and (13) are incrementally nonlinear, and thus, different stiffnesses for loading and unloading are modelled without a decomposition of the deformation into elastic and plastic parts. Factorsâ and a m in Eqs. (12) and (13) are related to stationary states or so-called critical states. In particular, the micropolar parameter a m is assumed to be a constant, [36], and for granular materialsâ can be adapted to the stress limit condition proposed by Matsuoka and Nakai [42]. Then functionâ depends on the so-called critical friction angle u c and the symmetric part of the normalized deviatoric stress, r ds ij , which can be represented as [7,10]: The influence of the current density and the pressure level on the incremental stiffness, the peak friction angle and dilatancy behavior is modelled with the stiffness factor f s and the density factor f d . In particular the density factor f d represents a relation between the current void ratio e, the critical void ratio e c and the minimum void ratio e d , i.e., where a is a material parameter. The critical void ratio e c and the minimum void ratio e d are pressure dependent and related to the maximum void ratio e i according to the postulate by Gudehus [26], i.e., where e i0 , e d0 , e c0 are the corresponding values for r kk % 0. With respect to the evolution Eq. (4) a decrease of the solid hardness h st also leads to a decrease of the limit void ratios and the critical void ratio.
The stiffness factor f s is the product of three parts, i.e., Herein the density dependent part f e is defined as the ratio of the the maximum void ratio to the current void ratio e, i.e., where b is a material parameter. f r takes into account a decrease of the incremental stiffness with an increase of the deviatoric stress, i.e., f b is called barotropy factor and obtained by a consistency condition which enables the embedding of the compression law (1) into the constitutive equation for the stress rate. The consistency condition required for the case of a constant solid hardness was first discussed by Gudehus [26]. In the present model the solid hardness changes according to the evolution Eq. (4), which must be taken into account for the derivation of the consistency condition. For consistency the response of the constitutive Eq. (12) for a monotonic isotropic compression starting from the maximum void ratio e i0 must coincide with Eq. (1). For the isotropic compression of an material element no polar terms appear in Eq. (12) and therefore the consistency condition leads to: with It can be noted that for the special case of coaxial and homogeneous element deformations starting from an initially symmetric stress tensor there are no polar effects and the micropolar hypoplastic model reduces to the non-polar version given by Gudehus [26] and Bauer [5], i.e., It can be concluded that with the exception of the polar constant a m and the evolution of d 50 all other material parameters of the micropolar hypoplastic model are the same as in the non-polar version of Eq. (23). Thus homogeneous element deformations are sufficient for the calibration of these material parameters. For instance the parameter h s0 , n, e i0 , e c0 , e d0 , a, b and u c can be calibrated based on data from oedometric and triaxial compression tests, and index tests as outlined in [5,31]. Parameters b p , n p and b d of the evolution equation for h st are related to grain fragmentation under isotropic and deviatoric tests and can be obtained by back calculation. Parameters b dp , b dd and n dp of the evolution equation for d 50 are related to the change of the grain size distribution curve. For the numerical simulations discussed in the next section the material parameters used are summarized in Table 1.

Results of numerical simulations
For the numerical simulations the micropolar hypoplastic model is implemented in the finite element code ABAQUS [1] and for plane strain conditions a four-node element with bilinear shape functions developed by Huang [34] was used to describe the displacements and Cosserat rotations within the element. The non-zero kinematic and static quantities relevant to plane strain conditions are illustrated in Fig. 3.
In the following the numerical results obtained from finite element simulations of plane strain oedometric and biaxial compression tests are discussed. The granular specimen with an initial width of b 0 ¼ 4 cm and a height of h 0 ¼ 12 cm is discretized by 3072 Cosserat elements, i.e., the initial size of each finite element is approximately 1:1mm Â 1:1mm. The corresponding boundary conditions are illustrated in Fig. 4. Along the whole boundary of the specimen the nodes can freely rotate. In both tests the compression is modelled by the vertical displacement s prescribed for the nodes at the top surface.

Plane strain oedometric compression
Plane strain oedometric compression tests are investigated for two different random distributions of the initial void ratio e 0 . The ranges of initial void ratios assumed are outlined in Fig. 5. The initial ranges are 0:55\e 0 \0:75 for the dense specimen and 0:92\e 0 \1:16 for the loose specimen. Figure 6 shows for a particular element that the compressibility of the initially dense specimen is much smaller than the one for the loose specimen. With an increase of the pressure level the distance between the two curves becomes smaller and for very high pressures close to the point of inflection the curves coincide. This indicates that under monotonic compression the memory of the material of the initial state is erased as a result of grain crushing and reorientation of the grain skeleton into a denser state. This is a general property of the proposed model and is also confirmed   Fig. 5 Histogram of the random distribution of the void ratio for an initially dense and loose specimen by experiments carried out with various granular materials and different initial densities. It can therefore be concluded that for monotonic compression the point of inflection is a well defined state and it is therefore appropriate to introduce the term solid hardness as a key parameter for the pressure level where grain crushing becomes dominant. As a result of the random distribution of the initial void ratio the change of the void ratio with increasing compression is not the same in all elements. Figure 7a and b shows the fluctuation of e over the cross section in the middle of the specimen for both initially dense and initially loose specimens. With respect to the assumed initial random distribution of the void ratio, the range of the fluctuation, De, is larger for the initially loose specimen. With monotonic increasing pressure a continued decrease of De can be observed. For higher compressions De tends towards zero for both initially dense and initially loose specimens as shown in Fig. 7c. This result again indicates that the memory of the material of the initial density is erased due to grain fragmentation under increasing stress.
The response of the constitutive relations for r ij and l ij is also influenced by a reduction of the mean grain diameter d 50 . The reduction of the mean grain diameter is larger for the initially dense material as shown in Fig. 8. According to the evolution Eq. (5) the value of d 50 is the same for the same stress state. However, the random distribution of the void ratio leads to an inhomogeneous evolution of the stresses within the specimen and consequently to a fluctuation of d 50 . Numerical simulations show that the fluctuation increases with increasing vertical compression. Although no microrotations are locked along the boundary of the specimen, polar quantities like couple stresses and microrotation can develop during compression. This is caused by the inhomogeneous distribution of stresses and deformation within the specimen. But no clear tendency of the evolution of couple stresses can be detected from the numerical simulations as shown for instance in Fig. 9 for three different elements.

Biaxial compression under constant lateral stress
To investigate plane strain biaxial compression under constant lateral stress the same random distribution of the initial void ratio is used as for the initially dense specimen in the oedometric compression test. Starting from an initial stress of À 100 Pa the vertical and horizontal stress is increased up to À 0:5 MPa. Then under the constant lateral stress of À 0:5 MPa the specimen is vertically compressed by a vertical displacement s prescribed for the nodes of the top surface (Fig. 4b). In contrast to oedometric compression a pronounced shear strain localization appears in the biaxial compression test, which leads to a strong inhomogeneous deformation. The shear band is clearly visible for a vertical compression of s ¼ 1:3 cm in Fig. 11. Figure 10 shows the development of the stress ratio, r 22 =r 11 , and volume strain, e v , for an element inside the shear band and an element outside the shear band until a vertical displacement of s ¼ 1:8 cm. It seems to be that the state where the curves bifurcate is a precursor of the onset of shear strain localization. However, the state the curves bifurcate is far before the peak states and thus contrary to what can be expected from experiments and the theoretical bifurcation analysis, i.e., the theoretical analysis usually predicts the state of the first possible bifurcation much closer to the stress peak, e.g., [6,11]. Therefore, the state the curves bifurcate is not coincide with the state where shear banding develops. This discrepancy between analytical investigation and numerical simulation may be explained by the fact that the initial random distribution of the void ratio also influences the evolution of the stresses within the specimen. In this context it can be noted that although the stress prescribed along the lateral boundaries of the specimen is constant, the horizontal stress r 11 within e 3p [MPa] Fig. 6 Reduction of the void ratio e with increasing mean pressure p the specimen is inhomogeneously distributed from the beginning of compression. This also causes the evolution of polar quantities in the specimen before the shear band fully develops. Figure 10a shows that after the curves bifurcate the increase of the stress ratio inside the shear band is larger than outside. This leads to a more pronounced reduction of the mean grain diameter d 50 inside the shear band (Fig. 11c), which has also been reported from experiments, e.g., [2]. At the beginning of loading the volume strain e v in both elements shows compaction, but with advanced vertical compression pronounced dilatancy only appears within the shear band (Fig. 10b). In the element outside the shear band the volume strain reaches an almost stationary value close to the stress peak. The pronounced dilatancy within the shear band is also visible in the contour plot of the void ratio in Fig. 11b. The deformation inside the shear band is also accompanied by pronounced microrotations, which are clearly visible in Fig. 11d. (a) (b) (c) Fig. 7 Fluctuation of the void ratio e in the vertical section in the middle of the specimen: a initially dense specimen; b initially loose specimen; c reduction of the range De against the vertical displacement s Fig. 8 Reduction of the mean grain diameter d 50 against the vertical displacement s Fig. 9 Evolution of couple stresses l 32 within three different elements against the mean pressure p

Conclusions
For modelling the influence of grain fragmentation on the mechanical response of a stressed granular body a micropolar continuum description is used. Herein the effect of a change in the grain size distribution is modelled in a simplified manner by a reduction of the corresponding mean grain diameter, which enters the constitutive equation as the internal length. Evolution equations are proposed to describe the reduction of the mean grain diameter and the incremental stiffness depending on increasing mean pressure and increasing deviatoric stress. Herein the solid hardness is defined in the sense of a continuum description and serves as a key parameter for modelling the densification of the granular material caused by grain fragmentation. The mean grain diameter and the solid hardness are embedded into a micropolar hypoplastic model as state variables. Numerical simulations carried out with an initial random distribution of the void ratio show that even under oedometric compression polar quantities like micropolar rotations and couple stresses can develop. With an increase of the stress level the memory of the initially inhomogeneous distribution of the void ratio is gradually erased, i.e., the amount of (c) (d) Fig. 11 Contour plot of the a random distribution of the void ratio in the undeformed specimen, b-d void ratio e, micropolar rotation x c 3 and mean grain diameter d 50 for a lateral stress of À 0:5 MPa and a vertical displacement of 1.3 cm fluctuation decreases with an increasing stress level. A similar behavior can also be detected under biaxial compression. When shear strain localization becomes dominant, the reduction of the mean grain diameter is more pronounced within the shear band. Although, the concept proposed for modelling the influence of grain fragmentation on the mechanical response is mathematically simple, further experimental and physical investigations are needed for a better interpretation of the material parameters involved.