For geotechnical engineering problems, it is unrealistic and unnecessary to conduct detailed and complete simulation for each solid particle on the micro scale. Establishing a new constitutive relationship under the continuous medium mechanics framework considering the influence of mesostructure characteristics of soil-rock mixture has become an effective way to simulate and predict the shear strength of soil-rock mixture. In fact, the macroscopic deformation and failure of the soil-rock mixture are mainly determined by the rock block (Li et al., 2016; Wang et al., 2020). Therefore, according to the physical and mechanical effects of the interaction of soil particles and rock blocks of the soil-rock mixture at different scales, the soil-rock mixture is divided into soil matrix and rock blocks, and a representative volume element (RVE) is constructed at the mesoscale, which is large enough compared with the local microstructure volume element size but small enough compared with the whole soil-rock mixture layer in the engineering scale. According to previous studies (Zhang et al., 2016; Du et al., 2017; Yang et al., 2021), the soil/rock threshold of 5 mm can be selected. Therefore, particles with diameters smaller than 5 mm are referred to as soil particles; those with diameters larger than 5 mm are referred to as rock blocks. In the RVE, the soil particles interact with the pore water to form a soil matrix, and the rock block is regarded as a rigid sphere, as shown in Figure 1.

It is assumed that the rock block is a rigid sphere and is wrapped by the soil matrix. Each rock block and its adjacent matrix soil form a cubic cell, as shown in Figure 1. Therefore, the volume of the RVE is L ³ and the volume of the rock block is πd ³/6. The side length L can be expressed as (Feng and Fang 2015c; 2016): L = π 6 α 1 3 d (1) where d is the particle size of the rock block; α is the rock block content, and denotes the ratio of the volume of the rock block to the total volume of the soil-rock mixture sample; Here, in order to make the soil-rock mixture sample have the soil-rock cell structure, the length L of the cube cell has to be larger than the diameter d of the equivalent spherical rock block. According to Eq. 1: L = π / 6 α 1 / 3 d > d , i.e., α < π / 6 ≈ 0.52 . Therefore, the applicability of the soil-rock cell model is α < π / 6 ≈ 0.52 .

Different structural levels of soil-rock mixture show different mechanical responses. The rock blocks mainly translate and rotate, while the matrix mainly produces microcracks (Wang et al., 2018; Zhao and Zhou 2020). Moreover, there are extremely complex phase to phase coordinated deformation and interface interaction between the soil matrix and the rock block, which makes the stress transfer, crack evolution and failure characteristics of the soil-rock mixture different from those of the pure soil and rock.

When subjected to external loads, the soil-rock mixture, whose skeleton is composed of mineral particles of cross grain group scale, shows that the contact bond and cementation bond between the matrix particles are broken and microcracks are generated at the microscale (Feng and Fang, 2016); at the mesoscale, it is shown that the plastic shape distortion of the soil matrix is induced by the incompatible deformation between the soil matrix and rock block in the soil-rock cell; at the macroscale, it shows complex engineering mechanical behaviors such as nonlinear and zigzag stress-strain relationship of the soil-rock mixture (Xu et al., 2011). Therefore, the mechanical behaviors of soil-rock mixtures at different scales are not only coupled and related to each other, but also have certain independence. As a meso RVE reflecting the basic mechanical effects and particle characteristics of soil-rock mixtures, the soil-rock cell is the key factor to establish a multiscale framework for the soil-rock mixture.

