# Generating Irregular Models for 3D Spherical-Particle-Based Numerical Methods

Authors: Gang-Hai Huang; Yu-Yong Jiao; Xiu-Li Zhang; Amoussou-Coffi Adoko; Shu-Cai Li
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The realistic representation of an irregular geological body is essential to the construction of a particle simulation model. A three-dimensional (3D) sphere generator for an irregular model (SGIM), which is based on the platform of Microsoft Foundation Classes (MFC) in VC++, is developed to accurately simulate the inherent discontinuities in geological bodies. OpenGL is employed to visualize the modeling in the SGIM. Three key functions, namely, the basic-model-setup function, the excavating function, and the cutting function, are implemented. An open-pit slope is simulated using the proposed model. The results demonstrate that an extremely irregular 3D model of a geological body can be generated using the SGIM and that various types of discontinuities can be inserted to cut the model. The data structure of the model that is generated by the SGIM is versatile and can be easily modified to match various numerical calculation tools. This can be helpful in the application of particle simulation methods to large-scale geoengineering projects.


---

## Body

## 1. Introduction

Although particle-based methods such as the discrete element method (DEM) [1–4] and discontinuous deformation analysis (DDA) [5, 6] have been used to solve numerous types of engineering-oriented problems [7–13], the following fundamental simulation issues should be addressed to ensure that these methods solve static and quasistatic problems in a precise and efficient manner [14]. The first simulation issue involves the acquisition of the particle-scale mechanical properties of a particle from the measured macroscopic mechanical properties of rocks. The second simulation issue is that fictitious time instead of physical time is used in the particle simulation of a quasistatic problem. The third issue is that the conventional loading procedure used in the distinct element method is conceptually inaccurate from a force propagation point of view. After these fundamental problems were successfully addressed, a new type of particle-based method, known as the particle simulation method [14, 15], was proposed and implemented to solve several types of geological and engineering problems associated with large-scale static, quasistatic, and dynamic systems [16–23].The first step in applying particle simulation methods to solve any practical engineering problem involves the generation of a particle model of discrete objects, which are packed in a form that realistically represents the physical model of the problem. Extensive research has been conducted on this topic [24–28], in which the shapes of the particles include disc (under 2D conditions), sphere, ellipsoid, and tetrahedron shapes. The use of simplified geometric entities to model discrete objects has been demonstrated to provide acceptable approximations of complex physical phenomena [25]. Due to the simple contact detection of spheres, it has been extensively used as a basic element to compute different types of discrete models. The packing of spherical particles in a certain domain has become the focus of research in recent years.The majority of studies on sphere generation algorithms focus on ways to efficiently generate spheres and pack them in a container [24, 25, 29–32]. These algorithms are usually applied in the DEM or DDA. The algorithm described in [31] generates small balls that are enlarged until the domain is completely filled. Under the action of gravity, the balls fall, roll, and come to rest in a stable and dense state. This sphere-generation algorithm is used extensively in the popular particle flow code (PFC) [33]. A shortcoming of ordinary algorithms is that they are time consuming. To overcome this disadvantage, an advancing front face (AFF) algorithm for particle packing was also proposed [32]. The application of a spherical particle model to simulate a geological body is challenging because geological bodies frequently contain extremely irregular surfaces. As the main component of a geological body, a rock mass is not a homogeneous material due to numerous discontinuities (joints and faults). The failure mode of a geological body highly depends on the distribution of discontinuities; the influence of discontinuities requires consideration. Research on this aspect has been conducted by Scholtès and Donzé [9] and John et al. [12]. A 3D numerical model that is based on the discrete element method, in which preexisting discontinuities are initially inserted into the particle model, was used to model the progressive failure of a jointed rock mass. However, because these models are extremely simple, they inadequately represent a geological body. Although several similar studies have attempted to consider discontinuities, they do not consider the complexity of a geological body, which is a basic requirement to successfully obtain realistic particle simulation models for geoengineering projects.Few studies focus on the implementation of an irregular body into spherical particle modeling for complex engineering, especially large-scale geological body models with extremely irregular surfaces. The majority of existing methods are predominantly applied to the secondary development of certain existing programs [7] or commercial software [9, 12], for which development of a complex 3D model to solve numerous modeling problems in real engineering projects is challenging.To overcome these difficulties, the platform of Microsoft Foundation Classes (MFC) in VC++ is used to develop a program, which is referred to as a 3D sphere generator for irregular models (SGIM) in this paper. Step-by-step visualization modeling operations can be implemented using this program. Different types of discontinuities can be inserted to cut the model to obtain a more realistic rock mass, which is significant to the application of particle simulation methods in numerical simulation of large-scale geoengineering problems.


