Analysis of cement-treated clay behavior by micromechanical approach

Experimental results show the significant influence of cement content on the mechanical properties of cement-treated clays. Cementation is produced by mixing a certain amount of cement with the saturated clay. The purpose of this paper is to model the cementation effect on the mechanical behavior of cement-treated clay. A micromechanical stress-strain model is developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding and debonding during mechanical loading are introduced in two ways: an additional cohesion in the shear sliding and a higher yield stress in normal compression. The model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Ariake clay and Singapore clay with various cement contents. The applicability of the present model is evaluated through comparisons between numerical and experimental results. The evolution of local stresses and local strains in inter-cluster planes are discussed in order to explain the induced anisotropy due to debonding at contact level under the applied loads. The numerical simulations demonstrate that the proposed micromechanical approach is well adapted for taking into account the main physical properties of cement-treated clay, including damage and induced anisotropy under mechanical loading.


Introduction
Stabilization of soft ground by deep cement mixing and jet grouting methods, which both introduce cement into the ground, has been widely used in geotechnical projects [1][2][3][4]. Due to the increase of their mechanical properties, the cement-treated soil layer helps to decrease the soil movements during construction. Thus, proper understanding of the cementation influence on the mechanical behavior of cement-treated clay is of great importance.
The improvement of the engineering properties of cement-treated soil has been attributed to the soil-cement reaction [5,6]: Hydration reaction and pozzolanic reaction between the hydrated lime and the clay minerals. This creates a type of microstructure different from the natural structured soils (soft sensitive clay and stiff clay [7,8]), for which soil structures are often formed by solute deposition at inter-particle contacts, charge deficiencies, and Van der Waal's forces [9].
Experimental results show an increase of shear strength and shear modulus with cement content, accompanied by a large post-peak reduction in the strength of the treated soil [3][4][5]10,11]. Cementation also makes the yield stress and swelling index higher than those of untreated soil. Bond degradation is often observed if the compression load exceeds the yield stress, as it has also been observed in natural structured soils [12][13][14]. The difference between natural structured clay and cement-treated clay lies in the properties of host clay modified by adding cement, which needs to be considered in constitutive modeling.
The behavior of cement-treated clay has been modeled by conventional elasto-plastic models, e.g., Lee et al. [15], Horpibulsuk et al. [16]. In their approaches, the material structure, including fabric and bonding, was accounted for by the initial size of the yield surface. A damage-type mechanism was introduced which produced the reduction of the yield surface size due to bond degradation. However, the bonding and debonding at inter-cluster contact level has not been explicitly considered in these models.
To tackle this problem, a micromechanical modeling appears to be a suitable tool. The approach is based on a homogenization technique in order to derive the stressstrain relationship of the granular assembly from forces and displacements at the particle level [17][18][19][20][21][22][23][24][25][26][27][28][29]. Along these lines, Chang and Hicher [27] have proposed a model for granular materials based on the development of micromacro links between strain and inter-particle displacements on one hand, and stress and inter-particle forces on the other hand, which allow both the micro-system (intergranular forces and displacements) and the macro-system (overall stress-strain) to be in balance. More recently, the model has then been extended to clays by considering the clayey microstructure as a collection of clusters which play the same role as the grains in granular assemblies for various features [7,8,[30][31][32][33].
Horpibulsuk et al. [3] have shown that the clay fabric consists of clay cluster (aggregation of clay particles) and that the role of the induced cementation is to weld clay clusters (see Fig. 1(a)). Kamruzzaman et al. [4] have made SEM photos of untreated and cement treated Singapore clay with various cement contents under various consolidation pressures (see Fig. 1(b)), which demonstrate that the size of the clay clusters (as shown by dash circle in Fig. 1(b)) can be considered the same for all cases. Therefore, the cement-treated clay can be seen as an assembly of clusters (aggregates of clay-cement particles). The deformation of an assembly can be obtained by integrating the movement of the inter-cluster contacts in all orientations. Thus, the effect of the inter-cluster bonding can be explicitly taken into account. As the orientation dependent properties are explicitly introduced, the induced anisotropy due to inter-cluster debonding can be modeled in a direct way.
In this study, the cementation effect on the mechanical properties of cement-treated clays is first investigated. Based on this, the dependency of the inter-cluster bonding and of the critical state line (CSL) location on the cementcontent are implemented into the model developed by Yin et al. [7] for sensitive clay. The extended model is then used to simulate isotropic compression and undrained triaxial tests under various confining stresses on cementtreated clays with different cement contents. The evolution of local stresses and strains in inter-cluster planes due to externally applied loads is also discussed in order to describe the mechanism of induced anisotropy due to debonding at the micro contact level. The overall applicability of the present model is evaluated based on the comparison between experimental results and numerical simulations.

Influence of cement content on mechanical properties
In this section, the cement-treated Ariake clay with cement content varying from 0% to 18% (0% representing untreated clay [3]) and the cement-treated Singapore clay with cement content varying from 0% to 50% [4] were selected to investigate the influence of cement content on clay mechanical properties. The percentage of cement content was calculated based on the ratio of dry weight of cement to dry weight of soil. Some physical properties of the untreated clays are summarized in Table 1. The type I ordinary Portland cement was used for treating both clays.
The initial water content of the clay slurry may influence the initial void ratio, the compression yield stress, and the destructuring rate during mechanical loading for samples of cement-treated clay. In this paper, the two clay slurry were prepared at a water content equal to 1.4-1.5 times the liquid limit (180% for Ariake clay, 120% for Singapore clay), and the influence of the initial water content was not studied. Furthermore, all the experimental results selected for this study were obtained on samples of cement-treated clay with a curing time of 28 days. The influence of curing time on the cementation was not studied at current stage.

Cement-content dependent behavior under isotropic and 1D compression
Isotropic compression and oedometer tests by Horpibulsuk et al. [3] and Kamruzzaman et al. [4] on two cement-treated clays show the significant influence of the cement content on the swelling index (see Fig. 2(a)) and on the compression yield stress (see Fig. 2(b)). This influence can be expressed by Eq. (1) for the swelling index κ and by Eq. (2) for the bonding pressure p b defined as the difference of the compression yield stress between cement-treated clay p yc and untreated clay p yu : p b = p ycp yu , as follows: where κ 0 is the swelling index for untreated clay; c is the cement content; β k , β p1 and β p2 are material constants (see Fig. 2(a), (b)). By comparing the post-yield stress-strain curves of cement-treated samples to those of untreated samples [3,4], a mechanical bond degradation process can be observed. This bond degradation is an irreversible phenomenon that,  experimentally, appears to be controlled by plastic strain accumulation [34]. The bonding and debonding can be linked directly to the inter-cluster bonds due to the cementation formed when adding cement into clay slurry.

Cement-content dependent shear behavior
The friction angle at failure was measured based on undrained triaxial tests by Horpibulsuk et al. [3] and Kamruzzaman et al. [4] on both cement-treated clays by extending the curves to an axial strain up to 15% according to the previous trend. As shown in Fig. 2(c), for Ariake clay, having an initially high friction angle, the cement content influences slightly the friction angle at failure (we obtained f μ = 40°measured from undrained stress paths, which is slightly different from f μ = 38°provided by Horpibulsuk et al. [3]). However, for Singapore clay with an initially small friction angle, the cement content influences significantly the friction angle at failure. This influence can be expressed by Eq. (3) as follows: where f μ0 is the friction angle at failure for the untreated clay; f μmax is the maximum friction angle at failure for the cement-treated clay; β f is a material constant. The apparent cohesion due to cement bonds is measured from the position of the failure line in the p'q plane. Figure 2(d) shows the significant influence of the cement content on the apparent cohesion c b for both cementtreated clays. This influence can be expressed by Eq. (4), having a form similar to Eq. (2) for the bonding pressure, as follows: where β c1 and β c2 are material constants. The stress-strain curves of the two cement-treated clays show also the influence of the cement content on the initial shear modulus, which is similar to its influence on the swelling index (or bulk modulus). Since the shear modulus cannot be accurately measured based on digitized data, we will not discuss it in this paper.
Upon increase of the applied shear stresses, the intercluster bonds are progressively damaged. As a result, a decrease of the deviatoric stress takes place after the peak for cement-treated samples.

Cement-content dependent critical state line (CSL)
The critical state lines are obtained from the undrained triaxial tests for untreated clays, as shown in Figs. 3(a) and (b). The void ratio corresponding to this state is e c . The CSL can be written as follows for clay: where p' is the mean effective stress; l is the slope of the critical state line in the e-logp' plane, e cr is the reference void ratio at critical state corresponding to p ref = 100 kPa. For cement-treated clays, due to the strain localization of samples during shearing, the critical state line is difficult to be obtained accurately. If the failure state in the p'-q plane (corresponding to 15% of axial strain) is considered as the critical state, the critical state lines in the e-logp' plane can be obtained for all cement-treated samples. For cementtreated Ariake clay (see Fig. 3(a)), the critical state line can be assumed to be parallel to that of the untreated clay. If this assumption is also adopted for cement-treated Singapore clay, the critical state lines can be obtained as shown in Fig. 3(b). Therefore, the influence of cement content on the location of the critical state line can be obtained by making the reference void ratio function of the cement content, as shown in Fig. 3(c), expressed as follows where e cr0 is the reference critical state void ratio for the untreated clay; e crmax is the maximum reference critical state void ratio for cement-treated clay; β e is a material constant.

Micromechanical model
The model developed in this study is an extension of the model developed by Yin et al. [7] for sensitive clay. The modifications made to account for the effect of cement content on the inter-cluster bonding and debonding and on the CSL location are discussed below. As stated in the previous section, the bonding and debonding will be considered for two different loading conditions: normal compression and shear sliding.

Inter-cluster behavior
To take into account the effect of the inter-cluster bonding, the plastic law of the shear sliding at inter-cluster contacts proposed by Yin et al. [7] is adopted. The yield function for shear sliding is assumed to be of Mohr-Coulomb type, written in a contact plane (e.g., , τ s , τ t defined in Appendix A, see coordinate system in Fig. 4) as follows where c b is the inter-cluster cohesion due to the intercluster bonding, H 1 ðγ p Þ is a hardening/softening function. The values of c b and f μ depend upon the cement content (see Eqs. (3) and (4)). The influence of the cement content on the shear behavior is thus introduced. For untreated clay, c b = 0, and the equation can be reduced to that used by Yin et al. [33]. The local dilatancy equation is derived from that proposed by Yin et al. [33], as follows D is a material constant which controls the dilatancy rate during shearing. The void ratio at critical state e c is calculated using Eq. (5) where the influence of the cement content is introduced by e cr (see Eq. (6)). The degradation of the inter-cluster bonds can be modeled as a damage of the bonded contacts. Therefore, the damage law proposed by Yin et al. [7] is adopted in the expression of the inter-cluster cohesion, as follows where c b0 is the initial inter-cluster cohesion; d is the factor of damage representing the influence of the tangential displacement in the damage law. Therefore, during loading each contact sliding produces a progressive damage of the cohesion.
To account for the inter-cluster bonding effect on the compressive behavior of cement-treated clay, we adopt the plastic law for the normal compression of inter-cluster contacts proposed by Yin et al. [7] where p ¼ pi ð1 þ χÞ is the compression yield stress; pi is the intrinsic compression yield stress corresponding to the yield stress measured from an isotropic compression test on a reconstituted sample of untreated clay; χ is the bonding ratio defined by χ ¼ p = pi -1. The value of p depends on the cement content ( p ¼ pi þ p b , see Eq. (2) for p b ). Thus, the influence of the cement content on the compression behavior is introduced. The hardening function controlling the evolution of pi is defined as follows where c p is the compression index determined from the compression curve plotted in the ε plog plane.
The damage law proposed by Yin et al. [7] is also adopted in the expression of the bonding ratio at intercluster level, where χ 0 is the initial bonding ratio depending of the cement content; n is the factor of damage representing the influence of the normal displacement in the damage law. Therefore, during loading each contact produces a progressive damage of the bonding ratio. For untreated clay χ 0 = 0 (or p ¼ pi ), the equation can be reduced to that used by Yin et al. [33].

Study of cementation related parameters
The parametric study of the cementation related parameters of the model is carried out under different loading conditions using the model constants given in Table 2 for the cement-treated Ariake clay. First, an isotropic compression test with p b = 25 kPa is considered. The numerical simulations using different values of the damage parameter n are shown in Fig. 5(a). In this figure, the destructuration process is controlled by the destructuration rate which depends on the damage parameter n . Based on this, different values of n can be determined for different amounts of cement.
Then, an undrained triaxial test on a normally consolidated sample is selected for the simulation, as shown in Figs

Test simulations
In this section, the model is used to simulate isotropic compression and undrained triaxial tests under various confining stresses on untreated and cement-treated Ariake clay and Singapore clay with different cement contents. All undrained triaxial tests were simulated by taking into account the first stage of isotropic compression before shearing as imposed on the specimen in the laboratory. The applicability of the model is investigated by comparing its predictions with experimental results.  Table 2).
1) Contact number per unit volume N/V and mean aggregate size d. Based on SEM photos, the mean size d of the clay aggregates is assumed to be 10 μm, and N/V is obtained by the mean size d and the void ratio e.
2) Inter-cluster elastic constants: k n0 , k rR and n. The exponent n = 1 is usually assumed in order to obtain a linear rebound line in the e-logp' plane. K rR = 0.5 can  usually be assumed for clay (see [7,8,[30][31][32][33]). The contact stiffness k n0 can be calculated by using the swelling index κ, measured from the swelling curve by using the relation k n0 ¼ 9ð1 þ e 0 Þ at =ð4κÞ. Therefore, only κ 0 and e 0 are needed as input parameters. For the Ariake clay, κ 0 = 0.086 and e 0 = 3.44 were measured.
3) Inter-cluster normal compression: σ pi0 and c p . The initial value of the intrinsic compression yield stress σ pi0 = 10 kPa is assumed for the untreated Ariake clay, since there is no measurement available. Using the intrinsic compression index l i and the swelling index κ of an assembly of clay aggregates, we can also determine the plastic inter-cluster compression index c p by the relation c p ¼ ðl i -κÞ=3ð1 þ e 0 Þ. Therefore, l i is used as an input parameter instead of c p . For the Ariake clay, l i = 0.446 was measured.
4) Inter-cluster friction angle: f μ and m. m = 1 is generally assumed for clay (see [7,8,[30][31][32][33]). f μ0 is assumed equal to the macroscopic friction angle of the normally consolidated reconstituted material. For Ariake clay, f μ0 = 40°was measured. 5) Inter-cluster plastic shear stiffness ratio k pR . According to the hardening rule used for shear sliding, a bigger value of k pR provides a stiffer shear stress-strain curve. For Ariake clay, k pR = 0.15 was obtained by fitting the deviatoric stress versus axial strain curve. 6) Dilatancy constant D. According to the local dilatancy relation, a bigger value of D gives bigger plastic volumetric changes, which results in higher excess pore pressure in undrained condition. Based on this, D = 6 was obtained for Ariake clay by fitting the effective stress path of the selected undrained test. 7) Critical state for the soil e cr0 . e cr0 is the reference void ratio at critical state corresponding to p ref ¼ 100 kPa. For Ariake clay, e cr0 = 1.911 was measured. 8) Cementation related parameters. Based on the analysis of the influence of cement content on the mechanical properties of cement-treated Ariake clay (see Section 2, Figs. 2 and 3), the cementation related parameters are summarized. 9) Inter-cluster debonding related parameters n and d . Based on two selected isotropic compression and undrained triaxial tests under a consolidated stress of 200 kPa on cement-treated Ariake clay (see Fig. 5b "c = 6%", Fig. 7(2) "cement content: 6%"), debonding related parameters can be determined (see Table 2). n controls the debonding rate of χ. For a given compression curve n = 10 was obtained (see discussion in Section 3.2). The parameter d = 3 was obtained by curve fitting of the selected undrained triaxial test (see discussion in Section 3.2).
All determined values of the parameters with the cement-dependent parameters (Figs. 2 and 3) were used to predict the other tests.

