#### Abstract

Gravity is the primary load source in geotechnical engineering, and its effect becomes extremely important in these cases such as lunar research station construction and in situ resource utilization. However, at the current stage, the study on the effect of gravity-induced stress gradients on the behavior of geomaterials is seriously lacking. In this work, the size-dependent behaviour of representative volume element (RVE) of granular matter under the gravity-induced stress gradient is studied. A homogenization theory is proposed firstly, which homogenizes the granular matter into stress gradient continuum through RVE. Then, numerical triaxial loading for RVEs with various heights and particle sizes is carried out under 0∼10 g gravity. The results show that the size effect and its origination appear to have no obvious influence on isotropic compression. However, the macroscopic shear strength, elastic modulus, and volumetric deformation all show size effect dependence, and the gravity-induced stress gradient leads to significant differences between different originations of size effect. The valence, orientation, and elongation degree of microscopic void cells are further analyzed, which suggests that applying of stress gradient results in the evolution of void cells demonstrates a dependence on both the size effect and its origination. The results indicate that the gravity-induced stress gradient and corresponding size effect should be paid sufficient attention to granular mechanics of terrestrial and off-earth engineering.

#### 1. Introduction

Gravity is the primary load source in geotechnical engineering, which results in a uniform body force inside geomaterials. The effect of gravity on the behaviour of geomaterials becomes extremely important in cases suffering from terrestrial and nonterrestrial gravities. One typical situation is the scaled physical modelling, in which significantly larger gravities were applied (e.g., [1–3]). On the other hand, with limited resource reserves and increased environmental threat on the Earth, off-earth base construction and in situ resource utilization (ISRU) have become urgent in recent decades [4–6]. The primary challenge in off-earth geotechnical engineering is the extraterrestrial gravities which are quite different from that on the Earth, for instance, 1.662 N/kg on the Moon [7] and 3.711 N/kg on the Mars [8]. A central premise of scaled physical modelling is that the constitutive behaviour of geomaterials depended only on the gravity-induced stress level (i.e., confining pressure) and not on the gravity-induced stress gradient. Workablility of this premise would also guarantee the current geomechanics being applicable on other planets and asteroids. However, studies performed under nonterrestrial gravities have showed that the constitutive behaviour of geomaterials is related to the stress gradient [9–12].

Most current geomechanical theories were developed under the framework of classical continuum mechanics, building on a concept of material point (e.g., [13–15]). On the one hand, employment of the material point results in higher order stress and strain (i.e., stress gradient and strain gradient) being theoretically excluded from the constitutive behaviour of geomaterials. On the other hand, as typical heterogeneous media, classical continuum mechanics are only applicable for geomaterials when the macroscopic wavelengths of relevant field quantities are significantly larger than the characteristic microstructural dimensions. Therefore, the homogenization theory was proposed defining the macrofield quantities as the volume average of the corresponding microscopic counterparts over a representative volume element (RVE) [16, 17]. The determined RVE should be (a) structurally entirely typical of the overall particle assembly; (b) sufficiently large so that the averaged macroscopic field quantities are uniform in the vicinity of it and significantly vary only in a much larger macroscopic scale [17–19]. The intrinsic size of particles and the minimum size of RVE together determine whether the RVE is “representative” when homogenizing the heterogeneous geomaterials to Cauchy continuum [20–22]. Since the higher order stress and strain were ignored, both the particle size independent and RVE size independent constitutive behaviour can be obtained when the RVE was sufficiently “representative.”

However, a limitation of classical geomechanics was manifested already when describing the cases with higher order stress or strain. Several studies have devoted to take the higher order field quantities into account by considering the concept of generalized continuum (e.g., [23–25]). The generalized continuum-based higher order/grade continuum mechanics was generally formulated by adding the higher order items into additional balance equation. Thereby, it further requires to formulate respective additional constitutive relation to account for the higher order field quantities [24–27]. The additional constitutive relation can be obtained by homogenizing the heterogeneous media to generalized continuum over a RVE [28–30]. Therefore, an issue arise that the macroscopic constitutive relation may depend on the size of RVE because of the higher order field quantities being included. Hütter [30] pointed out that the RVE size dependency of the constitutive relation is inevitably necessary since the potential generalized boundary condition on higher order field quantities require an interpretation with respect to the location of the macroscopic boundary relative to material heterogeneity. Li et al. [29, 31] proposed a homogenization methodology to transfer the RVE size dependency of the constitutive relation to intrinsic size dependency. The dependency of macroscopic constitutive behaviour for heterogeneous media on the intrinsic particle size and RVE size is an important issue under the framework of higher order/grade continuum mechanics.

The concept of generalized continuum provides an effective way to account for the effect of gravity-induced stress gradient on the constitutive behaviour of granular geomaterials. As mentioned above, both influences of the stress gradient and the size effect become important due to the adoption of RVE. In the previous studies of this work [11, 12], we have revealed the effect of gravity-induced stress gradient on the behaviour of RVE at both macroscopic and microscopic scales. Furthermore, the present contribution will study the dependency of the constitutive behaviour on the particle size and the RVE size under the gravity-induced stress gradient condition. The present work is organized as follows: Section 2 first derived the homogenization theory for the gravity-induced stress gradient continuum, based on which a series of numerical triaxial loading for the RVE were designed with various size and stress gradient in Sections 3–5 studied the effect of the particles size and the RVE size on the constitutive behaviour of granular materials under the uniform and gradient stresses from macro- and microperspectives, respectively; the conclusion remarks are finally given in Section 6.

