Generation of Numerical Models of Anisotropic Columnar Jointed Rock Mass Using Modified Centroidal Voronoi Diagrams

A columnar joint network is a natural fracture pattern with high symmetry, which leads to the anisotropy mechanical property of columnar basalt. For a better understanding the mechanical behavior, a novel modeling method for columnar jointed rock mass through field investigation is proposed in this paper. Natural columnar jointed networks lies between random and centroidal Voronoi tessellations. This heterogeneity of columnar cells in shape and area can be represented using the coefficient of variation, which can be easily estimated. Using the bisection method, a modified Lloyd’s algorithm is proposed to generate a Voronoi diagram with a specified coefficient of variation. Modelling of the columnar jointed rock mass using six parameters is then presented. A case study of columnar basalt at Baihetan Dam is performed to demonstrate the feasibility of this method. The results show that this method is applicable in the modeling of columnar jointed rock mass as well as similar polycrystalline materials.


Introduction
Columnar jointed rock is a typical geological structure formed by a lot of ordered colonnades like the world famous natural wonders of Fingal's Cave and Giant's Causeway [1,2].Being a miraculous natural phenomenon, the study on the geological structure of columnar jointed rock can be dated back to the 17th century [3].Nowadays, the formation of columnar jointing is reasonably understood as a result of cracks propagation in cooling lava flows [4][5][6][7].In situ observations and laboratory tests have also been employed in the study of columnar jointing [8][9][10].
With the increase of human activities, some engineering projects have encountered columnar jointing, and a large hydropower station is even founded on it [11].Because of its adverse geologic conditions, study on the mechanical properties of columnar jointed rock mass is very important in engineering projects.King et al. studied elastic-wave propagation in columnar joints with a series of cross-hole acoustic measurements made between four horizontal boreholes drilled from a near-surface underground opening situated in a basaltic rock mass [12].In connection with the Baihetan Hydropower Station project, several studies have been implemented for the mechanical property of columnar jointed rock.Meng systematically analyzed the anisotropic properties of columnar jointed basalt using analytical and numerical methods [13].Zheng et al. proposed a 3-D modeling method for columnar basalt using random Voronoi tessellation [14].In References [15][16][17], 2-D and 3-D discrete element simulation methods of columnar are implemented, and the representative elementary volume scale and equivalent mechanical parameters are obtained.Based on meso-structural analysis, macro-anisotropic constitutive model is developed and 3-D numerical simulation is performed for the dam foundation considering the effect of columnar joints [18].Laboratory tests are carried out using experimental analogs of jointing under uniaxial compression [19].These research results indicate that the characteristics of columnar jointed rock mass are very complex with anisotropic and nonlinear properties.
Most of the modelling methods of columnar jointed rock mass use the program developed by Zheng et al. [14]; however, this method has some shortcomings, such as a low efficiency and immutable columnar shapes, as pointed out in Reference [20].The objective of this paper is to overcome these shortcomings in the generation of columnar joints in rock mass.An efficient and controllable method for the generation of columnar joints is proposed based on a constrained centroidal Voronoi tessellation diagram.In Section 2, the characteristics of typical columnar rock mass are briefly introduced.The properties of Voronoi diagrams and the latest related research on this topic are discussed in Section 3. Based on these discussion and algorithms, a detailed procedure for the generation of columnar jointed rock mass is developed in Section 4. A case study on columnar jointed basalt generation for the Baihetan Hydropower Station is presented in Section 5. Finally, the conclusion is given in Section 6.

Typical Shapes of Columnar Jointed Rock Mass
With the development of hydropower station, columnar jointing is often encountered in southwestern China.As a result of cooling lava flows and ash-flow tuffs, it occurs in many types of volcanic rocks in the formation of a regular array of polygonal prisms or columns [21].Several famous typical columnar joints are shown in Figure 1.
Hydropower Station project, several studies have been implemented for the mechanical property of columnar jointed rock.Meng systematically analyzed the anisotropic properties of columnar jointed basalt using analytical and numerical methods [13].Zheng et al. proposed a 3-D modeling method for columnar basalt using random Voronoi tessellation [14].In References [15][16][17], 2-D and 3-D discrete element simulation methods of columnar are implemented, and the representative elementary volume scale and equivalent mechanical parameters are obtained.Based on mesostructural analysis, macro-anisotropic constitutive model is developed and 3-D numerical simulation is performed for the dam foundation considering the effect of columnar joints [18].Laboratory tests are carried out using experimental analogs of jointing under uniaxial compression [19].These research results indicate that the characteristics of columnar jointed rock mass are very complex with anisotropic and nonlinear properties.
Most of the modelling methods of columnar jointed rock mass use the program developed by Zheng et al. [14]; however, this method has some shortcomings, such as a low efficiency and immutable columnar shapes, as pointed out in Reference [20].The objective of this paper is to overcome these shortcomings in the generation of columnar joints in rock mass.An efficient and controllable method for the generation of columnar joints is proposed based on a constrained centroidal Voronoi tessellation diagram.In Section 2, the characteristics of typical columnar rock mass are briefly introduced.The properties of Voronoi diagrams and the latest related research on this topic are discussed in Section 3. Based on these discussion and algorithms, a detailed procedure for the generation of columnar jointed rock mass is developed in Section 4. A case study on columnar jointed basalt generation for the Baihetan Hydropower Station is presented in Section 5. Finally, the conclusion is given in Section 6.

