# Improved Block Element Method for Simulating Rock Failure
## 4 Generation of Random Voronoi Grids  
The fracture path of rock masses is primarily governed by weak structural planes and the microstructural characteristics of intact rock blocks. In the BEM, weak planes such as joints and bedding are typically simulated by pre-defined layer elements, while intact rock domains are discretized using regular grids combined with virtual joints. But, within the theoretical framework of the traditional BEM, damage only occurs along the layer elements between blocks, and the blocks are rigid. This inherent limitation causes the geometric configuration of regular grids to artificially constrain the crack propagation paths. This condition may result in significant deviations between the simulated failure patterns and the natural failure processes of actual rock masses, which then affects the accuracy and reliability of the computational results. To overcome this limitation, a discretization strategy based on the Voronoi diagram is proposed. As a spatial proximity-based stochastic discretization method, the Voronoi diagram has been widely applied in other discrete element analyses [7; 8; 9; 10; 11; 33; 34; 35; 36]. It partitions the computational domain into a set of regions according to a series of seed points, with each region containing all the locations closest to a given seed point [37]. Studies have shown that Voronoi polygons closely resemble the geometric patterns of rock crystal structures and can naturally simulate complex fracture paths, such as intergranular cracking. Moreover, the randomness of the Voronoi tessellation overcomes the limitationsof pre-defined crack paths and allows fractures to develop autonomously based on the actual stress state. This dual characteristic enables a more realistic simulation of intricate fracture patterns such as branching and tortuous crack propagation during rock failure. However, it has been found that purely random Voronoi grids can lead to excessively distorted elements, resulting in numerical instability and inaccuracy [38]. To address this, the present study introduces the area-perimeter fractal dimension as an iterative control index. An iterative algorithm is applied to optimize the initial random mesh, improving its geometric quality while preserving sufficient randomness to ensure computational stability. Area and perimeter are the basic parameters used to describe the geometric shapes. For elements with the same area, the perimeter can vary significantly due to differences in shape complexity, because it reflects the influence of shape irregularity on boundary characteristics. It is worth noting that in fractal geometry, area (_Area_) and perimeter (_Peri_) exhibit a linear relationship in a logarithmic coordinate system [39], as expressed in Equation (30). It can be seen the area-perimeter fractal dimension (_Dim_) is twice the slope of the line formed by the perimeter and area on a double-logarithmic coordinate system. When _Dim_ = 1, it corresponds to regular Euclidean geometric shapes, for which the perimeter is proportional to the square root of the area. \[\text{Lg}(Peri)=\frac{Dim}{2}\text{Lg}(Area)+\text{const.} \tag{30}\]  
As shown in Figure 5, the random Voronoi grid generation process used in this study is as follows:  
* Step 1: Random seed points are generated within and around the computational domain, as shown in Figure 5a;
* Step 2: The Voronoi diagram is generated using the seed points from the previous step within the computational domain, as shown in Figure 5b;
* Step 3: The Lloyd iteration algorithm [40] is applied for mesh optimization, in which the area (Area) and perimeter (Peri) of each element are extracted. The Lg(Area)-Lg(Peri) data is then subjected to the least squares linear regression method to calculate the fractal dimension Dim. If Dim stabilizes, the iteration process terminates; otherwise, the process returns to Step 2. A homemade MATLAB r2023a program is used to discretize a 50 mm \(\times\) 100 mm rectangular specimen.