#### 2. Homogenization from Geomaterials to Stress Gradient Continuum

##### 2.1. Gravity-Induced Stress Gradient Tensor

The stress gradient *R* in continuum mechanics belongs to the third rank tensor and defined asin which *σ* is the Cauchy stress tensor, ⊗ is the tensor product, *σ*_{ij,k} is components of *R*, *e*_{i} is the Cartesian orthonormal basis, and ▽ is the nabla operator. Since the components of stress gradient tensor are symmetric with respect to the first two indices, there are only 18 independent components [32].

When a gravity is applied, the balance equation of geomaterials is given bywhere *f*_{i} denotes the gravity-related body force and is given as follows:in which is the component of gravity acceleration in the direction *i* and *ρ* is density.

Meanwhile, relationships between the gravity-induced stress gradients of normal stresses and the components of gravity in respective normal directions are given by

A combination of equations (2)–(4) results in

Therefore, the number of independent components of the gravity-induced stress gradient tensor reduced to 15.

Furthermore, the gradients of shear stresses vanished when the gravity-induced stress gradient tensor was transferred to the space of principal stresses. Therefore, the number of independent components further reduced to nine. The nonzero components of gravity-induced stress gradient tensor in the space of principal stresses can be expressed as follows with matrix form:

If the consolidation of geomaterials under *k*_{0} condition was taken into account, the components in equation (6) can be calculated bywhen a predetermined gravity was applied, where *k*_{i} is the lateral pressure coefficients of geomaterials.

In a more simple case, when ignoring the *k*_{0} consolidation, only 3 nonzero components of the stress gradient tensor are remained, as follows:

##### 2.2. Link between Macroscopic and Microscopic Field Quantities

The RVE was employed to adapt the macrohomogeneity concept to the cases when considering the gravity-induced stress gradient. Assume the volume of the RVE is *V*, and the exterior boundary is Ω. The global coordinate system is *x*_{i}, while a local coordinate system *X*_{i} is adapted associated with the RVE and originated at its geometrical center.

Assume a macroscopic stress field *σ*_{ij} of the RVE is issued from smoothing the oscillating microscopic stress field . The equivalence between the macroscopic and microscopic stress field is guaranteed when their average gradients over a RVE under order *n* keep the equality, namely,where the average of a field quantity *ζ* over a RVE is defined by

Expanding the macroscopic stress field into a Taylor’s series at the geometrical center *X*_{c} of the RVE results inin which *O*(*x*^{4}) is the higher order remainder terms.

Therefore, the average of the macroscopic stress field over the RVE can be expressed as follows by using equation (11):in which , , and are given by

From the similar algebraic, the average of macroscopic stress gradient field over the RVE can be expressed as follows:

Consequently, by taking equations (9), (12), and (1) and (14) into account, equation (11) can be rewritten as

Thus, the link between the macroscopic stress field and the average of the microscopic stress field and its gradient field is established.

Furthermore, the traction *t*_{i} on the boundary of the RVE can be expressed as

On the other hand, a so-called microdisplacement degree Φ [32] was introduced to the stress gradient continuum in addition to the usual displacement components. The microdisplacement degree is a third rank tensor which is conjugate to the stress gradient tensor R. Introducing microdisplacement degrees results in a generalized stain measure which is conjugate to the stress tensor, namely,

From similar algebraic for equation (17), the macroscopic displacement field and macroscopic microdisplacement field can be expressed as follows, respectively, to establish a link to the average of corresponding microscopic quantities:

It is noted that the rigid body movement is omitted from equation (20).

##### 2.3. Macrohomogeneity Condition for RVE

Consider a RVE loaded with traction *t*_{i} on its exterior boundary Ω, and assume a virtual kinematically admissible displacement field occurred on Ω:

Since this work concentrates on the gravity-induced stress gradient, only the first order gradient of the stress, displacement, and microdisplacement was taken into account, namely, equations (17)–(21) will be truncated after the first-order terms in the following derivation.

The virtual work of the traction *t*_{i} on the virtual displacement field *δħ*_{i} on the boundary of the RVE can be expressed as

The virtual work in equation (23) can further be rewritten as follows by taking equations (20) and (23) into account:

By combining equation (18), the first term on the right hand of equation (24) can further be written as

It is noted that equation (24) further defined the average stress measure of the RVE with respect to the boundary traction, which is expressed as

The second term on the right hand of equation (24) can further be written as follows by combining equation (18):

Similarly, equation (27) further defined the average stress gradient measure of the RVE with respect to the boundary traction, which is expressed as

The third term on the right hand of equation (24) can further be written as follows by combining equation (18):

Consequently, by combining equations (25), (27), and (29), equation (24) becomes

The above equation further defines the average strain measure of the RVE with respect to the average displacement and average microdisplacement, namely,

It is notable that the right hand of equation (30) is actually the variation of the average strain energy density of the RVE, as expressed by

Therefore, equation (30) can further be rewritten aswhich can read as the virtual work principle applied on the RVE. Equation (33) gives the macrohomogeneity condition for the RVE, also known as the Hill–Mandel energy condition [16, 33].