For the soil-rock cell, the natural granular and structural characteristics of the soil-rock mixture at the microscale cause its internal deformation to show significant discontinuity and nonlinearity, which is manifested as the gradient phenomenon of shear strain of the soil matrix connecting the rock block (Feng and Fang 2016). Strain gradients occur when the microstructures of the material undergo plastic shape distortion. Therefore, the strain gradient theory is used to describe the plastic shape distortion inside the soil-rock cell in this study, and then the strain energy of the soil-rock cell can be expressed by the strain and the strain gradient (Fang 2014; Feng and Fang 2015c): U ∼ = a 2 G m τ s 2 = a 2 G m τ R V E m 2 + U η (2) and U η = b G η η 2 / 2 (3) τ R V E m = G γ R V E m (4) where U ∼ is the strain energy of the soil-rock mixture; γ s is the yield shear strain of the matrix; η is the strain gradient of the soil matrix adjacent to the rock block, and can be defined as the second order gradient of displacement of the RVE (Gao et al., 1999); U η is the strain energy corresponding to the strain gradient η generated by the rock block; τ s and τ R V E m are the shear yield stress of the soil-rock mixture and the soil matrix inside the cell respectively; G m is the shear modulus of the matrix; G η is the shear strain gradient modulus; b and a are the nonlinear coefficients of the soil-rock cell and soil matrix respectively, the physical meaning of a and b is shown in Figure 2, where A _arc is the area of the sector (OAB); A _tri is the area of the triangle (OAB); γ B is the yield shear strain; η B is the strain gradient at yielding; τ η is higher-order shear stress, and is work conjugate with strain gradient η ; τ A is the yield stress; and τ η A is the higher-order shear stress at yielding.

A concise theoretical model for calculating the shear yield stress of the meso representative volume elements can be obtained from Eq. 2–4: τ s = τ R V E m 1 + l e 2 ⋅ η 2 (5) and l e = G m G η b a τ R V E m 2 (6) and l e = τ s 2 τ R V E m 2 − 1 η 2 (7) where l e is the intrinsic length of the soil-rock mixture, which can be determined by specially designed tests. For the strain gradient η , considering the small size of cube soil-rock mixture representative volume element (RVE), it can be approximated by the first-order difference of mesoscopic representative volume element strain for simplicity (Fang 2014), namely: η = Δ γ L − d / 2 ≈ 2 L − d γ x = L / 2 − γ x = 0 ≈ 2 γ L − d (8) where L is the side length of the volume element represented by the meso cube of the soil-rock mixture; d is the particle size of the rock block; α is the rock block content; γ denotes the shear strain of the RVE. According to Eq. 8, the strain gradient η of the RVE of the soil-rock mixture can be further expressed by combining of Eqs 1, 8: η = 2 γ π / 6 α 1 / 3 − 1 d (9)

In the soil-rock mixture, the diameter of the rock blocks is varying, therefore, the average strain gradient of the soil-rock mixture can be calculated by the grading function of the rock blocks, one may find: η ¯ = ∫ d m d M g ′ d η δ d / ∫ d m d M g ′ d δ d (10) where g ′ d is the first derivative of g d , and g d is the gradation function of rock block. d M and d m are the maximum and minimum particle sizes of the rock blocks in the soil-rock mixture, respectively.

Eq 5 is a concise theoretical model of predicting the multiscale coupled shear strength of soil-rock mixture. However, the mechanical response of the soil matrix inside the soil-rock cell is different from that of the pure matrix due to the influence of the rock block, and the stress-strain relationship of the matrix inside the soil-rock cell is difficult to determine. According to the purpose of dividing the soil-rock mixture in to soil matrix and rock block, for the practical application of Eq. 5, the following two problems still have to be solved through theoretical analysis: (1) how the shear stress and shear displacement relationship τ R V E m = f 0 s of matrix in soil-rock mixture can be approximately expressed by the mechanical response τ c = f c s of pure matrix material, where τ c is the shear stress of pure matrix material, as shown in Figure 3; (2) The intrinsic length scale of the soil-rock mixture is a key parameter to describe scale characteristics of materials, which is related to the content and gradation of rock blocks, and the soil matrix properties; how to realize the relationship between the intrinsic length scale and the classical finite element equation? In the present work, based on the Eshelby-Mori-Tanaka equivalent inclusion and average stress principle, we realizes the quantitative correlation between the matrix shear stress in the soil-rock mixture and the mechanical response of the pure matrix material, and establishes the constitutive relationship and finite element equation of the soil-rock mixture with embedded intrinsic scale parameters through the couple stress theory.

