A multiscale finite element method for soil-rock mixture

Soil-rock mixture is a complex multi-phase composite geotechnical material, and its strength is determined by the physical properties of constituent multi-phase materials and their coupling mechanical response between different phases of materials. Based on the Eshelby-Mori-Tanaka equivalent inclusion average stress principle, a theoretical model of multi-scale coupled shear strength of soil-rock mixture considering the interaction effect of rock block and soil is established, and the rotational freedom reflecting the microscopic motion details of rock block is introduced. Moreover, a multi-scale coupled constitutive relationship of soil-rock mixture is derived and compiled into a multi-scale finite element program. Based on the large-scale direct shear test of soil-rock mixture, the model parameters of the multi-scale finite element method are determined, and then the multi-scale finite element program is used to simulate and predict the cross-scale deformation process of the soil-rock mixture slope. The results show that the multi-scale finite element method can effectively describe the influence of the mechanism of the micro motion characteristics of the soil-rock mixture on the macro mechanical response, and can effectively overcome the pathological mesh-dependency of the classical finite element method; the rotation displacement of the rock block is mainly concentrated within the shear zone of the slope. The maximum rotational displacement of rock blocks inside the soil-rock mixture slope is 40.7°, and the rotational displacement of rock blocks outside the shear zone is about 0°. The physical mechanism of the cross scale evolution of the shear band of the soil-rock mixture slope is that: the rotation of the rock blocks weakens the strain transmission ability between the rock block and the matrix soil, thus forming the concentration and development of the plastic strain, and finally leading to the penetration of the shear bands of the slope and the overall sliding failure.