Singapore clay
The determination of parameters for cement-treated Singapore clay were based on the analysis of the influence of cement content on the mechanical properties of cementtreated Singapore clay (see the section 2, Figs. 2 and 3), considering one one-dimensional test (see Fig. 5(c) "c = 0%"), one isotropic compression test (see Fig. 5(c) "c = 10%") and one undrained triaxial test under a consolidated stress of 300 kPa on cement-treated Singapore clay (see Fig. 8-1 "cement content: 10%"). The procedure of parameter determination is similar to the one described for Ariake clay and is not repeated here. All determined parameters, summarized in Table 2, were used to predict all the other tests. Figure 5(b) shows the comparison between experimental results and numerical simulations of isotropic compression tests on untreated and cement-treated Ariake clay with different cement contents. Figure 5(c) shows the same comparison for one-dimensional compression test on untreated Singapore clay and isotropic compression tests on cement-treated Singapore clay with different cement contents. For both clays, a satisfactory agreement is achieved by using for each clay the corresponding set of parameters in Table 2. Figure 7 shows the comparisons for all the selected undrained triaxial tests on untreated and cement-treated Ariake clay with different cement contents by using for all the numerical simulations the same set of parameters presented in Table 2. Similarly, Fig. 8 shows the comparisons for all selected undrained triaxial tests on cement-treated Singapore clay with different cement contents. In the whole, the numerical simulations agree with the experimental results, although discrepancies were found for some of the tests, which may be due to some variations of the sample physical characteristics. The model appears to be able to reproduce the mechanical behavior of cement-treated clay and the effects due to the cementation and the damage process during shearing with reasonable accuracy. Figure 9 shows the comparisons between experiments and simulations for undrained triaxial tests on samples of Ariake clay with different cement contents consolidated up to 200 kPa. The influence of the cement content on the behavior of cement-treated Ariake clay was well reproduced by the proposed model: 1) The increase of stiffness and shear strength with cement content, accompanied by a Fig. 7 Comparisons between experiments and simulations for undrained shearing tests on cement treated Singapore clay with cement content of (1) 0%; (2) 6%; (3) 9%; (4) 12%; (5) 18%. (a) Effective stress path; (b) deviatoric stress versus axial strain; (c) excess pore pressure versus axial strain large post-peak reduction in the strength of the treated soil; 2) samples with low cement content (c = 0%, 6%, 9%, 12%) exhibit a contractive behavior and samples with high cement content (c = 18%) a dilative one. Figure 10 shows the comparisons between experiments and simulations for undrained tests on Singapore clay with different cement contents consolidated up to 500 kPa. Again, the influence of cement content was well reproduced by the proposed model.