Hu et al. (2013) and Fang et al. (2013) have respectively derived the effective elastic modulus of pebble soil and the equivalent effect gradient of cohesive soil by using the Eshelby equivalent inclusion principle, and obtained good theoretical results. According to Eshelby equivalent inclusion principle (Eshelby 1957), the shear stress of the soil matrix and the rock block in the representative volume element (RVE) of the soil-rock mixture can be expressed as: τ R V E m = τ c + τ ∼ = G m γ c + γ ∼ (11) τ R V E r = τ c + τ ∼ + τ p t = G m γ c + γ ∼ + γ p t − γ * = G r γ c + γ ∼ + γ p t (12) where τ R V E r and τ R V E m are the shear stresses of rock block and soil matrix in the representative volume element of soil-rock mixture, respectively; τ c and γ c are the shear stress and shear strain of the pure soil matrix material (without any rock block, α = 0 ); τ ∼ and γ ∼ are respectively the perturbation shear stress and the perturbation shear strain of the soil matrix produced by the rock block; τ p t and γ p t are the disturbed shear stress and the disturbed shear strain of the rock block, respectively, which can be determined by the intrinsic strain γ * of the rock block (Hill 1963; Weng 1984); G m and G r are the material constants of the soil matrix and rock block, respectively. τ p t = T τ * (13) γ p t = S γ * (14) G m = 3 G 0 = 3 E 0 2 1 + υ 0 − 1 (15) G r = 3 G R = 3 E R 2 1 + υ R − 1 (16) where S is Eshelby constant, T is hill constant, and S + T = 1 (Hill 1963; Weng 1984); G 0 and G R are the shear modulus of the soil matrix and the rock block, respectively; E 0 and E R are the elastic modulus of the soil matrix and the rock block, respectively; υ 0 and υ R are the Poisson’s ratio of the soil matrix and the rock block, respectively. According to the Eshelby equivalent inclusion principle (Weng 1984), the relationship between the perturbation strain γ ∼ and the intrinsic strain γ * can be established by the Eshelby constant S, one may find: γ ∼ = S γ * (17)

For ellipsoidal rock block (spherical is its special case), Eshelby constant can be expressed as: S = 2 4 − 5 υ 0 / 15 1 − υ 0 (18) where υ 0 is the Poisson’s ratio of the soil matrix.

Combining Eqs. 11–14, one may find: τ * = − G m γ * (19)

Substituting Eqs 11, 17 into Eq. 12, yields: τ p t = G m S γ * − G m γ * (20)

According to the Mori Tanaka average stress method (Mori and Tanaka 1973), the average stress remains unchanged, the shear stress of the soil matrix can be expressed as: τ c = α τ c + τ ∼ + τ p t + 1 − α τ c + τ ∼ (21)

According to Eq. 21, one may find: τ ∼ = − α τ p t (22)

Substituting Eq. 22 into Eq. 20, yields: τ ∼ = − G m α S − 1 γ * (23)

Substituting Eq 12, 21 into Eq. 11, one may find: γ * = b a γ c (24) and a 0 = − G r − G m S 1 − α + α + G m (25) b 0 = G r − G m (26)

Substituting Eqs 23, 24 into Eq. 11, yield: τ R V E m = τ c + τ ∼ = G m γ c + G m α S − 1 γ * = G m γ c + G m α S − 1 b 0 a 0 γ c = G m γ c 1 + α S − 1 b 0 a 0 (27)

It can be seen from Eq. 27 that the occurrence of rock block causes stress concentration in the adjacent soil matrix, and the corresponding stress concentration coefficient χ can be expressed as (Zhao and Zhou 2020): χ = τ R V E m τ c = 1 + α S − 1 b 0 a 0 (28)

Combining Eqs 27, 28, the shear stress-shear displacement relationship between the soil matrix in the RVE and the pure soil can be established: f 0 s = χ f c s (29) where τ c = f c s is the shear stress-shear displacement relationship of pure matrix material, and τ R V E m = f 0 s and τ c = f c s .

Eq 28 gives the stress concentration coefficient of the soil matrix in the soil-rock mixture with an assumption that all the rock blocks have the same diameter. For the actual soil-rock mixture, the size of the rock block is distributed in the sample. In order to consider the influence of the size of rock blocks on the strength characteristics of soil-rock mixtures, the number density function f d of rock blocks is introduced (Feng and Fang 2016). Its physical meaning is the percentage of the number of rock blocks with equivalent diameter d in the total number of rock blocks in the unit volume of soil-rock mixture, which can be determined by the grading function g d of rock blocks. Therefore, f d can be expressed as (Feng and Fang 2015a; Feng and Fang, 2016): f d = 6 α g ′ d d d N ⋅ π d 3 (30) where N is the total number of rock blocks in the unit volume of soil-rock mixture; d d is the differential of size of rock block.