Introduction
The finite element method is an effective numerical method for solving boundary value problems in solid mechanics (Strouboulis et al., 2001), and has been widely applied in geotechnical engineering problems such as foundation pit design, tunnel engineering, foundation treatment and slope stability analysis (Collinswilliams 2015;Liu et al., 2015;Cai et al., 2020;Mikkelsen et al., 2022). However, the soil-rock mixture has the heterogeneous and discontinuous nature of multiphase of solid particles at different scales, water and air. When the classical elastic-plastic finite element method based on the principle of continuum mechanics is used to analyze and calculate geotechnical engineering problems, it often produces large deviations, such as the pathological mesh-dependency of the prediction results and the difficulty in convergence of the solution process (Broumand and Khoei 2013). The mechanical properties of the particle interface and the particle interior of the soil-rock mixture are completely different. Even under uniform stress condition, the deformation across the interface will be discontinuous. Therefore, the deformation and strength characteristics of soil-rock mixture are significantly different from those of the continuous medium. The discontinuity of the deformation of the soil-rock mixture is caused by the discontinuity of its internal structure and the physical properties of particle characteristics, that is, the natural granularity and structure of the soil-rock mixture make it a discontinuous medium. The mechanical properties of the discontinuous medium are usually highly nonlinear, with nonlinear characteristics such as deformation localization, parameter sensitivity and system catastrophe (Fang and Li 2016). The classical soil mechanics treats the soil-rock mixture as a homogeneous continuous medium, which embeds the assumptions of continuity, material homogeneity and physical mechanism consistency, and correspondingly neglects the natural granularity, structure, multiphase and mechanical mechanism couplings of the soil-rock mixture at different scale structure levels. Only for small uniform deformation, it is possible to obtain approximately reasonable theoretical results, and significant errors will occur in the analysis and prediction of the generation and development of large deformation shear bands closely related to the particle size and structure of soil-rock mixture (Tang et al., 2018). For example, when the classical elasticplastic finite element method based on the principle of macro continuum mechanics simulates the evolution of shear bands, its prediction results strongly depend on the size of the mesh and show significant pathological mesh-dependency (Soltani and Maekawa, 2015;Chang et al., 2021). Therefore, the simulation and prediction results of the mechanical response of the soil-rock mixture with natural multiscale coupling characteristics by the classical elastic-plastic finite element method cannot conform to the physical reality and lose objectivity. The reason is mainly reflected in two aspects. At the physical level, the classical elastic-plastic finite element method regards the soil-rock mixture as a uniform continuous medium, and neglects the influence of the mechanism of the physical details and motion characteristics of the soil-rock mixture at the micro level on the macro mechanical response; At the mathematical level, when the softening constitutive relation is used to describe the cross-scale deformation process of the soilrock mixture slope in the classical elastic-plastic finite element method, its control partial differential equation loses its positive definiteness (Arriaga et al., 2015).
With the development of the research in the interdisciplinary field of computer science and computational mechanics, the finite element method has made great progress in simulating and predicting the multiscale coupled mechanical response of the deformation failure process of soil-rock mixtures. Plenty of previous studies have focused on the use of finite element method for prediction of the mechanical responses of the soilrock mixture. Recently, some finite element methods that can describe the microstructure characteristics of the soil-rock mixture by dividing it into rock blocks and soil matrix have appeared (Guo and Zhao, 2015;Guo et al., 2021). The ideas of these macro and micro finite element analysis methods can be divided into two categories. One is to theoretically construct various representative elements that can reflect the microstructure characteristics of rock blocks and soil matrix. For example, Goodman et al. (1968) constructed a joint element that describes the discontinuity of rock mass and evaluated the impact of mechanical properties of joints on the strength and stability of rock mass; Ling and Ye (2005) constructed a plane coordinated manifold soil element to describe the consolidation characteristics of soil and predicted the settlement and pore water pressure distribution of the riverbed; Li and Wan (2011) constructed a micropolar element reflecting the rotational displacement of soil particles and analyzed the micromechanism of soil strain localization; Liang and Zhao (2019) constructed a representative element reflecting the micro structure characteristics of soil particles and simulated the multiscale coupled evolution process of large deformation of soil-rock mixture. The other is to directly capture the micro physical details of the soilrock mixture by various imaging methods and reconstruct them into files recognizable by the finite element program by digital image processing technology. For example, Xu et al. (2008) obtained the material information of the soil-rock mixture by using an optical camera and analyzed the influence of the micro structure characteristics of the soil-rock mixture on its macro strength characteristics through the combination of digital image processing technology and finite element method; Zheng et al. (2014) used nuclear magnetic resonance technology to characterize the material composition and combined with finite element method to simulate and predict the nonlinear mechanical response of materials; Bubeck et al. (2017) used CT scanning technology to determine the shape characteristics of pores in rocks, and combined with finite element method to study the influence of different pore shapes on the strength of geomaterials.
Although the above-mentioned researches have preliminarily explored the establishment of finite element method considering the micro motion details and structural characteristics of soil-rock mixture and achieved certain results, the model parameters involved do not have clear physical meanings. Meanwhile, the theoretical basis of these research results is still based on the phenomenological uniform continuum mechanics at the macroscale, which fails to break through the limitation that the classical constitutive model of soil-rock mixture cannot describe the coupling of mechanical responses of soil-rock mixture at different scale structure levels. Therefore, it is impossible to scientifically simulate and predict the cross-scale coupling correlation characteristics of the influence of the micro motion details of The representative volume element of soil-rock mixture. rock and soil particles at different scale on the macro mechanical response in the deformation process of the soil-rock mixture from the physical nature and physical mechanism. Therefore, Fang (2015a, 2015b) established a multiscale soil cell model of cohesive soil containing sand based on the coupled stress theory for the micro motion details of clay particles and sand particles, but only gave a complex theoretical solution, which has not been promoted to the engineering application level. In this paper, based on the idea of the soil cell model, the multiscale coupled constitutive relationship embedded the rotational freedom of soilrock mixture is derived. Using the user subprogram interface of ABAQUS finite element program UEL, a multiscale coupled finite element program considering the rotational freedom of the element is compiled to quantitatively analyze the influence of rock block content on the cross-scale evolution process of the soil-rock mixture slope, In order to lay a solid foundation for the direct application of multiscale coupled finite element method to solve geotechnical engineering problems closely related to particle size characteristics of the soil-rock mixture.
2 Theoretical basis of multi-scale mechanics of soil-rock mixture 2.1 Mesoscopic representative volume element (RVE) of soil-rock mixture 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 soilrock 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 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 soilrock mixture layer in the engineering scale. According to previous studies 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 3 and the volume of the rock block is πd 3 /6. The side length L can be expressed as (Feng and Fang 2015c;: 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.

Multi-scale coupled shear strength theoretical model of soil-rock mixture
Different structural levels of soil-rock mixture show different mechanical responses. The rock blocks mainly translate and rotate, while the matrix mainly produces microcracks 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 soilrock 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): and Frontiers in Materials frontiersin.org whereŨ 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 τ m RVE 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: and and 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: 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: In the soil-rock mixture, the diameter of the rock blocks is varying, therefore, the average strain gradient of the soil-rock  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 soilrock 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 τ m RVE 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:

Eshelby-Mori-Tanaka equivalent inclusion average stress principle
where τ r RVE and τ m RVE 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; τ pt and γ pt 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. τ pt Tτ* (13) 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γ* For ellipsoidal rock block (spherical is its special case), Eshelby constant can be expressed as: where υ 0 is the Poisson's ratio of the soil matrix. Combining Eqs. 11-14, one may find: Substituting Eqs 11, 17 into Eq. 12, yields: 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: According to Eq. 21, one may find: Substituting Eq. 22 into Eq. 20, yields: Substituting Eq 12, 21 into Eq. 11, one may find: and Substituting Eqs 23, 24 into Eq. 11, yield: Frontiers in Materials frontiersin.org 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): Combining Eqs 27, 28, the shear stress-shear displacement relationship between the soil matrix in the RVE and the pure soil can be established: where τ c f c (s) is the shear stress-shear displacement relationship of pure matrix material, and τ m RVE = 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): where N is the total number of rock blocks in the unit volume of soilrock mixture; dd 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: where τ f RVE 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: where N k is the number of rock block with particle size d k ; A r k is the area of rock block with the size of d k , A r k πd 2 k /4. Substituting Eq. 30 into Eq. 32, yields: 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: where C is a constant. Simultaneous Eq. 31 into Eq. 33, yields: Combination Eqs 12, 27, 29, 35, the shear strength of the soilrock mixture can be further expressed as: and where τ cf 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: 3 Multiscale coupled finite element method based on the soil-rock cell model 3.1 Stress-strain relationship of soil-rock cell model 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 . 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 σ ij and m ij , respectively. Therefore, the strain of soil-rock cell can be written as (Fang and Li 2016): Frontiers in Materials frontiersin.org The generalized equivalent strain of soil-rock cell (RVE) can be expressed by rotational gradient, namely: 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: Where: s ij σ ij − pδ ij is stress deviation; p σ ii /3 is spherical stress; m ij ′ m ij − mδ ij is couple stress deviation; m m ii /3 is spherical couple stress.
Eq 39 can be written in the following matrix form: Where: [B] is the flexibility matrix; ε { } is the strain array; u { } is the displacement array.

Elastic stress-strain relationship of soil-rock mixture
The elastic strain energy density of soil-rock mixture is defined as w(γ ij ′ , γ m , η ij ′ , η m ), where the stress deviation is s ij , the couple stress deviation is m i j ′ , the spherical stress p and the spherical couple stress m are respectively work conjugate with the strain deviation γ ij ′ , the curvature tensor deviation η ij ′ , the spherical strain γ m , and the spherical curvature η m , then the strain energy density increment of soil-rock mixture can be expressed as: The elastic stress-strain relationship (constitutive relationship) of the soil-rock mixture can be obtained from Eq. 46: For elastic strain energy, according to the definition of generalized equivalent strain given in Eq. 41, the strain energy of the soil-rock mixture can be simplified as w(γ e , γ m , η m ). Meanwhile, γ e , γ m and η m are work conjugate with q e , p and m, respectively. Therefore, the strain energy increment of the soil-rock mixture can be expressed as: δw γ e , γ m , η m q e δγ e + pδγ m + mδη m where G is the shear modulus; G m is the curvature shear modulus, which can be calculated and determined by Eq. 51 according to the intrinsic length scale l e that can be determined by Eq. 38. The relationship between q e and γ e can be determined by conventional triaxial compression test of soil-rock mixture. Accordingly, the relationship between elastic spherical stress and spherical strain of soil-rock mixture and the relationship between spherical couple stress and spherical curvature η m can be expressed as: where K is the elastic bulk modulus, which can be determined by conventional triaxial compression test; K m is the curvature bulk modulus, which can be calculated by Eq. 53 according to the intrinsic length scale l e . Write Eq. 47 in matrix form:  where [D e ] is the elastic modulus matrix of soil-rock mixture. It can be seen from Eq. 54-58 that the elastic stress-strain relationship of soil-rock mixture is a nonlinear stress-strain relationship.