On the other hand, the constitutive equation for the RVE can further be given according to equation (32), as expressed by

It is important to point that this constitutive equation is obtained from homogenization of the RVE and stands for the constitutive behaviour of the RVE. It can further be transferred to the constitutive behaviour of the stress gradient continuum by homogenization, which is beyond the scope of this work and will not be discussed here.

##### 2.4. Boundary Condition for RVE

The boundary condition for the RVE is essential to whether it is ‘representative’ when obtaining its constitutive behaviour. According to equation (18), the traction boundary condition for the rectangular RVE can be expressed aswhen only considering the first-order gradient of the stress. Here, and stand for the boundary surface on the positive and negative direction of the local coordinate system, respectively, and *l*_{i} stands for the length of the RVE in *x*_{i} direction.

If the consolidation of geomaterials under *k*_{0} condition was taken into account, the traction boundary condition for the RVE in the space of principal stresses can be simplified as follows by taking equation (7) into account:

Furthermore, a more simple traction boundary condition for the RVE in the space of principal stresses can be obtained when ignoring the *k*_{0} consolidation, namely,

#### 3. Numerical Experiment Schedule

##### 3.1. Model Setup

The size-dependent behaviour was studied by a series of numerical triaxial tests for a rectangular 2D RVE using the discrete element method (DEM). A simple case was considered in which the gravity was applied with its direction coinciding with the major principal stress. Meanwhile, the effect of *k*_{0} consolidation was ignored for the purpose of concentrating on the influence of a single stress gradient. According to these preconditions and equation (37), a physical model can be built for the rectangular RVE with its boundary condition showing in Figure 1. *σ*_{1} and *σ*_{3} denote the major and minor principal stress, and *ρ*, , and *l*_{1} are the bulk density, gravity, and RVE height, respectively.

The numerical triaxial simulations for the RVE were performed using the DEM code *PFC*^{2D}, which is based on the application of Newton’s conservation of momentum law to entities by employing an explicit time-stepping algorithm on the central-difference scheme [34]. The linear contact model—slip model—together with the viscous damping model [35] was adopted to regulate the behaviour of particles and contacts. The initial RVE was prepared with 2D rigid circular particles. Four rigid boundaries were employed to perform loading, and only the translation of them was allowed. Specific physical parameters are given in Table 1.

As shown in Figure 2, the numerical triaxial testing included four stages:(1)*Preparing Initial RVE* (Figure 2(a)). The initial RVE was prepared by free falling several layers of particles with equivalent height. A gravity of 100 g was applied at first to get a dense assembly. The particles in each layer were fixed immediately after completing deposition. All particles were set free after completing the deposition of the whole RVE and then permitted to achieve an equilibrium state. A top boundary was generated finally at the target height.(2)*Unloading Gravity* (Figure 2(b)). The gravity was altered from 100 g to the predetermined intensity. To minimize the possible disturbance on the microstructure of the RVE, the gravity was reduced in 1/1000 g steps, immediately followed by a quasistatic equilibration, to eliminate imbalance.(3)*Loading Confining Pressures* (Figure 2(c)). The predetermined isotropic confining pressures were applied to the RVE by fixing the bottom boundary while performing a servo-control mechanism on the others. This servo-control mechanism achieved a stress state on RVE as shown in Figure 2(c).(4)*Triaxial Loading* (Figure 2(d)). A velocity of 1.0 × 10^{−4} m/s was applied to the top boundary to perform triaxial loading. The bottom boundary remained fixed, while the lateral boundaries were moved by the servo-control mechanism to maintain a constant confining pressure.

##### 3.2. Testing Schedule

The ratio of the height to width of the RVE maintained a constant being 2.05, while the radius of particles distributed linearly and maintained a constant radius ratio (defined as the ratio of maximum to minimum radius) being 2.0. Therefore, the size of the RVE can be fully described by its height *l*_{1}, while the size of the particles can be fully described by their medium radius *R*_{m}. A size coefficient of RVE, *λ*, was defined by equation (38) in this study.

The triaxial loading schedule for the RVE is given in Table 2. Two groups of tests were designed with a same variation in size coefficient from 40.92 to 286.45, but with different height of RVE and particle medium radius. The variation of the size coefficient in subgroups numbered by “*λ*-Hn” was caused by the increase in height of the RVE, while that in other subgroups numbered by “*λ*-Rn” was caused by the decrease in particle radius. Each subgroup of RVEs were tested under the gravities of 0 g, 1 g, 2 g, 5 g, and 10 g, respectively. Furthermore, the RVEs under each gravity were tested under the confining pressures of 10 kPa, 20 kPa, 30 kPa, and 40 kPa, respectively.

A target porosity of 0.20 was predetermined for all initial RVEs for each group of test. In order to achieve that the preparation of the initial RVE in each group was repeated several times until an initial RVE with its porosity less than 0.5% of the target porosity is obtained. Statistics in porosity of all initial RVEs suggested an initial porosity of 0.20 ± 0.11% finally.

#### 4. Size-Dependent Macroscopic Behaviour of RVE

##### 4.1. Isobaric Consolidation