Typical Shapes of Columnar Jointed Rock Mass
With the development of hydropower station, columnar jointing is often encountered in southwestern China.As a result of cooling lava flows and ash-flow tuffs, it occurs in many types of volcanic rocks in the formation of a regular array of polygonal prisms or columns [21].Several famous typical columnar joints are shown in Figure 1.Although the characteristics of columnar jointed rock is typical, they vary for different locations.The column in some areas is regular and straight like Giant's Causeway, while some other places like Baihetan Dam, it is irregular.Besides, the diameters of columns vary from meters to centimeters.Columns are usually parallel and straight, and the length can be as much as several meters.Furthermore, the number of sides of an individual column also varies from three to eight.In order to give a more specific illustration of the columns, a colorized O'Reilly map is presented in Figure 2. In this figure, hexagons are green, heptagons are blue, octagons are purple, pentagons are orange, and Although the characteristics of columnar jointed rock is typical, they vary for different locations.The column in some areas is regular and straight like Giant's Causeway, while some other places like Baihetan Dam, it is irregular.Besides, the diameters of columns vary from meters to centimeters.Columns are usually parallel and straight, and the length can be as much as several meters.Furthermore, the number of sides of an individual column also varies from three to eight.In order to give a more specific illustration of the columns, a colorized O'Reilly map is presented in Figure 2.
In this figure, hexagons are green, heptagons are blue, octagons are purple, pentagons are orange, and squares are red.From the introduction of columnar jointed rock mass above, a diagrammatic drawing is proposed in Figure 3.It is appropriate to model columnar joints using Voronoi diagrams.

Classical Voronoi Tessellation
A Voronoi diagram is a classical domain partitioning method named after Georgy Voronoi, a Russian mathematician.Voronoi generation is the dual algorithm of Delaunay triangulation.A typical Voronoi figure is generated using the perpendicular bisector of the lines composed by a set of points called seeds.The Voronoi diagram is a topic in computational geometry and has already widely applied in some other research areas like hydrology and crystal mechanics.
As shown in Figure 4, a Voronoi diagram with 10 random seeds Pi is generated.It can be seen that the domain is divided into 10 patches and each patch has a single seed.For the generation of a

Classical Voronoi Tessellation
A Voronoi diagram is a classical domain partitioning method named after Georgy Voronoi, a Russian mathematician.Voronoi generation is the dual algorithm of Delaunay triangulation.A typical Voronoi figure is generated using the perpendicular bisector of the lines composed by a set of points called seeds.The Voronoi diagram is a topic in computational geometry and has already widely applied in some other research areas like hydrology and crystal mechanics.
As shown in Figure 4, a Voronoi diagram with 10 random seeds Pi is generated.It can be seen that the domain is divided into 10 patches and each patch has a single seed.For the generation of a

Classical Voronoi Tessellation
A Voronoi diagram is a classical domain partitioning method named after Georgy Voronoi, a Russian mathematician.Voronoi generation is the dual algorithm of Delaunay triangulation.A typical Voronoi figure is generated using the perpendicular bisector of the lines composed by a set of points called seeds.The Voronoi diagram is a topic in computational geometry and has already widely applied in some other research areas like hydrology and crystal mechanics.
As shown in Figure 4, a Voronoi diagram with 10 random seeds P i is generated.It can be seen that the domain is divided into 10 patches and each patch has a single seed.For the generation of a Voronoi diagram, a Delaunay triangulation is first implemented.Then the perpendicular bisector is plotted to partition the domain into different Voronoi cells.With the development of computational geometry, Voronoi tessellation is included in many software and packages like MATLAB, Mathematica, and SciPy.However, there are two obvious shortcomings of the classical Voronoi diagram: (a) The Voronoi cell is not closed.For the Voronoi diagram in Figure 4, only one cell is closed and the other nine cells are open.This brings inconvenience for the analysis.(b) The shape of the Voronoi cell is random and it is very hard to generate a Voronoi diagram with a specified statistical distribution.

