# Improved Block Element Method for Simulating Rock Failure
## 1 Introduction  
As discontinuous deformation methods [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17], the rigid-body spring method (RBSM) [7, 8, 9, 10, 11], interface element method (IEM) [12, 13, 14], and block element method (BEM) [15, 16, 17] all use the rigid-body displacements of block centroids as the primary unknown variables. Based on equilibrium conditions, compatibility conditions, and the constitutive relations of layer elements between adjacent blocks, these methods establish an implicit solution framework for structural analysis. The main difference lies in how block interactions are represented: the RBSM models the interactions as \"point contacts\" by a set of springs located at the center of the interface, connecting adjacent blocks. In contrast, the BEM and IEM use a \"surface contact\" model instead of \"point contacts\", capturing the mechanical behavior of materials more realistically through layer element constitutive relations and better reflecting the non-uniform stress distribution and deformation between blocks. Unlike the finite element method (FEM), the BEM does not require compatibility between blocks, allowing more flexible block geometries. This feature enables the BEM to discretize blocks of arbitrary shapes, making it suitable for complex rock masses [15, 17, 18]. Due to its unique advantages in handling discontinuous deformation problems, the BEM has drawn wide attention. [PERSON] et al. [15; 16] developed the fundamental theory and methodology of the BEM and applied it to the stability analysis of jointed rock masses. [PERSON] et al. [19] extended the BEM to elastoplastic and contact analysis of rock masses using parametric variational principles. With the development of the BEM, researchers have proposed various improvements to address its theoretical and application limitations. Regarding displacement modes, to overcome the limitations of the rigid-body displacement assumption, several new formulations have been developed, including the constant strain mode [20], the curvature mode [21], and high-order polynomial modes [22; 23], all of which have significantly improved numerical accuracy. To broaden its application scope, coupling strategies [17; 24; 25; 26], especially FEM-BEM hybrid methods, have also been explored. For instance, [PERSON] et al. [25] combined the FEM and BEM to enhance their ability to solve complex engineering problems. However, the BEM still has limitations in mesh discretization and constitutive modeling. In BEM simulations, the rigid blocks cannot fracture internally; failure is restricted to the layer elements between blocks. As a result, the initiation and propagation of cracks are highly dependent on the initial discretization. Presently, a dual strategy is typically adopted [15]. For jointed rock masses, real structural planes are used for meshing; for intact rock, artificial joints are introduced, often with triangular or regular meshes. While this approach captures failure along weak planes in jointed rocks, it imposes constraints on crack propagation directions in intact rocks, reducing simulation reliability. In rock mechanics, various failure criteria have been developed to characterize the complex failure behavior of rock materials under different loading conditions, such as the Mohr-Coulomb criterion [15; 17; 19; 20] and Hoek-Brown criterion [27; 28]. The existing BEM mainly uses the Mohr-Coulomb criterion for nonlinear analysis, describing the failure process of layer elements through a simple stress release mechanism. While this simplified model offers high computational efficiency, it fails to accurately capture the entire failure process, especially the nonlinear softening and damage accumulation behavior. These limit the method's applicability to complex fracture simulations. To overcome the limitations of the traditional BEM, this study introduces improvements in both the constitutive model and mesh discretization.