## 1 Introduction
They found that for the development of the overall shear strength of the joint, the cohesive strength of the joint played an  
Figure 1: Schematic diagram for a gravity-type anchorage in a ground-anchored suspension bridge. important role when the applied normal stress confinement was fairly low (0.2\(\sim\)0.5 MPa), while the friction was more important under higher normal stress confinement. Based on experimental observation and theoretical analysis, several constitutive models have also been developed. [PERSON] and [PERSON] [30] developed a constitutive model to predict the mechanism of the development of overall shear strength of the C-R joint. This constitutive model simulated important physical processes including asperity sliding, asperity shearing, and post-peak behavior of the concrete-rock joint. However, [PERSON] and [PERSON] [31] have pointed out the limitations of existing empirical and analytical approaches. The obtained results posted significant dependencies on specific experimental site and test conditions for the existing empirical models, while the requirements of complex input parameters and oversimplifications showed difficulty in terms of the applicability of existing analytical approaches. As a result, they proposed a new theoretical model that can describe the main mechanism and requires joint roughness, basic rock properties, and boundary conditions for the cases of R-R and C-R joints. However, these constitutive models still could not have investigated the process of asperity degradation and shear fracturing during a direct shear test, indicating that the empirical-based and theoretical-based approaches suffer from significant limitations in their applications. The viability of continuum and discontinuum-based numerical methods to numerically simulate the shear behavior of R-R and C-R joints has advanced alongside computational geomechanics technology. [PERSON] et al. [32] used a three-dimensional (3D) FEM to develop a cohesive interface model of cemented C-R joints and simulated a direct shear test with a smooth joint surface. Using the discrete element method (DEM), [PERSON] et al. [33] simulated the shear behavior of rock joints using the bonded particle model (BPM), which is available in the two-dimensional (2D) particle flow code (PFC2D). They modeled the shear behavior of a saw-tooth triangular bonded and unbonded R-R joint and considered the joint roughness coefficient (JRC) profiles [7] as they investigated the rock damage process during a direct shear test. Ultimately, they successfully simulate the propagation of tensile cracks and the asperity crushing. Using the BPM of the 3D particle flow code (PFC3D), [PERSON] and [PERSON] [34] modeled a series of direct shear tests with rock-rock joints and investigated the effect of the geometrical topology and the micro-properties of the joints on shear behavior. Although the DEM with the BPM successfully modeled the fracture asperity and gouges, this method required a more tedious calibration process to acquire the unknown micro-parameters until the response of the macro-scale model accorded with the laboratory test results rather than the FEM application. Many of the computational resources also accompany simulating the crack propagation of the rock and concrete than the FEM. Recently, the combined finite-discrete element method (FDEM), originally proposed by [PERSON] [35; 36], has attracted significant attention in the field of computational solid mechanics. Since then, various FDEM codes have been developed in several research groups across the world [37; 38; 39; 40; 41; 42; 43; 44; 45; 46; 47; 48; 49; 50; 51; 52; 53; 54]. The FDEM can model the three important aspects of solid fracture, i.e., the continuous deformation in the intact regime, the transition from continuum to discontinuum, and the very complex contact mechanics between solid surfaces, including newly created macroscopic fracture surfaces.