## 3 Mesoscale Simulation of Soil Rock Mixtures
### Implementation of the Mesoscale Computational Approach  
Based on the generated random SRM structure, the mesoscale computational approach can be performed by using the coupling calculation of COMSOL Multiphysics and Matlab. In this study, the SRM mesostructure generation, material discretization, material degradation, and material distribution matrix updates were performed by Matlab. The establishment of the initial condition and boundary condition, constitutive model, and calculation of stress and strain was performed by COMSOL. The coupling computation of Matlab and COMSOL was realized through the COMSOL livelink for Matlab. The communication between the two softwares was achieved through a generated txt file.  
Figure 4: Example of the area discretization.  
The general flowchart of the mesoscale computational approach is shown in Figure 6. To realize the initial parameters, a set of values of the gradation curve, volume content with different particle shapes, and fractal dimension is assigned to generate a random SRM structure. At the same time, the material discretization resolution and element number \(N\) of the whole simulation area is determined. The damage criterion and degradation rule are predefined as well in Matlab. In this study, the maximum tensile stress criterion (Equation (3)) and the Mohr-Coulomb criterion (Equation (4)) were chosen to determine whether the element is damaged.  
\[F_{1}=-\sigma_{3}-\sigma_{t0}=0 \tag{3}\]  
\[F_{2}=\sigma_{1}-\sigma_{3}\frac{1+\sin\varphi}{1-\sin\varphi}-\sigma_{c0}=0 \tag{4}\]  
where \(\sigma_{1}\) and \(\sigma_{3}\) are the maximum and minimum principal stress, \(\sigma_{t0}\) and \(\sigma_{c0}\) are the uniaxial tensile strength and uniaxial compressive strength, and \(\varphi\) is the internal friction angle.  
The constitutive model, initial condition, and boundary condition were selected in COMSOL. To save computational resources and improve computational efficiency, a linear elastic model was used in this study. At the same time, the geometry model and computational mesh were established in the initialization process.  
Then, a random SRM mesostructure was generated by the algorithm mentioned above. After that, the whole computational area was discretized into numbers of cubic areas, and a material distribution matrix was generated based on collision detection.  
Based on the material distribution matrix, geometry model, initial condition, and boundary condition, a loop calculation was performed. For each loop, the maximum and minimum principal stress, \(\sigma_{1}\) and \(\sigma_{3}\), and maximum and minimum principal strain, \(\epsilon_{1}\) and \(\epsilon_{3}\), were calculated and outputted as txt files by COMSOL. Then, the damage criterion was used to locate the damaged element. All damaged rock elements were transformed into soil elements, and the degradation law was used to reduce the modulus and strength of the damaged element. The material distribution matrix was updated after the degradation process. At the end of each loop computation, the calculation results were checked to  
Figure 5: Examples of SRM geometries with different rock content before and after discretization.  
determine whether the precision requirement was satisfied. If the requirements were met, the loop calculation would end. Otherwise, a new time step was added, and the loop calculation was performed again.  
However, it should be noted that, due to the inherent limitation of the proposed method, this method cannot simulate the contact behavior of soil matrix and rock particles, which might induce errors during the simulation.