Figure 3 gives the porosity of RVEs with various heights and particle sizes after applying 20 kPa confining pressure under conditions 0 g and 5 g, respectively. The results show that the size effect has no obvious influence on the variation of porosity during isobaric consolidation. The porosity of RVEs with different sizes, caused by the variation of both the RVE size and particle size, is reduced after applying 20 kPa consolidation pressure. At the same gravity condition, the decrement of porosity keeps nearly unchanged along with the increase of *λ*. On one hand, the 2D particle assembly is already in a dense state at the porosity of 0.20, so the decrement of porosity during isotropic consolidation is relatively small. On the other hand, since the microscopic structure of assembly is not affected by the size effect related to RVE size and particle size, the consolidation is not affected either.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 3 demonstrates that gravity has a slight influence on the consolidation of RVEs of the same size. For the RVEs with various heights, an average reduction of 0.10% in porosity is observed at 0 g gravity, while that is 0.08% at 5 g gravity. The corresponding percentages for the RVEs with different particle sizes is 0.14% and 0.08%, respectively. This suggest that when the RVEs with the same size are isotopically compressed, the reduction in porosity at the 0 g condition is slightly larger than that at 5 g gravity. This suggests that the compression of RVE under uniform stress is more significant than that at gradient stress.

##### 4.2. Strength Characteristics

Figures 4(a) and 4(b) give the relationships between deviatoric stress and axial strain of RVEs with different size coefficients under the 0 g gravity and 20 kPa confining pressure. Both the size effects caused by the variation of RVE height and particle size are demonstrated. At the smaller *λ* cases in Figure 4(a), the deviatoric stress shows a significant local fluctuation along with the increasing axial strain. This fluctuating tendency reduces gradually with the increase of *λ*. For the size variation caused by RVE height, the deviatoric stress increases at first followed by a stable state with the increasing *λ*; the corresponding critical *λ* is 204.604.

**(a)**

**(b)**

**(c)**

**(d)**

For RVEs with different size coefficients caused by particle size in Figure 4(b), the evolution of deviatoric stress with respect to *λ* is similar to that of RVEs at various heights. A comparison between Figures 4(a) and 4(b) suggest that, despite the different origins of size effect, the evolution of deviatoric stresses and their peak values of RVEs with equivalent *λ* is basically consistent under the 0g gravity. This indicates that the RVEs with unique *λ* own similar deviatoric stress-axial strain relationships under uniform stress conditions.

Figures 4(c) and 4(d) give the relationship between deviatoric stress and axial strain of RVE at 5 g gravity and 20 kPa confining pressure, in which both the size effects related to RVE height and particle size are demonstrated. When considering the size effect caused by RVE height in Figure 4(c), the deviatoric stress-axial strain relationship shows a similar tendency to that under the 0g gravity. There is larger deviatoric stress with significant local fluctuation when *λ* is smaller than 122.762, but both the peak deviatoric stress and its fluctuation reduce gradually with the increasing *λ* and become stable finally when *λ* exceeds 204.604. A comparison with Figure 4(a) shows that, under 5g gravity, the deviatoric stress of RVE is less than that of RVE with the same *λ* at 0g gravity. This indicates that the deviatoric stress of RVE decreases along with increasing gravity-induced stress gradient.

For the group with various particle sizes Figure 4(d), the deviatoric stress - axial strain relationship shows a similar tendency to the RVEs with different heights. However, a comparison with Figure 4(c) suggests that the deviatoric stress in the group with various particle sizes is generally smaller than that of RVE with the same *λ* in the group with various RVE heights. By further comparison with Figures 4(a) and 4(b), it can be found that the decrement of deviatoric stress from 0g to 5g in the former group is greater than that in the latter group. The significant reduction in deviatoric stress occurs within the RVE whose *λ* increases from 40.921 to 122.76. When considering the absolute particle size in RVE with *λ* of 40.921∼122.76, as given in Table 2, it can be found that the diameter in the former group is slightly larger than that in the latter group. Therefore, when a gravity-induced stress gradient is applied, the larger particle size has a more significant influence on the deviatoric stress.

Figures 5(a) and 5(b) give the relationship between the principal stress ratio and axial strain of RVEs with *λ* of 40.921 and 286.445 at 0 g gravity, in which both the RVE height and particle size which caused size effects are demonstrated. For the group with various heights in Figure 5(a), when the confining pressure increases from 10 kPa to 30 kPa, the principal stress ratios of RVEs with the same *λ* show good consistency under small axial strain. With further increase of strain, relatively larger differences occur in the RVEs with *λ* of 40.921; but the basic synchronization still exists. These RVEs with *λ* of 286.445, however, keep good consistency in the whole process. Both the magnitude and local fluctuation of the principal stress ratio decrease significantly along with the increasing *λ*, and the decrements are basically the same at different confining pressures.

**(a)**

**(b)**

**(c)**

**(d)**

For RVEs with various particle sizes in Figure 5(b), the evolution of the principal stress ratio is similar to that of the group with different RVE heights. A comparison with Figure 5(a) shows that, although the origin of the size effect is different, the principal stress ratios are basically the same at 0 g gravity when *λ* of RVE keeps consistent.

Figures 5(c) and 5(d) give the relationship between the principal stress ratio and axial strain of RVEs with a *λ* of 40.921 and 286.445 at 5 g gravity, in which the influence of RVE height and particle size are demonstrated. When considering the size effect related to RVE height in Figure 5(c), the evolution of the principal stress ratio at 5 g gravity is similar to that at 0 g. By further comparison with Figure 5(a), it can be found that when the gravity increases from 0 g to 5 g, the principal stress ratio shows a decreasing tendency under all *λ* and confining pressures. This suggests that the principal stress ratio, regardless of RVE size and confining pressure level, decreases along with increasing gravity-induced stress gradient.

