## 1 Introduction  
Natural defects such as joints, fissures and holes, are commonly included in rock masses because of weathering evolution, geological structure and human actions [1]. Under loading, cracks may initiate from the corners or tips of these defects due to stress concentration, and crack propagation and coalescence may eventually lead to instability and failure of rock structures [2]. Moreover, the excavated openings in geotechnical engineering, such as tunnels, roadways and shafts, can also be regarded as artificial defects formed in rock masses. Nowadays, it is generally believed that the mechanical properties of rocks or rock masses are significantly influenced by the size, number, shape and distribution of these defects [3; 4; 5; 6; 7; 8; 9]. Virtually, the defects in natural rock are usually filled with sands, ooze, calcareous cements, etc. The filling can also be manually filled with mortar, concrete and other cementing materials during the construction process, such as grouting, shotcrete, etc., to arrest the crack propagation and improve the strength and stability of rock mass [10; 11]. It is well known that the rock strength is governed by the crack propagation during loading [12; 13], and the crack behavior around filled defects is complicated. Therefore, it is of great engineering significance to study the effect of filling on mechanical behavior and fracture evolution of rocks. Until now, some studies on the cracking behavior of rock or rock-like materials containing various inclusions have been conducted. Based on [PERSON]'s solution [14], [PERSON] [15] first derived a closed form solution to calculate the stress concentration around a spherical particle in the matrix under compressive stress. [PERSON] and [PERSON] [16] applied the sliding crack model to study the crack propagation and coalescence around polygonal inclusions in a heterogeneous matrix subjected to far field compression. The results showed that the crack initiated from the inclined interface will continue to slip along the interface in the form of shear, but will propagate into the matrix in the form of tension when it reaches the ends of the interface. Compression tests on prismatic concrete with two circular openings and inclusions were first carried out by [PERSON] and [PERSON] [17]. It was observed that interface cracking always initiated at the interface of inclusions, while matrix cracking always initiated at the roof and floor of the circular openings. [PERSON] and [PERSON] [18] performed uniaxial compression tests concrete specimens containing stone-matrix interface oriented at four different angles, and studied the mixed-mode crack initiation initiated from the aggregate-matrix interface in concrete using fracture mechanics approach, finite element modeling and holographic interferometry techniques. The analytical results were later verified by [PERSON] et al. [19]. They found that tensile debonding at the interface of inclusion always precedes the initiation of tensile cracks in the matrix. Moreover, the mode of crack initiated at the interface of rectangular inclusion was not pure tensile, but a mixed-mode (tensile behavior dominated, however). [PERSON] et al. [20] and [PERSON] et al. [21] compared the crack behavior of unfilled and filled flaws by conducting uniaxial compression tests on rock-like specimens, and found that the crack initiation time, initiation location and propagation behavior of cracks were all influenced by the filling. Recently, the effect of different filling materials on rock properties has been the focus of many researchers. [PERSON] and [PERSON] [22] studied the crack process in prismatic gypsum specimens containing either one or two inclusions with different sizes, shapes, stiffness and strength, and the experimental tests were further discussed and verified by [PERSON] and [PERSON] [23] through numerical simulation with the enriched numerical manifold method.