## 2. Program Structure

The structure of the SGIM is illustrated in Figure1. The model parameters are inputted from a keyboard and mouse and subsequently translated into the basic model, the excavating faces and the discontinuities by the basic-model-setup function, the excavating-faces-generating function, and the discontinuities-generating function, respectively. The basic model is excavated by excavating faces via the excavating function and is cut by discontinuities via the cutting function. Thus, the basic-model-setup, excavating, and cutting functions are key elements in the program. A model with spherical particles, which can perform numerical calculations, is obtained. By utilizing the OpenGL-drawing function, various components of the model are displayed on a computer screen. By employing the model-showing parameter-operation function, the model can be moved, rotated, and displayed in different scales and in animation mode. All or part of the components can be hidden or displayed translucently.Figure 1
Program structure of the SGIM.


## 3. Algorithms of the Three Key Functions

### 3.1. Basic-Model-Setup Function

The first function, which is referred to as the basic-model-setup function, facilitates the development of the basic model. It comprises a cuboid with edges parallel to the coordinate axes. In the current version of the SGIM, the equal spheres packing method is applied. All balls in the model contain the same radii and are packed by row, column, and layer. Note that selection of the particle size is a fundamental problem in modeling with spherical particles. An arrangement of spheres with equivalent sizes may result in a lattice effect on the mechanical behavior of the particle assembly. For this reason, the random packing method is usually recommended. However, implementation of this method is time consuming for a large-scale geological body. In the proposed model, a small particle size was adopted, which can reduce the lattice effect. The optimal size was based on a trial-and-error method. The number of rows, columns, and layers is obtained as(1)I=z2-z12r,J=y2-y12r,K=x2-x12r,(I,J,andKroundingtointegers),
where I, J, and K are the number of rows, columns, and layers, respectively; x2-x1, y2-y1, and z2-z1 comprise the intervals that define the model in the Cartesian coordinates x, y, and z, respectively; and r is the radius of the sphere.The center coordinates of a sphere are defined as(2)x=x1+r(2k-1),y=y1+r(2j-1),z=z1+r(2i-1),
where x, y, and z are the three components of the center coordinates of the sphere, and i, j, and k are the layer, the row and the column, respectively, in which the sphere is located (i=1,2,…,I; j=1,2,…,J; and k=1,2,…,K).


### 3.2. Excavating Function