Constrained Voronoi Diagram Generation
To overcome the first drawback of classical Voronoi tessellation, a constrained Voronoi diagram algorithm is employed to generate a closed Voronoi tessellation with seed points P i and domain D. Taking a model with 16 points (x i , y i ) and a square domain in Figure 5a as example, this algorithm can be described as follows [22]: (a) A Voronoi diagram is generated using a classical tessellation method (Figure 5b).

Constrained Voronoi Diagram Generation
To overcome the first drawback of classical Voronoi tessellation, a constrained Voronoi diagram algorithm is employed to generate a closed Voronoi tessellation with seed points Pi and domain D. Taking a model with 16 points (xi, yi) and a square domain in Figure 5a as example, this algorithm can be described as follows [22]: (a) A Voronoi diagram is generated using a classical tessellation method (Figure 5b).

Constrained Voronoi Diagram Generation
To overcome the first drawback of classical Voronoi tessellation, a constrained Voronoi diagram algorithm is employed to generate a closed Voronoi tessellation with seed points Pi and domain D. Taking a model with 16 points (xi, yi) and a square domain in Figure 5a as example, this algorithm can be described as follows [22]: (a) A Voronoi diagram is generated using a classical tessellation method (Figure 5b).

Random and Centroidal Voronoi Diagram
The Voronoi diagram created using a constrained Voronoi tessellation method with randomly generated seeds is shown in Figure 6a.It can be seen that the Voronoi cell is irregular and there is a clear deviation between the Voronoi seeds and the centroids.In order to obtain more regular columnar as shown in Figure 2, an algorithm for centroidal Voronoi diagrams (Figure 6b) is proposed to tackle this problem.Compared with the random Voronoi diagram, the seed of a cell is coincident with its centroid in centroidal Voronoi tessellation (CVT).It is widely used in related fields like mesh generation and data compression.Previous study has indicated that some natural patterns like Giant's Causeway can be represented by the centroidal Voronoi tessellation [5].
There are a number of methods for the generation of centroidal Voronoi diagram, among which are the algorithms by Lloyd and by MacQueen [23,24], which are widely used, and some other methods are variations of these two methods.Considering that MacQueen's algorithm needs twice as many Monte Carlo simulations in each iteration, which consumes a large amount computation time, Lloyd's centroidal method is employed for modelling columnar jointed rock mass in this study.

Random and Centroidal Voronoi Diagram
The Voronoi diagram created using a constrained Voronoi tessellation method with randomly generated seeds is shown in Figure 6a.It can be seen that the Voronoi cell is irregular and there is a clear deviation between the Voronoi seeds and the centroids.In order to obtain more regular columnar as shown in Figure 2, an algorithm for centroidal Voronoi diagrams (Figure 6b) is proposed to tackle this problem.

Random and Centroidal Voronoi Diagram
The Voronoi diagram created using a constrained Voronoi tessellation method with randomly generated seeds is shown in Figure 6a.It can be seen that the Voronoi cell is irregular and there is a clear deviation between the Voronoi seeds and the centroids.In order to obtain more regular columnar as shown in Figure 2, an algorithm for centroidal Voronoi diagrams (Figure 6b) is proposed to tackle this problem.Compared with the random Voronoi diagram, the seed of a cell is coincident with its centroid in centroidal Voronoi tessellation (CVT).It is widely used in related fields like mesh generation and data compression.Previous study has indicated that some natural patterns like Giant's Causeway can be represented by the centroidal Voronoi tessellation [5].
There are a number of methods for the generation of centroidal Voronoi diagram, among which are the algorithms by Lloyd and by MacQueen [23,24], which are widely used, and some other methods are variations of these two methods.Considering that MacQueen's algorithm needs twice as many Monte Carlo simulations in each iteration, which consumes a large amount computation time, Lloyd's centroidal method is employed for modelling columnar jointed rock mass in this study.Compared with the random Voronoi diagram, the seed of a cell is coincident with its centroid in centroidal Voronoi tessellation (CVT).It is widely used in related fields like mesh generation and data compression.Previous study has indicated that some natural patterns like Giant's Causeway can be represented by the centroidal Voronoi tessellation [5].
There are a number of methods for the generation of centroidal Voronoi diagram, among which are the algorithms by Lloyd and by MacQueen [23,24], which are widely used, and some other methods are variations of these two methods.Considering that MacQueen's algorithm needs twice as many Monte Carlo simulations in each iteration, which consumes a large amount computation time, Lloyd's centroidal method is employed for modelling columnar jointed rock mass in this study.

