# Improved Block Element Method for Simulating Rock Failure
## 2 Fundamental Equations of the Block Element Method  
In the modeling of the block element method (BEM), the computational domain is discretized into several rigid block elements, which interact with each other through thin-layer elements [15; 16], as illustrated in Figure 1. These layer elements can represent actual geological discontinuities (such as joints and fractures) or be defined as virtual contact interfaces to meet specific modeling requirements. Inspired by the finite element analysis framework, the shape functions are employed to construct the rigid-body displacement mode in the BEM. For two-dimensional problems, the displacement field at any point on the block boundary can be expressed based on coordinate transformation theory by assuming the block undergoes only small rotations. Taking block \(i\) as an example (see Figure 2), the displacement field at a boundary point is given as follows:  
\[\mathbf{U}= \begin{bmatrix}u&v\end{bmatrix} ^{T}=\mathbf{N}_{i}\mathbf{U}_{i} \tag{1}\]\[\mathbf{N}_{i}=\left[ \begin{array}{ccc}1&0&y_{gi}-y\\ 0&1&x-x_{gi}\end{array} \right] \tag{2}\]  
\[\mathbf{U}_{i}=\left[ \begin{array}{ccc}u_{gi}&v_{gi}&\theta_{gi}\end{array} \right] ^{T} \tag{3}\]  
where \(u\) and \(v\) denote the translational displacements at an arbitrary point on the block element. \(\mathbf{N}_{i}\) represents the shape function matrix, and \(\mathbf{U}_{i}\) is the nodal displacement vector at the centroid of block \(i\). The \(u_{gi}\), \(v_{gi}\), and \(\theta_{gi}\) components correspond to the translational and rotational displacements of the centroid of the block, respectively. \(x_{gi}\) and \(y_{gi}\) are the global coordinates of the block centroid, while \(x\) and \(y\) are the global coordinates of an arbitrary point on the block. Given that the predefined thickness \(t\) of the thin-layer element is very small, the problem can be considered as a plane strain problem in the local coordinate system. Therefore, based on the geometric relations in elastic mechanics, the relative displacement and stress of the layer element surface can be calculated as follows:  
\[\delta=\left[ \begin{array}{cc}\delta_{n}&\delta_{s}\end{array} \right]^{T}= \mathbf{L}(\mathbf{U}_{j}-\mathbf{U}_{i}) \tag{4}\]  
\[\mathbf{\epsilon}=\left[ \begin{array}{cc}\varepsilon_{n}&\varepsilon_{s}\end{array}  \right]^{T}=\frac{\delta}{t} \tag{5}\]  
\[\mathbf{\sigma}=\mathbf{D}^{\star}\mathbf{\epsilon}=\mathbf{D}\delta \tag{6}\]  
\[\mathbf{L}=\left[ \begin{array}{cc}l&m\\ -m&l\end{array} \right] \tag{7}\]  
where \(\delta_{n}\), \(\delta_{s}\) and \(\varepsilon_{n}\), \(\varepsilon_{s}\) are the relative displacement and strain of the layer element surface in the local coordinate system, respectively. \(\mathbf{L}\) denotes the coordinate transformation matrix  
Figure 1: Schematic diagram of computational domain discretization for BEM. Figure 2: Relative displacement between adjacent blocks. between the global coordinate system and the local coordinate system, while \(l\) and \(m\) are the components. For isotropic elastic materials, the constitutive relationship in the local coordinate system can be expressed as follows:  
\[\mathbf{D}^{\star}= \begin{bmatrix}\frac{E}{1-\mu^{2}}&0\\ 0&\frac{E}{2(1+\mu)}\end{bmatrix}  \tag{8}\]  
\[\mathbf{D}= \begin{bmatrix}k_{n}&0\\ 0&k_{s}\end{bmatrix} = \begin{bmatrix}\frac{E}{2t(1-\mu^{2})}&0\\ 0&\frac{E}{2t(1+\mu)}\end{bmatrix}  \tag{9}\]  
where \(E\) is the Young modulus of the material, \(\mu\) is Poisson's ratio, and \(k_{n}\) and \(k_{s}\) are, respectively, the stiffness coefficients of the normal and tangential direction of the local coordinate system. Under the initial assumption, the blocks are considered rigid and cannot store strain energy. All strain energy is accumulated within the layer elements.