The second function is used to excavate the basic model. Excavating faces are composed of a series of interconnected triangular facets. These facets are obtained by inputting the coordinates of the points located along the edges of the excavating faces. The principle of an excavating function is to delete the spheres of the basic model if their centers are located in excavating spaces.As shown in Figure2, three noncollinear points L0, M0, and N0 compose a triangular facet (denoted by ΔL0M0N0). ΔL0M0N0 denotes the intersection of two spaces: space I and space II. Space I (blue color) is defined as the region located inside the regular infinite-length prism in which the distance to ΔL0M0N0 is smaller than r. Space II (pink color) is defined as the region situated inside the vertical prism and above ΔL0M0N0. The excavating space associated with ΔL0M0N0 (denoted by Ω) is defined as Ω=I∪II. The condition that enables the center of one sphere to be located in excavating space Ω is calculated by the following seven steps.Figure 2
Excavating space of a triangular facet.Let the coordinates of the three points beL0(x01,y01,z01), M0(x02,y02,z02), and N0(x03,y03,z03) and the center coordinates of the sphere, O(x0,y0,z0).(1) The equation of the planeL0M0N0 can be written as
(3)A0x+B0y+C0z+D0=0,
where
(4)A0=|y01z011y02z021y03z031|,B0=-|x01z011x02z021x03z031|,C0=|x01y011x02y021x03y031|,D0=-|x01y01z01x02y02z02x03y03z03|,
and n→1={A0,B0,C0} is the normal vector of the plane.(2) The coordinates of the three pointsL1(x11,y11,z11), M1(x12,y12,z12), and N1(x13,y13,z13), which determine space I (see Figure 2), are obtained by
(5)x11=x01+A0,y11=y01+B0,z11=z01+C0,x12=x02+A0,y12=y02+B0,z12=z02+C0,x13=x03+A0,y13=y03+B0,z13=z03+C0.(3) The equations of the planes that correspond to the three sides of the prismL0M0N0-L1M1N1 can be obtained using the method described in step (1) as follows:
(6)planeL0L1M0:A1x+B1y+C1z+D1=0,planeM0M1N0:A2x+B2y+C2z+D2=0,planeL0L1N0:A3x+B3y+C3z+D3=0.
(4)  O(x0,y0,z0) is located in space I if
(7)|A1x0+B1y0+C1z0+D1|≤|D1-D10|,|A1x0+B1y0+C1z0+D10|≤|D1-D10|,|A2x0+B2y0+C2z0+D2|≤|D2-D20|,|A2x0+B2y0+C2z0+D20|≤|D2-D20|,|A3x0+B3y0+C3z0+D3|≤|D3-D30|,|A3x0+B3y0+C3z0+D30|≤|D3-D30|,|A0x0+B0y0+C0z0+D0|A02+B02+C02<r,
where
(8)D10=-A1x03-B1y03-C1z03,D20=-A2x01-B2y01-C2z01,D30=-A3x02-B3y02-C3z02.(5)  The coordinates of the three pointsL2(x21,y21,z21), M2(x22,y22,z22), and N2(x23,y23,z23), which determine space II (see Figure 2), are obtained by
(9)x21=x01,y21=y01,z21=z01+1,x22=x02,y22=y02,z22=z02+1,x23=x03,y23=y03,z23=z03+1.(6) The equations for the planes that correspond to the three sides of the prismL0M0N0-L2M2N2 can be written as follows:
(10)planeL0L2M0:A4x+B4y+C4z+D4=0,planeM0M2N0:A5x+B5y+C5z+D5=0,planeL0L2N0:A6x+B6y+C6z+D6=0.(7)O(x0,y0,z0) is located in space II if
(11)|A4x0+B4y0+C4z0+D4|≤|D4-D40|,|A4x0+B4y0+C4z0+D40|≤|D4-D40|,|A5x0+B5y0+C5z0+D5|≤|D5-D50|,|A5x0+B5y0+C5z0+D50|≤|D5-D50|,|A6x0+B6y0+C6z0+D6|≤|D6-D60|,|A6x0+B6y0+C6z0+D60|≤|D6-D60|,z0≥-A0x0+B0y0+D0C0,
where
(12)D40=-A4x03-B4y03-C4z03,D50=-A5x01-B5y01-C5z01,D60=-A6x02-B6y02-C6z02.Note that, whenC0=0, ΔL0M0N0 is a vertical plane, and space II becomes a face. Thus, C0≠0 is a prerequisite for the calculation of space II.The excavating space of one triangular facet can be calculated with these seven steps. The total excavating spaces are the union of the excavating spaces of all triangular facets.


### 3.3. Cutting Function