Micromechanical analysis of induced anisotropy
In this section, we investigate the local stress-strain behavior for selected contact planes with different orientations θ = 0°, 18°, 28°, 45°, 55°, 72°and 90°(see Fig. 4, x representing vertical direction). The undrained triaxial tests on samples of Ariake clay with cement contents of 0%, 9% and 18% under a consolidation stress of 200 kPa (see Fig. 9) were chosen as an example for illustrating the local behavior in contact planes (see Figs. 11 and 12). The influence of the cement content on the local stress-strain behavior is thus discussed. To study the stress-induced anisotropy during undrained shearing, we selected for all the tests three steps for plotting the variable distribution in the rose diagrams (see Fig. 13): the beginning of loading, the peak deviatoric stress and the end of loading.  Figure 11(a) presents the local stress paths for all selected contact orientations (from 11-1 to 11-5) for all cases (cement contents of 0%, 9% and 18%). For the 0°and 90°c ontact planes, the shear stress is null. Thus, the stressstrain behaviors are not plotted. For the others, the stress state with the highest stress ratio is located in the 70°c ontact plane. Figure 11(b) shows the local shear stressstrain curves which clearly indicate that every contact plane is mobilized to a different degree. The planes with the largest movements are located around the orientation of 72°(close to π/4 + f μ /2 = 65°). These active contact planes contribute largely to the overall deformation of the specimen. Among all selected contact planes, only the 72°contact plane has a local shear stress-strain curve similar to the global one. Figure 11(c) presents the local normal stress-strain curves, in which normal strains are mostly composed by elastic strains due to normal stresses and plastic strains induced by the dilatancy law during shearing. The results show that, although the global volumetric strain is null during undrained shearing, local normal strains develop in every contact plane to different  For each sample with different cement content, each contact plane behaves differently, and the influence of the cement content on the local stress-strain behavior is significant, which agrees with the global behavior (see Fig. 9).
The degradation of the bonding ratio and of the bonding cohesion for all selected contact planes are plotted versus the global axial strain in Fig. 12 for samples with cement contents of 9 and 18%. For the bonding cohesion (Figs. 12 (a)-(b), most of the degradation took place in the 72°c ontact plane due to its largest local displacements. The other contact planes have only a slight debonding because their mobilized displacements are small. For the bonding ratio ( Fig. 12(c)), a very small reduction took place because the maximum normal strain remained smaller than 1% (see Figs. 11(c)).

Local stresses, strains and bond distributions
The distributions of the local normal stiffness and bond in contact planes of various orientations are plotted in Fig. 13  (Fig. 13-1 for untreated clay, Fig. 13-2 for cement content of 9% and Fig. 13-3 for cement content of 18%).
1) The normal stiffness is a function of the normal stress and of the amount of cement. The normal stiffness distribution at the end of the isotropic consolidation has a circular shape (see the bold line in Figs. 13(1-a), (2-a) and (3-a) which corresponds to an isotropic distribution of the normal stress for all contact planes. During the undrained shearing, different cement contents give different evolutions of the stiffness distribution: for 0% cement content, the distribution becomes anisotropic with an elliptical shape having the long axis in the vertical direction; for 9% cement content, the distribution becomes also anisotropic with higher stiffness in the vertical direction and lower stiffness in the horizontal direction, and then shrinks during shearing due to bond damage; the same evolution is found for 18% cement content with higher stiffness values due to stronger bonds.
2) During the undrained shearing, the shear stress distribution expands while maintaining a similar shape for all cases (Figs. 13(1-b), (2-b) and (3-b)). This distribution shows a concentration of the maximum local shear stresses between the 45º and 55º orientations, in agreement with the local stress-strain curves in Fig. 11.
3) During the undrained shearing, the shear strain distribution for all cases (Figs. 13(1-c), (2-c) and (3-c)) shows that very large strains occurred at the end of the loading in contact planes around the 70º orientation, which agrees with the local stress-strain curves in Fig. 11.
4) The distribution of shear bonding cohesion changes slightly from step 1 to 2 due to small mobilized shear plastic strains for both 9% and 18% cement contents. Then, from step 2 to 3, the bonding cohesion reduces significantly in the 55º-70º contact planes, but does not change very much in the other planes, which agrees with the evolution plotted in Figs. 11(a) and (b). This also agrees with the shear strain distribution in Figs. 13(2-c) and 13(3-c), because the damage of the shear bonding cohesion is controlled by the amount of plastic shear strain.
In the present model, stresses and bonds in each plane are considered as internal state variables and their different evolution in each individual plane introduces in a natural manner the stress-induced anisotropy. Since many soil properties are stress-dependent (e.g., the plastic hardening parameter k p and the shear modulus k r are related to k n which is stress-dependent), the induced anisotropic behavior during undrained shearing is linked to the stress distribution in each plane. The analyses of local stressstrain behaviors also indicate that for real sample, the debonding process is more located in the direction of shear band.