For RVEs with various particles in Figure 5(d), when the confining pressure increases from 10 kPa to 30 kPa at 5 g gravity, the relationship between the principal stress ratio and the axial strain shows a similar tendency to the RVEs with different heights. However, the principal stress ratio of RVEs in Figure 5(d) is slightly smaller than that of RVEs with the same *λ* in Figure 5(c). This indicates that the attenuation of the principal stress ratio caused by the decreasing particle size is greater than that caused by the increasing RVE height under the gravity-induced gradient stress condition. A further comparison with Figure 5(b) shows that when the gravity increases from 0 g to 5 g, the principal stress ratio of RVE with the same *λ* demonstrates a slightly decreasing tendency.

Figure 6 gives the relationship between peak principal stress ratio and size coefficient of RVEs at 0 g and 5 g gravities, respectively. The peak principal stress ratio in figures is the average of the principal stress ratio under confining pressures of 20∼40 kPa, and the corresponding error bars are also given. It can be found from Figure 6(a) that the peak principal stress ratio decreases at first followed by a stable state with the increase of *λ* at 0 g gravity. The critical value of *λ* relative to the stable peak principal stress ratio is about 204.604. Meanwhile, the origins of the size effect have no obvious influence on the evolution of the peak principal stress ratio.

**(a)**

**(b)**

For the RVEs at 5 g gravity, the peak principal stress ratio shows a similar tendency to that at 0 g, but a smaller critical *λ* is observed to be 160. At the same time, the peak principal stress ratio of the RVEs in the group with various particle sizes is significantly smaller than that in the group with various RVE heights. When the size coefficient *λ* increases from 40.921 to 286.445, the difference of peak principal stress ratio between the groups with different origins of size effect reduces gradually. This indicates that the origins of the size effect have a relatively large influence on the peak principal stress ratio under the gravity-induced gradient stress. Meanwhile, this influence is strengthened by the increasing particle size. A comparison with Figures 6(a) and 6(b) suggests that when the gravity increases from 0 g to 5 g, the peak principal stress ratio shows a decreasing tendency in the whole range of *λ*.

##### 4.3. Internal Friction Angle

Figure 7 gives the relationship between the internal friction angle and size coefficient of RVE under the gravity of 0 g and 5 g and also shows the size effects caused by RVE heights and particle sizes. For the case of 0 g in Figure 7(a), with the increase of *λ*, the internal friction angle decreases at first followed by a stable state. The critical value of *λ* relative to the stable internal friction angle is about 204.604. Regardless of the origins of size effects, the internal friction angle of RVEs with the same *λ* at 0 g gravity are basically equal. On the other hand, the error variation of the internal friction angle is obviously larger at small *λ*. This error variation is essentially caused by the local fluctuation of microscopic interparticle force and therefore is also the basic property of RVE. When the size of RVE is small, the macroscopic average of stress and other quantities will be affected by the local fluctuation of microscopic quantities. However, this local disturbance is gradually dismissed by the increasing RVE size.

**(a)**

**(b)**

For the RVEs under 5 g gravity in Figure 7(b), the internal friction angle shows a similar tendency to that at 0 g gravity. However, the reduction of internal friction angle in the group with various particle sizes is smaller than that in the group with various RVE heights. The difference in internal friction angle between these two groups is significant when *λ* is smaller than 160, and it decreases gradually along with the increasing *λ*. The above observations suggest that the size effect has an obvious influence on the internal friction angle at the larger particle size cases. A comparison with Figure 7(a) indicates that when the gravity increases from 0 g to 5 g, the internal friction angle decreases slightly, regardless of the size effects and their origins.

##### 4.4. Elastic Modulus

Figure 8 gives the relationship between the elastic modulus and size coefficient of RVEs under the gravity of 0 g and 5 g and also shows the influence of size effects caused by RVE height and particle size, respectively. In the case of 0 g in Figure 8(a), with the increase of *λ*, the elastic modulus firstly decreases and then tends to stabilize. However, the elastic modulus of RVEs in the group with various heights shows a slightly increasing tendency after *λ* is larger than 204.604. By combining Figure 3, one can easily find out that this abnormal phenomenon is caused by the fluctuation of initial porosity, which has nothing to do with the size effect and thus will not be discussed further. The elastic modulus of RVEs with different size coefficients and their origins under 0 g gravity are quite close in magnitude when *λ* is small than 160.

**(a)**

**(b)**

As shown in Figure 8(b), regardless of the origins of the size effect, the elastic modulus of RVEs under the gravity of 5 g decreases first and tends to be stable along with increasing *λ*. However, the elastic modulus of RVEs in the group with various particle sizes is smaller than that in the group with various RVE heights. Meanwhile, the reduction of elastic modulus in the former group is obviously less than that in the latter group. The difference in elastic modulus between the two groups decreases along with the increasing *λ* and becomes insignificant when *λ* exceeds 204.604. These observations suggest that when the gravity-induced stress gradient is applied, the size effect has a larger influence on the elastic modulus at the larger particle sizes. A comparison with Figure 8(a) shows that, regardless of the origins of the size effect, the elastic modulus at 5 g gravity is statistically smaller than that at 0 g. This indicates that the elastic modulus decreases with the increase of gravity.