The cutting function, which is employed to cut the basic model with discontinuities, is based on certain assumptions. Several studies indicate that nearly equivalent lengths of discontinuities exist in the strike direction and the dip direction. As a result, the hypothesis that the discontinuity in 3D space can be described as a thin disc is often adopted [34, 35]. In this study, this same hypothesis is employed, and the discontinuities are divided into two types: the first type is zero-thickness discontinuity, and the second type is non-zero-thickness discontinuity. The objective is to obtain a realistic rock mass by cutting intact rock mass with discontinuities. The spheres located in the domain in which discontinuities cut through the rock mass are marked in the SGIM with specific colors, such as red, blue, and green. For a discontinuity with thickness, the spheres are shown in red if the distances from their centers to the discontinuity are smaller than half of the thickness, whereas the spheres adhering to the two sides of the discontinuity are shown in green and blue for zero-thickness discontinuities. When the spherical particle model is used for numerical calculation in the DEM or DDA, the bond strength between two spheres is dependent on the colors of the spheres. No bond forms if the spheres are green blue or if at least one of the spheres is red. The influence of discontinuities on the model is reflected in the mechanical calculations.Figure3 provides a spatial representation of a discontinuity using a disc plane. As seen in Figure 3, discontinuities with radius R, thickness d, coordinates of center M(x1,y1,z1), and normal vector S→={A,B,C} exist. The condition that one sphere is marked with red, green, or blue is calculated by the following steps, where the center of the sphere is assumed O(x0,y0,z0).(1)
The center of the sphere must be located inside the infinite-length cylinder with one of its normal sections as the discontinuity. This is guaranteed by (13) as follows:
(13)|MO→×S→||S→|≤R.(2)
Whend≠0, the sphere is marked with red if
(14)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤d2.(3)
Whend=0, the sphere is marked with green if
(15)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤2r,A(x1-x0)+B(y1-y0)+C(z1-z0)≥0.(4)
Whend=0, the sphere is marked with blue if
(16)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤2r,A(x1-x0)+B(y1-y0)+C(z1-z0)<0.Figure 3
A disc plane that represents discontinuity in 3D space.


## 3.1. Basic-Model-Setup Function

The first function, which is referred to as the basic-model-setup function, facilitates the development of the basic model. It comprises a cuboid with edges parallel to the coordinate axes. In the current version of the SGIM, the equal spheres packing method is applied. All balls in the model contain the same radii and are packed by row, column, and layer. Note that selection of the particle size is a fundamental problem in modeling with spherical particles. An arrangement of spheres with equivalent sizes may result in a lattice effect on the mechanical behavior of the particle assembly. For this reason, the random packing method is usually recommended. However, implementation of this method is time consuming for a large-scale geological body. In the proposed model, a small particle size was adopted, which can reduce the lattice effect. The optimal size was based on a trial-and-error method. The number of rows, columns, and layers is obtained as(1)I=z2-z12r,J=y2-y12r,K=x2-x12r,(I,J,andKroundingtointegers),
where I, J, and K are the number of rows, columns, and layers, respectively; x2-x1, y2-y1, and z2-z1 comprise the intervals that define the model in the Cartesian coordinates x, y, and z, respectively; and r is the radius of the sphere.The center coordinates of a sphere are defined as(2)x=x1+r(2k-1),y=y1+r(2j-1),z=z1+r(2i-1),
where x, y, and z are the three components of the center coordinates of the sphere, and i, j, and k are the layer, the row and the column, respectively, in which the sphere is located (i=1,2,…,I; j=1,2,…,J; and k=1,2,…,K).


## 3.2. Excavating Function