Lloyd's Algorithm
In this work, Lloyd's algorithm, named after Stuart P. Lloyd, is employed.It is an iteration method to partition the domain into well-shaped and uniformly-sized convex cells [25].The convergence of Lloyd's algorithm to a centroidal Voronoi diagram has been proven.The random Voronoi diagram can be iterated to a centroidal Voronoi diagram using Lloyd's algorithm and the algorithm can be simplified as follows: (1) For an initial seeds y i , generate a Voronoi diagram using constrained Voronoi tessellation; (2) Compute the centroid z i of the Voronoi diagram of y i ; (3) Move the generating point y i to its centroid z i ; (4) Repeat Steps 1 to 3 until all generating points converge to the centroids.

Estimation of the Centroid
The centroidal Voronoi algorithm is very simple but the calculation of the centroid of each Voronoi cell is time-consuming.In each iteration of Lloyd's method, the centroidal positions of all shapes of the Voronoi diagram are calculated.Thus, the estimation of the centroid is the most time-consuming task.For determining the centroid of a Voronoi region, a simple formula is: where χ is the position, V i is the region area, and ρ(χ) is the density function with ρ(χ) = 1 being the default.
In order to calculate the position of centroid fast, two efficient methods for evaluating the integrals are used in this work.One is based on an integration algorithm on triangle partitions and the other is a sampling method.An integration algorithm can calculate the coordinates of centroids exactly but not as fast as a sampling method; whereas a sampling method is fast but the results are not as accurate as that of an integration algorithm.

(1) Integration method on triangle partitions
For an arbitrary polygon with n vertices, it can be divided into n − 2 triangles with a clockwise or anti-clockwise direction (Figure 7).For each triangle with vertex A i (x i , y i ) (i = 1, 2, 3), the coordinates of centroid are given by: The area of the triangle is calculated as: Considering that the polygon can be discretized into n − 2 triangles and the centroid and area of each triangle can be expressed as G i (x gi , y gi ) and S i , the coordinates of the centroid of the polygon can be obtained as: Symmetry 2018, 10, x FOR PEER REVIEW 7 of 15 (

2) Sampling method
Although the integration on triangle partitions is a good solution for the estimation of the centroid of a polygon, it takes a long time to complete each step in Lloyd's iteration.For a quick calculation, a simple and efficient sampling method is proposed as an alternative.
For a certain area with different Voronoi tessellations, a set of random points are distributed in this area (Figure 8).For each Voronoi cell, if there are n points Pi(xi, yi) in this area, then the centroid of this cell can be approximated as the average of these points: In this way, the determination of centroids changes to the partition of areas into different Voronoi cells, i.e., the "N-D nearest point search" problem.Luckily, the nearest point search problem has a mature mathematical solution [26] and this algorithm has been written in a number of software libraries, such as Python and MATLAB.Using these libraries, all the random points can be divided efficiently and the centroid of Voronoi diagram can be computed quickly.With more sampling points, the calculation of centroid will be more accurate.If the number of random points is selected properly, this algorithm can be both fast and relatively accurate.(

2) Sampling method
Although the integration on triangle partitions is a good solution for the estimation of the centroid of a polygon, it takes a long time to complete each step in Lloyd's iteration.For a quick calculation, a simple and efficient sampling method is proposed as an alternative.
For a certain area with different Voronoi tessellations, a set of random points are distributed in this area (Figure 8).For each Voronoi cell, if there are n points P i (x i , y i ) in this area, then the centroid of this cell can be approximated as the average of these points: In this way, the determination of centroids changes to the partition of areas into different Voronoi cells, i.e., the "N-D nearest point search" problem.Luckily, the nearest point search problem has a mature mathematical solution [26] and this algorithm has been written in a number of software libraries, such as Python and MATLAB.Using these libraries, all the random points can be divided efficiently and the centroid of Voronoi diagram can be computed quickly.With more sampling points, the calculation of centroid will be more accurate.If the number of random points is selected properly, this algorithm can be both fast and relatively accurate.