##### 4.5. Volumetric Strain

Figures 9(a) and 9(b) give the relationship between the volumetric strain and axial strain of RVEs with various heights and particle sizes under the gravity of 0 g, and all these RVEs are tested under the confining pressure of 20 kPa. For the RVEs with various heights in Figure 9(a), there is significant volumetric dilatation at 0 g gravity and small *λ*. Meanwhile, the volumetric strain shows an obvious local fluctuation at the large axial strain. With the increase of *λ*, the volumetric dilatation gradually tends to be volumetric compression and finally becomes stable when *λ* exceeds 204.604. At the same time, the local fluctuation of volumetric strain is also eliminated by the increasing *λ*.

**(a)**

**(b)**

**(c)**

**(d)**

When considering the size effect caused by particle size in Figure 9(b), the RVE with *λ* of 40.921 also shows obvious volumetric dilatation at 0 g gravity, and a local fluctuation at large axial strain is also observed. A comparison with Figure 9(a) shows that the volumetric dilatation of RVE in the group with various particle sizes is slightly greater than that in the group with various RVE heights. The volumetric dilatation also tends to be volumetric compression and finally stabilizes at the critical *λ* of 204.604.

Figures 9(c) and 9(d) give the relationship between the volumetric strain and axial strain of RVEs under the gravity and confining pressure of 5 g and 20 kPa, respectively. When considering the size effect caused by the RVE height in Figure 9(c), it can be found that the RVE with *λ* of 40.921 shows a slight volumetric dilatation. Under the large axial strain, the local fluctuation of the volume strain is also observed. With the increase of *λ*, the volume of RVE tends to be compressed and stabilizes finally when *λ* exceeds 204.604. A comparison with Figure 9(a) shows that when the gravity increases from 0 g to 5 g, both the volumetric change and its local fluctuation are decreased significantly. The above observations suggest that both the volumetric dilatation and compression of RVE are weakened when applying a gravity-induced stress gradient.

For the RVEs in the group with various particle sizes, as shown in Figure 9(d), the significant volumetric dilatation is observed under 5 g gravity when *λ* is equal to 40.921. Different from other cases, the volumetric strain does not show any obvious local fluctuation at large axial strain. With the increase of *λ*, the volumetric dilatation tends to become compression and finally stabilizes when *λ* is larger than 204.604. A comparison with Figure 9(c) indicates that the volumetric dilatation in the group with various particle sizes is greater than that in the group with various RVE heights at small *λ*. However, with the increase of *λ*, the influence of the origin of size effect becomes less obvious. A further comparison with Figure 9(a) shows that when the gravity increases from 0g to 5 g, both the volumetric dilatation and its local fluctuation are dismissed. At the same time, the volumetric compression also seems to be weakened by the gravity-induced stress gradient when *λ* exceeds 122.762.

#### 5. Size-Dependent Microscopic Behaviour of RVE

To further analyze the size dependence of RVE under gravity-induced stress gradient, the evolution of microscopic behaviour is studied. The void-based microscopic representation [36] is employed here (as shown in Figure 10), which describes the RVE by a particle graph that consists of void cells. The void cell is the minimum element that can define both microscopic and macroscopic quantities, and thus, its valence, orientation, and elongation degree are important microscopic properties of RVE. The average valence of void cells is defined bywhere is the side number of the void cell and is the number of void cells. The orientation angle of void cell is defined bywhere the fabric tensor is given byin which is the set of interparticle contacts consisting of void cell and and are the contact normal vectors. The elongation degree in the direction of the gravity is defined by

In this section, the evolution of the average valence, orientation, and elongation degree of void cells is analyzed.

##### 5.1. Valence

The valence is defined as the average side number of void cells. Figure 11gives the evolution of valence in RVEs with various heights and particle sizes under the gravity of 0 g and 5 g, respectively. The initial valence in Figure 11(a) suggest that, with the increase of *λ*, the valence increases at first followed by a stable state during isotropic consolidation at 0 g gravity. For RVEs with various particle sizes in Figure 11(b), the initial valence shows a continuous increasing tendency along with increasing *λ*, but its increment is reducing. When consolidated at the small *λ*, the valence in the group with various particle sizes is smaller than that in the group with various RVE heights. However, with the increase of *λ*, the valence in the former group increases faster, which results in a larger initial valence than the latter group finally.

**(a)**

**(b)**

**(c)**

**(d)**

During the isotropic consolidation under 5 g gravity, as shown in Figures 11(c) and 11(d), a similar evolution of valence is observed. A comparison of Figures 11(c) and 11(a) suggests that the initial valence at 5 g gravity is smaller than that at 0 g. Comparison of Figures 11(d) and 11(c) shows that the initial valence of RVEs in the group with various particle sizes is larger than that in the group with various RVE heights at small *λ*. However, with the increase of *λ*, the final valence tends to be very close. This indicates that gravity has a greater influence on the valence of RVE with the larger particle size. A comparison between Figures 11(d) and 11(a) shows that the initial valence increases along with the increasing gravity at small *λ*, but the opposite trend is observed at large *λ*. This indicates that the particle size has a more significant influence on the gravity-related evolution of valence at small *λ*.