The second function is used to excavate the basic model. Excavating faces are composed of a series of interconnected triangular facets. These facets are obtained by inputting the coordinates of the points located along the edges of the excavating faces. The principle of an excavating function is to delete the spheres of the basic model if their centers are located in excavating spaces.As shown in Figure2, three noncollinear points L0, M0, and N0 compose a triangular facet (denoted by ΔL0M0N0). ΔL0M0N0 denotes the intersection of two spaces: space I and space II. Space I (blue color) is defined as the region located inside the regular infinite-length prism in which the distance to ΔL0M0N0 is smaller than r. Space II (pink color) is defined as the region situated inside the vertical prism and above ΔL0M0N0. The excavating space associated with ΔL0M0N0 (denoted by Ω) is defined as Ω=I∪II. The condition that enables the center of one sphere to be located in excavating space Ω is calculated by the following seven steps.Figure 2
Excavating space of a triangular facet.Let the coordinates of the three points beL0(x01,y01,z01), M0(x02,y02,z02), and N0(x03,y03,z03) and the center coordinates of the sphere, O(x0,y0,z0).(1) The equation of the planeL0M0N0 can be written as
(3)A0x+B0y+C0z+D0=0,
where
(4)A0=|y01z011y02z021y03z031|,B0=-|x01z011x02z021x03z031|,C0=|x01y011x02y021x03y031|,D0=-|x01y01z01x02y02z02x03y03z03|,
and n→1={A0,B0,C0} is the normal vector of the plane.(2) The coordinates of the three pointsL1(x11,y11,z11), M1(x12,y12,z12), and N1(x13,y13,z13), which determine space I (see Figure 2), are obtained by
(5)x11=x01+A0,y11=y01+B0,z11=z01+C0,x12=x02+A0,y12=y02+B0,z12=z02+C0,x13=x03+A0,y13=y03+B0,z13=z03+C0.(3) The equations of the planes that correspond to the three sides of the prismL0M0N0-L1M1N1 can be obtained using the method described in step (1) as follows:
(6)planeL0L1M0:A1x+B1y+C1z+D1=0,planeM0M1N0:A2x+B2y+C2z+D2=0,planeL0L1N0:A3x+B3y+C3z+D3=0.
(4)  O(x0,y0,z0) is located in space I if
(7)|A1x0+B1y0+C1z0+D1|≤|D1-D10|,|A1x0+B1y0+C1z0+D10|≤|D1-D10|,|A2x0+B2y0+C2z0+D2|≤|D2-D20|,|A2x0+B2y0+C2z0+D20|≤|D2-D20|,|A3x0+B3y0+C3z0+D3|≤|D3-D30|,|A3x0+B3y0+C3z0+D30|≤|D3-D30|,|A0x0+B0y0+C0z0+D0|A02+B02+C02<r,
where
(8)D10=-A1x03-B1y03-C1z03,D20=-A2x01-B2y01-C2z01,D30=-A3x02-B3y02-C3z02.(5)  The coordinates of the three pointsL2(x21,y21,z21), M2(x22,y22,z22), and N2(x23,y23,z23), which determine space II (see Figure 2), are obtained by
(9)x21=x01,y21=y01,z21=z01+1,x22=x02,y22=y02,z22=z02+1,x23=x03,y23=y03,z23=z03+1.(6) The equations for the planes that correspond to the three sides of the prismL0M0N0-L2M2N2 can be written as follows:
(10)planeL0L2M0:A4x+B4y+C4z+D4=0,planeM0M2N0:A5x+B5y+C5z+D5=0,planeL0L2N0:A6x+B6y+C6z+D6=0.(7)O(x0,y0,z0) is located in space II if
(11)|A4x0+B4y0+C4z0+D4|≤|D4-D40|,|A4x0+B4y0+C4z0+D40|≤|D4-D40|,|A5x0+B5y0+C5z0+D5|≤|D5-D50|,|A5x0+B5y0+C5z0+D50|≤|D5-D50|,|A6x0+B6y0+C6z0+D6|≤|D6-D60|,|A6x0+B6y0+C6z0+D60|≤|D6-D60|,z0≥-A0x0+B0y0+D0C0,
where
(12)D40=-A4x03-B4y03-C4z03,D50=-A5x01-B5y01-C5z01,D60=-A6x02-B6y02-C6z02.Note that, whenC0=0, ΔL0M0N0 is a vertical plane, and space II becomes a face. Thus, C0≠0 is a prerequisite for the calculation of space II.The excavating space of one triangular facet can be calculated with these seven steps. The total excavating spaces are the union of the excavating spaces of all triangular facets.