Numerical Implementation and Discussion
In the iterations of Lloyd's algorithm, the Voronoi generators and centroids are known, and energy can be employed to describe the Voronoi diagram [23].For the Voronoi tessellation { } 1 n i i= Ω , the energy is defined as: ( ) where χ is the position of the Voronoi generator, g χ represents the position of the centroid, and ( ) ρ χ is the density function with ( )=1 ρ χ being the default.
However, if the generator of Voronoi diagram is not known beforehand, energy cannot be estimated.For example, it is not easy to calculate the energy of Voronoi cells in Figure 2. Indirectly, the coefficent of variation, which has been proven as an effective parameter in the estimation of Voronoi diagram [27,28], is employed for the description of the properties of Voronoi cells.In the process of smoothing energy from a randomly-generated Voronoi diagram to a centroidal Voronoi diagram, the coefficient of variation value is recorded with the following formula: where SD is the standard deviation and m is the mean of Voronoi cell areas.The coefficient of variation gives a quantitative indication on the spatial distribution: the higher the coefficient of variation, the higher the tendency of cells to aggregate into clusters.
Based on the related algorithms described, numerical implementation is done in MATLAB, in which some good programming techniques from the open-source program of Burkardt are referenced [29].The 50 seeds for Voronoi diagram are generated randomly.After 50 iterations, the points are well-spaced and Lloyd's algorithm for computing centroidal Voronoi tessellation converged overall (Figure 9).In this way, the coefficient of variation has the same evolution trend as energy (E), and it is reasonable to use coefficient of variation to describe the distribution property of the Voronoi tessellation.

Numerical Implementation and Discussion
In the iterations of Lloyd's algorithm, the Voronoi generators and centroids are known, and energy can be employed to describe the Voronoi diagram [23].For the Voronoi tessellation {Ω i } n i=1 , the energy is defined as: where χ is the position of the Voronoi generator, χ g represents the position of the centroid, and ρ(χ) = 1 is the density function with ρ(χ) = 1 being the default.However, if the generator of Voronoi diagram is not known beforehand, energy cannot be estimated.For example, it is not easy to calculate the energy of Voronoi cells in Figure 2. Indirectly, the coefficent of variation, which has been proven as an effective parameter in the estimation of Voronoi diagram [27,28], is employed for the description of the properties of Voronoi cells.In the process of smoothing energy from a randomly-generated Voronoi diagram to a centroidal Voronoi diagram, the coefficient of variation value is recorded with the following formula: where SD is the standard deviation and m is the mean of Voronoi cell areas.The coefficient of variation gives a quantitative indication on the spatial distribution: the higher the coefficient of variation, the higher the tendency of cells to aggregate into clusters.
Based on the related algorithms described, numerical implementation is done in MATLAB, in which some good programming techniques from the open-source program of Burkardt are referenced [29].The 50 seeds for Voronoi diagram are generated randomly.After 50 iterations, the points are well-spaced and Lloyd's algorithm for computing centroidal Voronoi tessellation converged overall (Figure 9).In this way, the coefficient of variation has the same evolution trend as energy (E), and it is reasonable to use coefficient of variation to describe the distribution property of the Voronoi tessellation.

Modeling of Columnar Jointed Rock Mass
Analysis of typical shapes of columnar jointing shows that a columnar jointed rock mass can be modelled through the extrusion of a 2-D Voronoi diagram.Hence, the generation of 2-D columnar jointing, which is consistent with the site condition is a key task.In Section 3, it is shown that coefficient of variation can describe the distribution property in the iterations of Lloyd's algorithm.The centroidal Voronoi diagram is applicable in the description of a natural phenomenon such as Giant's Causeway.In fact, natural columnar jointing is not an exactly centroidal Voronoi diagram, and it lies between random and centroidal Voronoi diagrams.In the process of transforming a Voronoi diagram from totally random to centroidal, the value of CV indicates that the Voronoi cells change from heterogeneous to homogeneous.Therefore, it is feasible to describe the characteristics of columnar joints by employing two main parameters: the columnar density representing the scale and coefficient of variation reflecting the variation property.
A procedure for the generation of columnar jointed rock mass is proposed (Figure 10).Based on field investigation of geological conditions and site images, the columnar jointing at the site of interest is characterized by six parameters: columnar density (CD), coefficient of variation (CV), dip direction (DD), dip angle (DIP), transverse joint distance (TD), and probability (TP).The number of random seeds for generating a Voronoi diagram is estimated according to the calculated density.Based on the idea of the bisection method, a novel modified constrained centroidal Voronoi smooth algorithm is implemented to make the CV converge to the specified value as follows: (1) For a Voronoi diagram with CV larger than the specified value, calculate the centroid PC.
(2) The new generator Pnew,g is set at the midpoint of the old generator Pold,g and the centroid PC.
Calculate the coefficient of variation CV of the new Voronoi diagram.(3) Repeat Steps 1 and 2 to obtain the generator until the coefficient of variation CV converges to the specified value with a prescribed accuracy.