When considering the RVEs with various heights under 0 g gravity, as shown in Figure 11(a), the valence of RVE with small *λ* increases firstly and tends to be stable during triaxial loading. At the same time, a large local fluctuation is observed. With the increase of *λ*, the evolution of the valence relative to axial strain tends to be stable gradually, and the local fluctuation is also weakened. When the axial strain reaches a certain level, the valence of RVEs with different *λ* reaches a size-independent critical state.

For the RVEs with various particle sizes in Figure 11(b), a similar evolution of the valence is observed at small *λ*. However, with the increase of *λ*, the original evolution of valence tends to be stable in the whole loading and finally tends to an evolution of decreasing at first followed by stabilizing. When the axial strain is large enough, all RVEs reach a size-independent critical state. A comparison with Figure 11(a) suggests that there is no obvious difference in the critical valence of RVEs between the two groups. This indicates that, regardless of size coefficients and their origin, the valence of different RVEs finally reaches a unified critical state.

As shown in Figure 11(c), under the 5 g gravity, the valence of RVEs with various heights shows a similar evolution to that at 0 g at small *λ*. At the larger *λ*, the valence also tends to an evolution of decreasing at first followed by stabilizing. A comparison with Figure 11(a) shows that the critical valence under 5 g gravity is slightly smaller than that under 0 g. Therefore, the final critical state is gravity-related.

For RVEs with various particle sizes in Figure 11(d), the similar evolution of valence is observed. The two groups of RVEs with different origin of size effect finally reach the same size-independent critical state. Again, the observed critical valence in Figure 11(d) is smaller than that in Figure 11(b), which further proves the gravitational dependence of the critical state.

##### 5.2. Orientation

Figure 12 gives the evolution of orientation degree of void cells in RVEs during triaxial loading under 0 g gravity and 10 kPa confining pressure, and both the size effects caused by the variation of RVE height and particle size are demonstrated. The figure is marked by the parameters of gravity-*λ*-axial strain. It can be found from Figure 12(a) that the orientation degree shows a large local fluctuation at small *λ* during isotropic consolidation under 0 g gravity. At the same time, the number of void cells orientated in the vertical direction is slightly more than that in the horizontal direction. With the increase of *λ*, the number of vertically orientated void cells decreases gradually, and the local fluctuation is also dismissing. This indicates the orientation of void cells tends to be more uniform along with increasing RVE height during isotropic consolidation.

**(a)**

**(b)**

During triaxial loading under 0g gravity, the number of vertically orientated void cells increases significantly. However, the increment reduces with the increase of *λ*. This suggests that the increase of RVE height weakens the vertically orientating tendency of void cells.

After the isotropic consolidation of RVEs with various particle sizes at 0 g gravity, as shown in Figure 12(b), the orientation degree shows a significant local fluctuation and slight vertical distribution at small *λ*. Compared with Figure 12(a), the local fluctuation is smaller, but the vertical distribution is more significant. With the increase of *λ*, both the local fluctuation and vertically orientated void cells decreases. This indicates that the decreasing particle size weakens the vertically orientation of void cells during isotropic consolidation. At the same time, the number of vertically orientated void cells in the group of RVEs with various particle sizes is larger than that in the group with various RVE heights at larger *λ*.

During triaxial loading under 0g gravity, the orientation degree of void cells in RVEs with various particle sizes shows a similar evolution to that in the group with various RVE height. However, compared with Figure 12(a), there is less number of void cells orientated vertically in the former group. Meanwhile, the decreasing rate of vertically orientated void cells relative to *λ* is also lower than that in Figure 12(a), which leads to more void cells being orientated vertically at large *λ* in the group with various particle sizes. The above observations suggest that the orientation of void cells depends on the initial state after isotropic consolidation.

Figure 13 gives the evolution of orientation degree of void cells in RVEs during triaxial loading under 0g gravity and 10 kPa confining pressure. The first row in Figure 13(a) demonstrates the orientation of RVEs at various heights after isotropic consolidation. At the small *λ*, a significant local fluctuation is observed, and the number of void cells orientated vertically is slightly greater than that in the horizontal direction. A comparison with Figure 12(a) shows that when the gravity increases from 0 g to 5 g, the vertically orientated void cells decrease gradually. With increase of *λ*, the vertically orientated void cells show a significantly decreasing, and the number of horizontally orientated void cells finally becomes prevalent. Compared with the first row in Figure 12(a), it can be found that the vertically orientated void cells also decrease with the increase of gravity during consolidation.

**(a)**

**(b)**

During triaxial loading under 5 g gravity, as shown in Figure 13(a), the orientation of void cells shows a similar evolution to that under 0 g gravity. An RVE height dependence of the orientation of void cells is observed. A comparison with Figure 12(a) shows that when the gravity increases from 0g to 5 g, the number of vertically orientated void cells decreases gradually. This indicates that the vertical orientation of void cells is inhibited by gravity during triaxial loading.