The strength of the soil-rock mixture sample is provided by the soil matrix and rock blocks. Therefore, its shear strength can be expressed as: τ R V E f = A r τ R V E r + A m τ R V E m A = A r τ R V E r + A − A r τ R V E m A = τ R V E m + A r A τ R V E r − τ R V E m (31) where τ R V E f is the shear strength of the soil-rock mixture; A r is the area of rock block in the RVE; A m is the area of soil matrix in the RVE; A is the area of the RVE, A = L 2 = π / 6 α 2 / 3 d 2 . According to Eq. 30, the area of rock block in the soil-rock mixture sample can be expressed as: A r = ∑ N k A k r = ∑ N ⋅ N k N A k r = ∑ N ⋅ f d A k r (32) where N k is the number of rock block with particle size d _k; A k r is the area of rock block with the size of d _k, A k r = π d k 2 / 4 .

Substituting Eq. 30 into Eq. 32, yields: A r = 3 2 α ∫ d min d max g ′ d d d d (33)

For the case where the distribution range of rock block size is small, the grading curve can be approximately expressed as a linear function of the size of rock block (Feng and Fang 2015a), one may find: g ′ d = C (34) where C is a constant.

Simultaneous Eq. 31 into Eq. 33, yields: τ R V E f = τ R V E m + 3 α ∫ d min d max g ′ d d d d 2 π / 6 α 2 / 3 d 2 τ R V E r − τ R V E m (35)

Combination Eqs 12, 27, 29, 35, the shear strength of the soil-rock mixture can be further expressed as: τ R V E f = χ τ c f + 3 α ∫ d min d max g ′ d d d d 2 π / 6 α 2 / 3 d 2 ⋅ G r G m 1 + S + α − S α b a − χ τ c f = Ω τ c f (36) and Ω = χ + 3 α ∫ d min d max g ′ d d d d 2 π / 6 α 2 / 3 d 2 ⋅ G r G m 1 + S + α − S α b 0 a 0 − χ (37) where τ c f is the yield stress of pure matrix soil; Ω is the equivalent strength enhancement coefficient of the soil matrix around the rock blocks.

Combination Eqs 7, 36, the intrinsic length scale of soil-rock mixture can be further expressed as: l e = τ s 2 Ω 2 τ c 2 − 1 η 2 (38)

The classical finite element method treats the soil-rock mixture as a uniform continuous medium. However, at the microscale, the soil-rock mixture has significant inhomogeneous and discontinuous mechanical characteristics (Zhao et al., 2021). Based on the micro motion characteristics of rock blocks, a shear strength theoretical model considering the rotation effect of rock blocks is established based on the soil-rock cell model. According to the coupled stress theory (Toupin 1962; Borst 1991), under the action of shear stress, the soil-rock cell generates three translational displacements u x , u y , u z , and three rotational displacements ω x , ω y , ω z , and the corresponding work conjugate stress and coupled stress are σ i j and m i j , respectively. Therefore, the strain of soil-rock cell can be written as (Fang and Li 2016): γ i j = ∂ u i ∂ x j + e i j k ω k = u i , j + e i j k ω k (39) γ i j = ∂ u i ∂ x j + e i j k ω k = u i , j + e i j k ω k (40)

The generalized equivalent strain of soil-rock cell (RVE) can be expressed by rotational gradient, namely: γ e = 2 3 γ i j ′ γ i j ′ + l e 2 η i j ′ η i j ′ 1 2 (41) where γ i j ′ = γ i j − γ m δ i j is strain deviation; γ m = γ i i / 3 is spherical strain; η i j ′ = η i j − η m δ i j is curvature tensor deviation; η m = η i i / 3 is spherical curvature; η i j = ω i , j is curvature tensor; The intrinsic length scale l e of the soil-rock mixture can be determined by Eq. 38.