Modeling of Columnar Jointed Rock Mass
Analysis of typical shapes of columnar jointing shows that a columnar jointed rock mass can be modelled through the extrusion of a 2-D Voronoi diagram.Hence, the generation of 2-D columnar jointing, which is consistent with the site condition is a key task.In Section 3, it is shown that coefficient of variation can describe the distribution property in the iterations of Lloyd's algorithm.The centroidal Voronoi diagram is applicable in the description of a natural phenomenon such as Giant's Causeway.In fact, natural columnar jointing is not an exactly centroidal Voronoi diagram, and it lies between random and centroidal Voronoi diagrams.In the process of transforming a Voronoi diagram from totally random to centroidal, the value of CV indicates that the Voronoi cells change from heterogeneous to homogeneous.Therefore, it is feasible to describe the characteristics of columnar joints by employing two main parameters: the columnar density representing the scale and coefficient of variation reflecting the variation property.
A procedure for the generation of columnar jointed rock mass is proposed (Figure 10).Based on field investigation of geological conditions and site images, the columnar jointing at the site of interest is characterized by six parameters: columnar density (CD), coefficient of variation (CV), dip direction (DD), dip angle (DIP), transverse joint distance (TD), and probability (TP).The number of random seeds for generating a Voronoi diagram is estimated according to the calculated density.Based on the idea of the bisection method, a novel modified constrained centroidal Voronoi smooth algorithm is implemented to make the CV converge to the specified value as follows: (1) For a Voronoi diagram with CV larger than the specified value, calculate the centroid P C .
(2) The new generator P new,g is set at the midpoint of the old generator P old,g and the centroid P C .
Calculate the coefficient of variation CV of the new Voronoi diagram.(3) Repeat Steps 1 and 2 to obtain the generator until the coefficient of variation CV converges to the specified value with a prescribed accuracy.
Having obtained a 2-D Voronoi diagram with specified CD and CV, it can be extruded to a 3-D columnar shape in accordance with DD and DIP.Then for each extruding Voronoi cell and its length, a transverse joint can be generated by the determination of the intersection between the extruding Voronoi cell and the transverse plane.Sometimes, the transverse joint is non-penetrating, and the transverse joint can be selected with the specified TD and TP.In this way, columnar jointed rock mass with the site properties can be generated.Having obtained a 2-D Voronoi diagram with specified CD and CV, it can be extruded to a 3-D columnar shape in accordance with DD and DIP.Then for each extruding Voronoi cell and its length, a transverse joint can be generated by the determination of the intersection between the extruding Voronoi cell and the transverse plane.Sometimes, the transverse joint is non-penetrating, and the transverse joint can be selected with the specified TD and TP.In this way, columnar jointed rock mass with the site properties can be generated.

Engineering Geological Investigation
The Baihetan Hydropower Station is a multi-purpose project for harnessing the Yangtze River, developed mainly for power generation, flood control, sediment flushing, and improving the navigation conditions in the reservoir and downstream.The Baihetan Hydropower Station is located on the downstream reaches of the Jinsha River.At the dam site, the most apparent stratigraphic lithology is basalt that belongs to the Emeishan (Emei Mt.) formation of the Permain System (P2β).In the construction of the plant's hydraulic structures, columnar jointed basalt is found to be widely distributed at the arch dam foundation, underground caverns, abutment slopes, and other water conservancy tunnels (Figure 11).

Engineering Geological Investigation
The Baihetan Hydropower Station is a multi-purpose project for harnessing the Yangtze River, developed mainly for power generation, flood control, sediment flushing, and improving the navigation conditions in the reservoir and downstream.The Baihetan Hydropower Station is located on the downstream reaches of the Jinsha River.At the dam site, the most apparent stratigraphic lithology is basalt that belongs to the Emeishan (Emei Mt.) formation of the Permain System (P 2 β).In the construction of the plant's hydraulic structures, columnar jointed basalt is found to be widely distributed at the arch dam foundation, underground caverns, abutment slopes, and other water conservancy tunnels (Figure 11).Having obtained a 2-D Voronoi diagram with specified CD and CV, it can be extruded to a 3-D columnar shape in accordance with DD and DIP.Then for each extruding Voronoi cell and its length, a transverse joint can be generated by the determination of the intersection between the extruding Voronoi cell and the transverse plane.Sometimes, the transverse joint is non-penetrating, and the transverse joint can be selected with the specified TD and TP.In this way, columnar jointed rock mass with the site properties can be generated.