For isotropic consolidation of RVEs with various particle sizes under 5 g gravity, as shown in the first row of Figure 13(b), a particle size dependence of the orientation of void cells is observed. Almost an equal number of void cells are orientated in the vertical and horizontal directions at small *λ*. A comparison with Figure 12(a) shows that the number of vertically orientated void cells decreases slightly with the increase of gravity. With the increase of *λ*, the vertically orientated void cells also decrease. This suggests that the decrease of particle sizes results in a reduction in the number of vertically orientated void cells. A further comparison with Figure 13(a) shows that, after consolidated under 5 g gravity, the vertically orientated void cells are prevalent in the group with various RVE height at small *λ*. However, an opposite trend is observed at large *λ*.

Figure 13(b) shows that the size effect caused by particle size has a significant impact on the orientation of void cells during triaxial loading. A similar evolution of the orientation is observed compared with that in Figure 13(b). However, with the increase of gravity, the number of vertically orientated void cells decreases gradually at the same axial strain. This suggests that gravity inhibits the vertical orientation of void cells of RVEs in the group with various particle sizes. A further comparison with Figure 13(a) shows that, when applying a gravity-induced stress gradient, the void cells tend to be oriented more vertically at small *λ* and more horizontally at large *λ*. These observations indicate that the void cells with bigger particles are more sensitive to gravity.

##### 5.3. Elongation Degree

Figure 14 gives the evolution of the elongation degree of void cells in RVEs under the 0 g and 5 g gravity, respectively. For the RVEs with various heights in Figure 14(a), all of the elongation degrees are larger than 1.00 after isotropic consolidation, which suggests an initial anisotropy along the vertical direction. With the increase of *λ*, the anisotropy decreases gradually and finally tends to be stable. The RVEs with various particle sizes in Figure 14(b) show a similar evolution. Regardless of the origin of size effects, the elongation degree shows basically the same value at the same *λ*.

**(a)**

**(b)**

**(c)**

**(d)**

When RVEs with different heights are subjected to triaxial loading, as shown in Figure 14(a), the elongation degree of void cells increases at first followed by a stable state at small *λ*. With the increase of *λ*, the growth rate of the elongation degree reduces gradually. This indicates that the void cells quickly achieve a stable state during triaxial loading at smaller *λ*. All of the elongation degrees reach a unified stable state finally, regardless of the size coefficients.

For the RVEs with various particle sizes in Figure 14(b), a similar change in elongation degree is observed during triaxial loading. A comparison with Figure 14(a) suggests that all RVEs with various heights and particle sizes finally reach a similar elongation degree of about 1.35. This indicates that a size effect and its resource independent critical state for the elongation degree is reached during the triaxial loading under 0 g gravity.

When consolidated under 5 g gravity, as shown in Figures 14(c) and 14(d), the initial elongation degree shows similar evolution to that under 0 g gravity. A comparison with Figure 14(a) shows that, when the gravity increases from 0 g to 5 g, the initial elongation degree shows a decreasing tendency. Furthermore, the RVEs in Figure 14(d) with smaller *λ* show a larger decrease in elongation degree, which indicates that the RVEs with large particles are more sensitive to gravity.

As shown in Figure 14(c), the elongation degree of RVE with *λ* of 40.921 increases at first, decreases then, and becomes stable finally (regardless of local fluctuation) during triaxial loading. With the increase of *λ*, the evolution of elongation degree turns into the increasing-stabilizing style. Under large axial strain, all RVEs finally reach a unified critical state of elongation. A comparison with Figure 14(a) indicates that when the gravity increases from 0 g to 5 g, the critical elongation degree shows a decreasing tendency.

For the RVEs with various particle sizes in Figure 14(d), the elongation degree shows a similar evolution as in Figure 14(c). However, a more significant decrease is observed for the RVE with *λ* of 40.921. A comparison with Figure 14(b) shows that the elongation degree decreases along with the increasing gravity. By further comparing with Figure 14(c), it can be found that the critical elongation degree at the unified sable state is independent of the size effects and their origins.

#### 6. Conclusions

In this work, the granular matter was homogenized into stress gradient continuum by adopting RVE. Two origins of the size effect of RVE were considered numerically, and the size-dependent behaviour of RVE under gravity-induced stress gradient was investigated. The following conclusions are drawn:(1)The gravity-induced stress gradient can be naturally introduced into the higher-order continuum theory by adopting RVE, which provides a consistent framework to study the size-dependent behaviour of granular material under gravity-induced stress gradient.(2)The isotropic consolidation under gravity-induced stress gradient shows no obvious dependency on the size coefficients and their origins. However, the increasing stress gradient leads to an increase in compression.(3)During triaxial loading, the size-dependent evolution of deviatoric stress, internal friction angle, and elastic module show a gravity-induced stress gradient dependency. A large influence of stress gradient is observed in those RVEs with large particle sizes.(4)The gravity-induced stress gradient inhibits volumetric dilatation during triaxial loading, and the origin of the size coefficients has a significant influence on volumetric strain at a small size coefficient.(5)Both the gravity-induced stress gradient and the origin of the size effect affect the evolution of microscopic void cells (including valence, orientation, and elongation degree), and the influence of the stress gradient is closely related to the particle size.

#### Data Availability

The data are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 41902273 and 41772338); the China Postdoctoral Science Foundation (no. 2019M661986); the Natural Science Foundation of Jiangsu Province (no. BK20190637); and the Jiangsu Planned Projects for Postdoctoral Research Funds (no. 2019K194). The corresponding author would like to thank the financial support by the State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology (nos. Z19007 and Z19009).