Plastic stress-strain relationship of soil-rock cell model
In the soil-rock cell model, the rotation of the rock block will cause significant plastic shape distortion of its adjacent matrix. Therefore, the elastic constitutive relationship established by Eq. 54-58 is now extended to the elastoplastic stress-strain relationship. For simplicity, this paper introduces the Von Mises yield criterion based on the hardening (softening) law: f q e − σ y 0 where σ y is the yield stress of the soil-rock mixture. The yield stress of soil-rock mixture σ y can be linear or hyperbolic softening function of equivalent plastic shear strain γ where a 0 and b 0 are constants. The elastic stress-strain relationship expressed by Eq. 54 is rewritten in an incremental form: Introducing associated flow laws: where A is the hardening modulus of soil-rock mixture, which can be expressed as: Substituting Eq. 65 into Eq. 62, the plastic strain expression of the soil-rock cell model can be expressed as: Combining Eqs 61, 67, the elastoplastic stress-strain relationship of soil-rock mixture can be obtained: where [D ep ] is the elastic-plastic modulus matrix of soil-rock mixture; the plastic modulus matrix of soil-rock mixture [D p ] can be expressed as:

FIGURE 4
The programming framework of the multiscale finite element method based on ABAQUS software.

Element finite element matrix based on soilrock cell model
The multiscale coupled finite element method based on the soilrock cell model takes the soil-rock cell (RVE) as the basic element, and considers the translational displacement u i and rotational displacement ω i of the element, as well as the curvature tensor η ij and the couple stress tensor m ij . Therefore, the displacement u { }, stress σ { }, strain [D ep ] of the RVE based on the soil-rock cell model are consistent with the framework of the general classical elasticplastic finite element method, but the physical meanings and constitutive model are essentially different. In this paper, the secondary development platform provided by ABAQUS finite element software is used to redefine and encode the multiscale coupled finite element, and the multiscale coupled finite element program that meets the requirements of the multiscale theoretical framework of the soil-rock cell model and can be read and solved by ABAQUS finite element software is programmed.
The element displacement function of the multiscale finite element method based on the soil-rock cell model can be written as: where [N e ] is the shape function matrix; u e { } is the translational displacement and rotational displacement array of nodes. Using the geometric equation Eq. 43, the strain expression can be obtained:  According to the principle of virtual work, one may find: where f e and p e are the volume force and area force of the RVE, respectively. By substituting the stress-strain relationship of Eq. 68 and the stress-strain relationship of elastic (elastic region) or plastic (plastic region) into Eq. 73, the stiffness equation of the element can be expressed as: For the step loading (incremental method), the stiffness equation of the RVE can be expressed as: where F e { } is the element node load array; [K e ] is the element stiffness matrix.

Stress-strain relationship of soil-rock cell model
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.

Larger scale direct shear test for determining the intrinsic length scale l e
There are eight material property parameters in the multiscale coupled finite element method. The elastic modulus E, Poisson ratio, shear modulus G and bulk modulus K can be determined by the triaxial compression test of the soil matrix; The micro parameter l e can be determined by the large-scale triaxial compression test or large-scale direct shear test conducted on the soil-rock mixture samples; The density of soil-rock mixture can be determined by soil density test. The curvature

Test materials
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 largescale 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.

FIGURE 10
The model of the simulated soil-rock mixture slope.
Frontiers in Materials frontiersin.org

Test device
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.

Test scheme
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

Test results and analysis
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.

Comparison between theoretical prediction and experimental results
It can be seen from Eqs 5, 9; Eq. 38 that the theoretical model (soil-rock cell model) of multiscale coupled shear strength of soilrock mixture is embedded with the intrinsic scale l e and strain gradient η. Therefore, it can simulate and predict the multiscale shear strength of the soil-rock mixture.
The matrix equivalent stress enhancement coefficient Ω is related to the shear modulus, Poisson's ratio and the content of rock blocks, and can be determined by Eq. 37.
The intrinsic length scale of soil-rock mixture is related to the equivalent stress enhancement coefficient Ω and rock block content, and can be determined by Eq. 38.
The strain gradient of the soil-rock mixture is related to the content and size of the rock blocks, which can be determined by Eq. 10.
The elastic modulus of the soil matrix is determined by conventional triaxial compression test of the pure matrix    (Zhang et al., 2015;Zhao et al., 2021). The specific calculation of the model parameters of the multiscale coupling soil-rock cell model of soil-rock mixture are shown in Table 3. Among them, the particle size of the rock blocks of soil-rock mixture is 5-50 mm. According to Eq. 5; Table 3, the relationship between shear stress and shear displacement of soilrock mixture can be predicted. The comparison between theoretical results and experimental results is shown in Figure 9. It can be seen from Figure 9 that the multiscael soil-rock cell model of soil-rock mixture based on the meso physical mechanism can simulate and predict the shear stress-shear displacement relationship and the shear strength of soil-rock mixture. The relative error (|τ th − τ te |/τ te ) between the theoretical results τ th and the experimental results τ te of the shear strength of the soil-rock mixture are less than 10%, which is acceptable.
6 Multiscale finite element method for simulation of soil-rock mixture slope Figure 10 shows the soil-rock mixture slope with a rock block content of 45%, with a slope height of 10 m and a slope ratio of 1: 1. In the present work, the main physical and mechanical parameters of the soil-rock mixture are: elastic modulus E = 50MPa, Poisson's ratio υ = 0.3, and linear softening parameter H = 30 KPa. The size of rock particles used in the slope stability analysis is the same as those used in the larger scale direct shear test, namely, 5-50 mm.
The classical elastic-plastic finite element method and the multiscale coupled finite element method are respectively used to analyse the simulation of the stability of a soil-rock mixture slope. The finite element model is divided into four types of meshes, namely, 2,508, 2020, 1,660 and 1,286 elements. The displacement control method is adopted in loading process. The total vertical displacement loading is s = 0.5 m, the vertical uniformly distributed load on the top of the soil-rock mixture slope is p, and the width of the shear band of the soil-rock mixture slope is d 0 .