Accordingly, the generalized equivalent stress q e (Mindlin 1963) of the RVE can be expressed as: q e = 3 2 s i j s i j + l e − 2 m i j ′ m i j ′ 1 2 (42)

Where: s i j = σ i j − p δ i j is stress deviation; p = σ i i / 3 is spherical stress; m i j ′ = m i j − m δ i j is couple stress deviation; m = m i i / 3 is spherical couple stress.

Eq 39 can be written in the following matrix form: ε = B u (43) ε = γ 11 , γ 22 , γ 33 , γ 12 , γ 13 , γ 21 , γ 23 , γ 31 , γ 32 , η 11 , η 22 , η 33 , η 12 , η 13 , η 21 , η 23 , η 31 , η 32 T (44) u = u x , u y , u z , ω x , ω y , ω z T (45)

Where: B is the flexibility matrix; ε is the strain array; u is the displacement array.

Since the translation and rotation of the RVE of the soil-rock mixture are considered in the multiscale coupled finite element method, for the plane problem, three degrees of freedom need to be set, namely U1, U2 and R3, where U1 and U2 represent the translational displacement of the element and R3 represents the rotational displacement of the element; the material parameters of the multiscale coupled finite element method include the multiscale coupling parameters such as l e , G m and K m in addition to the conventional macroscopic physical and mechanical parameters such as the density ρ of the soil, the elastic modulus E , Poisson’s ratio υ , the shear modulus G and the bulk modulus K . There are eight material property constants in total. The specific material property constant value shall be given in the INP file of ABAQUS finite element program to be transmitted to the user subprogram of UEL to realize the identification and solution of the multiscale coupled finite element program prepared by ABAQUS finite element software. The programming framework is shown in Figure 4.

The soil-rock mixture used in this paper is taken from the north of Guangzhou, China and mainly consists of silty clay and sandy mudstone. The soil-rock mixture was taken out from the excavation area first and fully turned over by a 20-t excavator to make the rock blocks in the soil-rock mixture evenly distributed in the soil matrix. The mixed soil and rock blocks were compacted in layers by a 14-t road roller with a compaction ratio of 92%. The compacted soil-rock mixture forms a test embankment. As shown in Figure 5, the large-scale direct shear test samples of soil-rock mixture were prepared inside the embankment. The physical and mechanical properties of soil-rock mixture were shown in Table 1. The maximum particle size of rock block is 50 mm, and the gradation curve of rock blocks is shown in Figure 6. For the water content, specific gravity and liquid plastic limit test of the soil-rock mixture, the tests were carried out after the solid particles with diameters greater than 2 mm were removed from the soil-rock mixture.

The large-scale direct shear test device was mainly composed of a shear frame (the inner dimension of the shear frame is 500 mm × 500 mm × 250 mm), reaction force supplying and loading application systems and stress and displacement acquisition system. The photographs of the equipment installation of the large-scale direct shear test were shown in Figure 7.

The normal stress of each group of tests was 36, 72, 108 and 144kpa respectively, and each level of normal stress was loaded in three levels. The large-scale direct shear test scheme was shown in Table 2. The tests were conducted according to the Specification of soil test (GB/T50123-2019) Ministry of Water Resources of the People's Republic of China, 2019.

The shear stress-shear displacement curves of the samples were shown in Figure 8.

It can be seen from Figure 8 that the shear stress-shear displacement curves of pure matrix samples were relatively smooth; for the soil-rock mixture samples, the shear stress-shear displacement curve had serrated characteristics, and there was a discontinuous jump phenomenon of stress, and the overall mechanical response was strain hardening. Moreover, the shear strength of the soil-rock mixture samples increases with an increase in the content of the rock blocks under the all levels of normal stress conditions.

The above phenomenon that the shear strength and deformation characteristics of the soil-rock mixture changed with the change of the rock block content reflected the key influence of the rock blocks on the mechanical responses of soil-rock mixture.

The parameters of the proposed multiscale finite element method are shown in Table 4. Herein, G is equal to G ^m , G = E / 2 1 + υ 0 , K = E/[3×(1-2ν ⁰)]; G _m and K _m are determined by Eqs 51, 52, respectively.