Engineering Geological Investigation
The Baihetan Hydropower Station is a multi-purpose project for harnessing the Yangtze River, developed mainly for power generation, flood control, sediment flushing, and improving the navigation conditions in the reservoir and downstream.The Baihetan Hydropower Station is located on the downstream reaches of the Jinsha River.At the dam site, the most apparent stratigraphic lithology is basalt that belongs to the Emeishan (Emei Mt.) formation of the Permain System (P2β).In the construction of the plant's hydraulic structures, columnar jointed basalt is found to be widely distributed at the arch dam foundation, underground caverns, abutment slopes, and other water conservancy tunnels (Figure 11).A typical figure of columnar jointed basalt is drawn by the Huadong Engineering Cooperation Limited (ECIDI), as shown in Figure 12 with related parameters listed in Table 1.The orientation of columnar is regular with Dip = 72 • and DD = 145 • , and the length of columnar is about 1.5 m.Based on the skeleton figure, a digital image processing technique is employed in the analysis.The area of each columnar is calculated and the coefficient of variation is obtained.Using a sample of 10 typical columnar basalt drawings, the statistical mean of coefficient of variation was about 44.18% and it implied that the discreteness of columnar basalt was quite large.Furthermore, the density of columnar cells could also be obtained from skeleton drawings, which gave a density of about 20 in a 1 m × 1 m square.A typical figure of columnar jointed basalt is drawn by the Huadong Engineering Cooperation Limited (ECIDI), as shown in Figure 12 with related parameters listed in Table 1.The orientation of columnar is regular with Dip = 72° and DD = 145°, and the length of columnar is about 1.5 m.Based on the skeleton figure, a digital image processing technique is employed in the analysis.The area of each columnar is calculated and the coefficient of variation is obtained.Using a sample of 10 typical columnar basalt drawings, the statistical mean of coefficient of variation was about 44.18% and it implied that the discreteness of columnar basalt was quite large.Furthermore, the density of columnar cells could also be obtained from skeleton drawings, which gave a density of about 20 in a 1 m × 1 m square.

Columnar Jointing Model Generation
Considering the scale effect of heterogeneous materials, a proper size has to be selected.For convenience, the size of the numerical model was taken as 1 m × 1 m.In order to ensure that the 3-D model could be generated, an initial area of 5 m × 5 m was chosen for the generation of a 2-D Voronoi diagram.From the density of site conditions, it was estimated that about 500 seeds were needed.A random Voronoi diagram was generated in Figure 13a with a CV of 56.45%.With a specified CV value of 44.18%, a Voronoi diagram could be generated (Figure 13b) by applying the algorithm proposed in Section 4. Then an extrusion at a certain direction with parameter DIP and DD was made to extend the rock from 2-D to 3-D (Figure 13c).Afterwards, horizontal hidden joints was implemented to cut the rock (Figure 13d).Using these operations, the coordinates information of all rock blocks could be calculated.By trimming the solid with a certain boundary, the block information of the columnar jointed rock is shown in Figure 13e.It was easier for the particle model in Figure 13f, where the model could be generated by importing the DFN (discrete fracture networks) in Figure 13d to a cubic particle model, which could be done easily in DEM packages like PFC (Particle Flow Code), YADE (Yet another dynamic engine), and LIGGGHTS (LAMMPS improved for general granular and granular heat transfer simulations).

Columnar Jointing Model Generation
Considering the scale effect of heterogeneous materials, a proper size has to be selected.For convenience, the size of the numerical model was taken as 1 m × 1 m.In order to ensure that the 3-D model could be generated, an initial area of 5 m × 5 m was chosen for the generation of a 2-D Voronoi diagram.From the density of site conditions, it was estimated that about 500 seeds were needed.A random Voronoi diagram was generated in Figure 13a with a CV of 56.45%.With a specified CV value of 44.18%, a Voronoi diagram could be generated (Figure 13b) by applying the algorithm proposed in Section 4. Then an extrusion at a certain direction with parameter DIP and DD was made to extend the rock from 2-D to 3-D (Figure 13c).Afterwards, horizontal hidden joints was implemented to cut the rock (Figure 13d).Using these operations, the coordinates information of all rock blocks could be calculated.By trimming the solid with a certain boundary, the block information of the columnar jointed rock is shown in Figure 13e.It was easier for the particle model in Figure 13f, where the model could be generated by importing the DFN (discrete fracture networks) in Figure 13d to a cubic particle model, which could be done easily in DEM packages like PFC (Particle Flow Code), YADE (Yet another dynamic engine), and LIGGGHTS (LAMMPS improved for general granular and granular heat transfer simulations).Columnar jointed rock mass with different CV is shown in Figure 14.It can be seen that the volume of columnar became more and more uniform with the decrease of coefficient of variation.Columnar jointed rock masses with different directions are shown in Figure 15.The result shows that the direction was also an important parameter affecting the morphology of columnar jointed rock masses.From all these cases, it can be concluded that this method is applicable and effective in the Columnar jointed rock mass with different CV is shown in Figure 14.It can be seen that the volume of columnar became more and more uniform with the decrease of coefficient of variation.Columnar jointed rock masses with different directions are shown in Figure 15.The result shows that the direction was also an important parameter affecting the morphology of columnar jointed rock masses.From all these cases, it can be concluded that this method is applicable and effective in the modeling of columnar jointed rock mass with complex structures, which is extremely helpful in the analysis of physical and mechanical properties of such materials.