Parameter determination of the multiscale coupled finite element method
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] 0 )]; G m and K m are determined by Eqs 51, 52, respectively.
6.2 Finite element simulation results of deformation characteristics of soil-rock mixture slope 6.2.1 Classical elastoplastic finite element simulation results of shear band width of soil-rock mixture slope Figure 11 shows the simulation results of the shear band width of the soil-rock mixture slope by the classical elastic-plastic finite element method. The scale in the figure is the equivalent plastic strain of the soil-rock mixture slope, and the unit is %.
It can be seen from Figure 11 that the shear band width d 0 of mesh one to mesh four is 0.62m, 0.76m, 0.98m and 1.31 m respectively, and the shear band width changes significantly with the change of the mesh size. The width of shear band predicted by mesh four is 111.3% higher than that of mesh 2, when the shear band of soil-rock mixture slope is completely interconnected. Therefore, the simulation results of the classical elastoplastic finite element method display significant pathological mesh-dependency. Figure 12 shows the p~s curve predicted by the classical elasticplastic finite element method. Figure 12 shows that the denser is the mesh of the finite element model, the larger is the softening section of the p~s curve, the stronger is the strain softening behaviour and the earlier is the calculation termination of the finite element program due to the difficulty in model solution. The calculation termination displacement s of mesh 1-4 are 0.168 m, 0.201 m, 0.244 m and 0.319 m, respectively. The difference of the calculation termination displacement caused by different mesh division is up to 89.9%, showing a significant pathological mesh-dependency.