## 3.3. Cutting Function

The cutting function, which is employed to cut the basic model with discontinuities, is based on certain assumptions. Several studies indicate that nearly equivalent lengths of discontinuities exist in the strike direction and the dip direction. As a result, the hypothesis that the discontinuity in 3D space can be described as a thin disc is often adopted [34, 35]. In this study, this same hypothesis is employed, and the discontinuities are divided into two types: the first type is zero-thickness discontinuity, and the second type is non-zero-thickness discontinuity. The objective is to obtain a realistic rock mass by cutting intact rock mass with discontinuities. The spheres located in the domain in which discontinuities cut through the rock mass are marked in the SGIM with specific colors, such as red, blue, and green. For a discontinuity with thickness, the spheres are shown in red if the distances from their centers to the discontinuity are smaller than half of the thickness, whereas the spheres adhering to the two sides of the discontinuity are shown in green and blue for zero-thickness discontinuities. When the spherical particle model is used for numerical calculation in the DEM or DDA, the bond strength between two spheres is dependent on the colors of the spheres. No bond forms if the spheres are green blue or if at least one of the spheres is red. The influence of discontinuities on the model is reflected in the mechanical calculations.Figure3 provides a spatial representation of a discontinuity using a disc plane. As seen in Figure 3, discontinuities with radius R, thickness d, coordinates of center M(x1,y1,z1), and normal vector S→={A,B,C} exist. The condition that one sphere is marked with red, green, or blue is calculated by the following steps, where the center of the sphere is assumed O(x0,y0,z0).(1)
The center of the sphere must be located inside the infinite-length cylinder with one of its normal sections as the discontinuity. This is guaranteed by (13) as follows:
(13)|MO→×S→||S→|≤R.(2)
Whend≠0, the sphere is marked with red if
(14)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤d2.(3)
Whend=0, the sphere is marked with green if
(15)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤2r,A(x1-x0)+B(y1-y0)+C(z1-z0)≥0.(4)
Whend=0, the sphere is marked with blue if
(16)|A(x1-x0)+B(y1-y0)+C(z1-z0)|A2+B2+C2≤2r,A(x1-x0)+B(y1-y0)+C(z1-z0)<0.Figure 3
A disc plane that represents discontinuity in 3D space.


## 4. Modeling Example

The simulated geological body and its discontinuities are illustrated in this example. This model tries to visualize the body in three dimensions using the proposed algorithm. The spherical particle calculation model of an open-pit slope in the Zijinshan gold-copper mine, which is the largest gold mine of the Zijin Mining Group in China, was constructed. The assumption that the geological body contains two sets of joints and a fault was applied. The dip and strike directions of these discontinuities are arbitrary in this study.


### 4.1. Setup of the Basic Model

Figure4 illustrates a basic model with a sphere of radius r=0.5 that is located inside the three coordinate intervals x=-150 to 150, y=-125 to 125, and z=0 to 140.Figure 4
The basic model.


### 4.2. Excavating the Basic Model

The basic model, which is excavated by adding excavating faces, is shown in Figure5. An irregular model of a geological body with man-made open work in a hills landform is obtained. Figure 5(a) includes the site picture, and Figure 5(b) includes the spherical particle model. The model is similar to the site picture.The model after excavating.
(a)
The site picture
(b)
The spherical particle calculation model


### 4.3. Insertion of Discontinuities

Two groups of discontinuities with zero thickness and one group of discontinuity with a thickness of 4.0 are inserted into the model, as shown in Figure6. The reason for using zero-thickness groups is to simulate two perpendicular sets of rock mass joints, whereas the reason for considering the group with thickness is to model the fault. The colored spheres that are produced by insertion of the discontinuities are presented in Figure 7. In this figure, the transversal discontinuity (fault) is represented by red spheres, and the joints are represented by green and blue spheres.Figure 6
Model after insertion of discontinuities.Figure 7
Resulting spheres with different colors that adhere to discontinuities.