Conclusions and Discussion
Aimed at the analysis of columnar jointed rock mass with complex structures, this paper proposed a novel method for the generation of numerical model.From the results obtained, the following conclusions can be drawn: (1) The coefficient of variation is an effective parameter for representing the deviation between the

Conclusions and Discussion
Aimed at the analysis of columnar jointed rock mass with complex structures, this paper proposed a novel method for the generation of numerical model.From the results obtained, the following conclusions can be drawn: (1) The coefficient of variation is an effective parameter for representing the deviation between the

Conclusions and Discussion
Aimed at the analysis of columnar jointed rock mass with complex structures, this paper proposed a novel method for the generation of numerical model.From the results obtained, the following conclusions can be drawn:

Symmetry 2018 ,
10, x FOR PEER REVIEW 3 of 15 squares are red.From the introduction of columnar jointed rock mass above, a diagrammatic drawing is proposed in Figure3.It is appropriate to model columnar joints using Voronoi diagrams.

Figure 3 .
Figure 3. Diagrammatic drawing of the joints in columnar basalt: (a) vertical columnar joint, and (b) horizontal transverse joint.
(b) All the open cells with a vertex outside the domain D are identified (shown in Figure 5c).(c) A set of new seeds symmetric to the seeds of open cells with respect to the domain boundary are created (Figure 5d).(d) New Voronoi diagram is generated with a classical tessellation (Figure 5e).(e) By removing the open cells as well as the related seeds, the final Voronoi diagram is shown in Figure 5f.Voronoi diagram, a Delaunay triangulation is first implemented.Then the perpendicular bisector is plotted to partition the domain into different Voronoi cells.With the development of computational geometry, Voronoi tessellation is included in many software and packages like MATLAB, Mathematica, and SciPy.However, there are two obvious shortcomings of the classical Voronoi diagram: (a) The Voronoi cell is not closed.For the Voronoi diagram in Figure 4, only one cell is closed and the other nine cells are open.This brings inconvenience for the analysis.(b) The shape of the Voronoi cell is random and it is very hard to generate a Voronoi diagram with a specified statistical distribution.
(b) All the open cells with a vertex outside the domain D are identified (shown in Figure 5c).(c) A set of new seeds symmetric to the seeds of open cells with respect to the domain boundary are created (Figure 5d).(d) New Voronoi diagram is generated with a classical tessellation (Figure 5e).

Figure 4 .
Figure 4.An illustration of Voronoi tessellation with 10 generators: (a) a set of generators, and (b) the resulting Voronoi diagram.

Figure 4 .
Figure 4.An illustration of Voronoi tessellation with 10 generators: (a) a set of generators, and (b) the resulting Voronoi diagram.
(b) All the open cells with a vertex outside the domain D are identified (shown in Figure 5c).(c) A set of new seeds symmetric to the seeds of open cells with respect to the domain boundary are created (Figure 5d).(d) New Voronoi diagram is generated with a classical tessellation (Figure 5e).

Figure 8 .
Figure 8. Illustration of the sampling method: (a) Voronoi tessellation, and (b) random sampling points.

Figure 8 .
Figure 8. Illustration of the sampling method: (a) Voronoi tessellation, and (b) random sampling points.

Figure 11 .
Figure 11.Site conditions of the Baihetan Hydropower Station: (a) construction site, and (b) typical columnar jointed basalt.

Figure 11 .
Figure 11.Site conditions of the Baihetan Hydropower Station: (a) construction site, and (b) typical columnar jointed basalt.

Figure 11 .
Figure 11.Site conditions of the Baihetan Hydropower Station: (a) construction site, and (b) typical columnar jointed basalt.

Table 1 .
Parameters of columnar jointed rock mass in dam foundation of Baihetan Hydropower Station.

Table 1 .
Parameters of columnar jointed rock mass in dam foundation of Baihetan Hydropower Station.