FIGURE 12
The p-s curves simulated by classical elastic-plastic finite element method.
Frontiers in Materials frontiersin.org 6.2.2 Multiscale coupled finite element simulation results of shear band width of soil-rock mixture slope Figure 13 shows the simulation results of the shear band width of the soil-rock mixture slope predicted by the multi-scale coupled finite element method. The scale in the figure is the equivalent plastic strain of the soil-rock mixture slope, and the unit is %.
It can be seen from Figure 13 that the predicted results of the shear band width d 0 of the soil-rock mixture slope by using the proposed multiscale coupled finite element program of mesh 1-4 are 3.32m, 3.43m, 3.48m and 3.51 m respectively. Especially, when the number of elements in mesh four is 48.7% less than that in mesh 1, the prediction result of shear band width is only increased by 5.7%, while under the same condition, the prediction result of classical elastic-plastic finite element method is increased by 111.3%. It can be seen that the multiscale coupled finite element method can effectively overcome the pathological mesh-dependency of numerical simulation prediction results. Figure 14 shows the p~s curves predicted by the multiscale coupled finite element method.
It can be seen from Figure 14 that the p~s curves predicted by the multiscale coupled finite element method completely coincide at the stage of uniform linearity, and also approximately coincide at the stage of plastic and softening deformation. The curves have only a Especially in the softening deformation stage, the multiscale coupled finite element method has no difficulty in solving that occurs in classical elastic-plastic finite element method. It can reproduce the entire nonlinear and softening deformation process of the soil-rock mixture slope. The displacement loading is s = 0.5 m, and the softening characteristics remain unchanged with the mesh changes. Therefore, the prediction results of the p~s curves of the multiscale coupled finite element method depend on the physical and mechanical parameters of the soil-rock mixture and have nothing to do with the mesh density. It can simulate and predict the cross-scale evolution of the strain localization of the soil-rock mixture slope more realistically compared with that of the classical elastic-plastic finite element method.
For the multiscale coupled finite element method, the rotation of rock blocks induces the occurrence of stress concentration at the soil matrix around the soil-rock interface. When the stress concentration intensifies to a certain extent, plastic slip occurs at the interface first, and the plastic strain is constrained in the matrix between the rock blocks. The introduction of the rotational freedom of rock blocks can make the concentrated plastic strain transfer smoothly through the rotation of rock blocks. Moreover, the multiscale coupled finite element method introduces the intrinsic length scale l e with clear physical meanings and scale characteristics to regularize the solution of the governing equations of the multiscale coupled finite element program, and avoid the occurrence of pathological mesh-dependency of prediction results.
6.2.3 Simulation results of rock block rotation displacement of soil-rock mixture slope predicted by multi-scale coupled finite element method Figure 15 shows the distribution of rotational displacement of rock blocks in the soil-rock mixture slope. The scale in the figure is the rotational displacement, and the unit is radian.
It can be seen from Figure 15 that the rotational displacement of rock blocks in the soil-rock mixture slope is mainly concentrated within the shear band, wherein the maximum rotational displacement of rock blocks is 40.7°, and the rotational displacement of particles outside the shear band is about 0°. The rotational displacement of the rock blocks in the shear band of the soil-rock mixture slope changes the mesostructures of the soil-rock mixture, which produces unrecoverable plastic deformation (Ren and Zhao 2021). In this process, the directional alignment of the rock blocks caused by their rotational displacement destroys the occlusion effect between rock blocks, and weakens the ability of rock blocks to transfer shear deformation to the surrounding soil matrix, thus leading to the plastic deformation being restrained and the occurrence of plastic strain localization. The plastic strain localization gradually evolves into a shear band, and the gradual penetration of the shear band causes the sliding failure of the soil-rock mixture slope. Therefore, the plastic strain concentration of the soil-rock mixture caused by the

FIGURE 14
The p~s curves simulated by multiscale finite element method.

Frontiers in Materials
frontiersin.org rotation and directional arrangement of the rock blocks in the soil-rock mixture slope is the most critical factor for the formation and development of the shear band of the soil-rock mixture slope.

Conclusion
According to the physical and mechanical effects generated by the interaction between soil matrix and rock blocks in the soil-rock mixture, the multiscale coupled elastic-plastic constitutive relationship of the soil-rock mixture is derived and the corresponding multiscale coupled finite element procedure is independently programmed. The cross-scale deformation characteristics of the soil-rock mixture slope are successfully simulated and predicted by the proposed multiscale finite element method. The main conclusions are as follows: (1) The multiscale coupled constitutive relationship established based on the micro motion characteristics of soil-rock mixture is embedded with the rotational freedom of rock blocks, which can effectively simulate and predict the crossscale shear strength of the soil-rock mixture.
(2) The multiscale finite element program can quantitatively simulate the rotational displacement of the element and effectively realize the cross-scale coupling relationship between the micro rotation of the rock blocks and the macro mechanical response of the soilrock mixture. Moreover, the intrinsic length scale that regularizes the stiffness matrix of the finite element method is embedded in the multiscale finite element program, which can overcome the pathological mesh-dependency of the classical finite element method, and thus objectively and completely describe the progressive deformation of soil-rock mixture slope from particle rotation at the microscale to the sliding failure of the soil-rock mixture at the macroscale. (3) The maximum rotational displacement of rock blocks inside the soil-rock mixture slope is 40.7°, mainly concentrated in the shear band, and the rotational displacement of rock blocks outside the shear zone is about 0°. The plastic strain concentration of the soil-rock mixture caused by the rotation and directional arrangement of the rock blocks in the soil-rock mixture is the most critical factor for the formation and development of the shear band of the soil-rock mixture slope.
Finally, it is noted that the proposed soil-rock cell element model is a preliminary investigation of the derivation of multiscale coupling theoretical model of the S-RM, only three groups of the in situ large-scale direct shear tests were carried out to calculate the parameters of the proposed multiscale soil-rock cell element model and evaluate the feasibility of the proposed model. Meanwhile, the shape and breakage of rock blocks have not been considered in the present work. Therefore, specific and systematic experimental studies and theoretical enhancements are expected to improve and validate the proposed multiscale FEM model in the future work.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions
JL: Conceptualization, Methodology, Software, Test. DF: Data curation, Writing-Original draft preparation Writing-Reviewing and Editing.