### 4.4. The Final Spherical Particle Model

The final spherical particle model is shown in Figure8. This model represents a realistic geological body with different types of discontinuities, which can be used for numerical calculations in the DEM or DDA. The discontinuities with thickness are suitable for modeling large-scale faults in a geological body, whereas zero-thickness discontinuities can be used to model random joints in a jointed rock mass. In addition, for the purpose of effectively and efficiently simulating the infinite extension of a geological system, dynamic and transient infinite elements [36] can be added to the outer boundaries of this model in the future research.Figure 8
Spherical particle model for numerical calculation.


## 4.1. Setup of the Basic Model

Figure4 illustrates a basic model with a sphere of radius r=0.5 that is located inside the three coordinate intervals x=-150 to 150, y=-125 to 125, and z=0 to 140.Figure 4
The basic model.


## 4.2. Excavating the Basic Model

The basic model, which is excavated by adding excavating faces, is shown in Figure5. An irregular model of a geological body with man-made open work in a hills landform is obtained. Figure 5(a) includes the site picture, and Figure 5(b) includes the spherical particle model. The model is similar to the site picture.The model after excavating.
(a)
The site picture
(b)
The spherical particle calculation model


## 4.3. Insertion of Discontinuities

Two groups of discontinuities with zero thickness and one group of discontinuity with a thickness of 4.0 are inserted into the model, as shown in Figure6. The reason for using zero-thickness groups is to simulate two perpendicular sets of rock mass joints, whereas the reason for considering the group with thickness is to model the fault. The colored spheres that are produced by insertion of the discontinuities are presented in Figure 7. In this figure, the transversal discontinuity (fault) is represented by red spheres, and the joints are represented by green and blue spheres.Figure 6
Model after insertion of discontinuities.Figure 7
Resulting spheres with different colors that adhere to discontinuities.


## 4.4. The Final Spherical Particle Model

The final spherical particle model is shown in Figure8. This model represents a realistic geological body with different types of discontinuities, which can be used for numerical calculations in the DEM or DDA. The discontinuities with thickness are suitable for modeling large-scale faults in a geological body, whereas zero-thickness discontinuities can be used to model random joints in a jointed rock mass. In addition, for the purpose of effectively and efficiently simulating the infinite extension of a geological system, dynamic and transient infinite elements [36] can be added to the outer boundaries of this model in the future research.Figure 8
Spherical particle model for numerical calculation.


## 5. Conclusions and Discussion

The SGIM exhibits the following advantages: the SGIM does not rely on existing software, and the data structure of the model generated by the SGIM can be easily modified to match numerical simulation tools, such as PFC3D. By utilizing the SGIM, spheres can be generated for extremely irregular 3D models to reasonably consider discontinuities in complex engineering problems.Because the focus of this study is to develop a sphere generator for the 3D model for irregular geological bodies, some aspects are not considered in the present version of the SGIM and should be addressed in future studies. First, random packing of equal and unequal spheres should be adopted. In certain numerical simulation tools, such as PFC3D, ideal numerical simulation results can be obtained when the ratio of particle diameters is located within a certain range. Second, an algorithm to calculate the inner excavating space is required. Numerous excavating spaces frequently exist inside a geological body in practical engineering problems. However, the excavating function applied in this study is only suitable for excavation on the upper surface of the model. Third, the joint in a geological body is usually described as a thin disc. However, certain large-scale faults with complex shapes may exist in a geological body. A curved surface may be more suitable for this situation. Because some distribution laws of joints exist for a jointed rock mass, individually inputting all discontinuities is time consuming and impractical. Therefore, the incorporation of alternate ways to describe discontinuities, such as the Monte Carlo method, is required to simulate the network of discontinuities.



---

*Source: 302680-2013-09-03.xml*