Conclusions
Cementation among clay clusters is formed when a certain amount of cement is mixed with saturated clay. The influence of the cement content on the mechanical properties of cement-treated clays was investigated. The debonding effects on the stress-strain behavior of cementtreated clay were also discussed. Based on experimental results, a micromechanical stress-strain model was developed considering explicitly the cementation at inter-cluster contacts. The inter-cluster bonding is introduced by considering two aspects: an additional cohesion in shear sliding and a higher yield stress in normal compression.
Isotropic compression and undrained triaxial tests under various confining stresses on cement-treated Ariake clay and Singapore clay with different cement contents were The predicted behavior in contact planes was examined in the case of undrained triaxial tests with different cement contents. The local stress-strain responses in contact planes showed that every contact plane was mobilized to a different degree. A few active contact planes contributed mainly to the deformation of the assembly, while most contact planes contributed only slightly to the specimen straining. Therefore, the local strains are highly nonuniform. It has been shown from the rose diagrams that the shape of the contact stresses, strains and bond distributions, controlling the evolution of local soil properties, changes throughout the mechanical loading, which clearly indicates the development of anisotropy induced by the externally applied loads.
where p' is the mean effective stress of the packing; l is the slope of critical state line in elogp' plane, e cr0 is the reference void ratio at critical state corresponding to p ref = 100 kPa. The resistance can be related to the state of packing void ratio e by where m is a material constant (m = 1 can be generally assumed).
3) Overall stress-strain relationship During the integration process, the static hypotheses are used to obtain the relation between the strain of assembly and inter-cluster strain where _ γ j is the local strain between two contact clusters; n k is the unit vector of the branch joining the centers of two contact clusters, and N is the total number of contacts, over which the summation is carried out. The tensor B α ik in Eq. (A12) is defined as B α ik ¼ A -1 ik ðl α Þ 2 , where the f abric tensor Using the principle of energy balance, which states the work done in a representative volume element equal to the work done on all inter-cluster planes within the element and using Eq. (A12), the local stress on the αth contact plane is derived as follows: The stress increment _ ij can be obtained by the contact forces and branch vectors for all contacts. In terms of local stress, it is Using Eqs. (A9), (A12), and (A15), the following relationship between stress and strain can be obtained: A numerical calculation with an iterative process is used to carry out the summation in Eq. (A18).