Digital image-based numerical modeling method for prediction of inhomogeneous rock failure

ARTICLE IN PRESS

International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

Digital image-based numerical modeling method for prediction of inhomogeneous rock failure S. Chen, Z.Q. Yue*, L.G. Tham Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, People’s Republic of China Accepted 8 March 2004

Abstract This paper presents a two-dimensional digital image-based numerical modeling method for prediction of inhomogeneous rock failure behavior under loadings. Actual inhomogeneities of granitic rocks are extracted from color images of the granite crosssections. They are represented as the internal spatial distribution of three main granite minerals (quartz, feldspar and biotite). The actual mineral spatial distribution on granite cross-section is then incorporated into conventional numerical software packages to examine the rock failures under loading. Some digital image processing algorithms are presented to isolate and identify the main internal minerals and their distribution from color digital images. A simple method is proposed to transform the actual image data into vector data for generation of finite meshes or grids. The vector data are used directly as uniform square element meshes or grids that can be inputted into the existing software packages. The finite difference software package FLAC is used as an example for the present investigation. The conventional Mohr–Coulomb and tensile stress failure criteria are used to examine the failure behavior of a circular granite cross-section under the conventional Brazilian indirect tensile test loading conditions. The numerical results indicate that the vertical tensile crack initiates in a biotite located near the geometrical center of the granite cross-section and the actual spatial distribution of the three minerals plays an important role in modifying the propagation pattern of the tensile crack from its theoretical position at the central vertical diameter of a homogeneous circular cross-section. The numerically predicted failure load and tensile strength values for inhomogeneous granite are much lower than the expected values. r 2004 Elsevier Ltd. All rights reserved. Keywords: Rock mechanics; Granite; Minerals; Inhomogeneity; Digital image processing; Mesh generation; Finite difference method; Failure prediction

1. Introduction Rocks are inhomogeneous natural materials and widely encountered in geotechnical and mining engineering projects. Mechanical response of rock materials and particularly their failure predictions under external loading are of the most importance in the design and construction of various problems encountered in geotechnical and mining engineering. During last few decades, numerical tools have been proved to be an effective approach to predict the mechanical responses or failure process of rock materials under loading by taking into account their internal inhomogeneities. In many numerical models, the inhomogeneities and *Corresponding author. Tel.: +852-2859-1967; fax: +852-25595337. E-mail address: [email protected] (Z.Q. Yue). 1365-1609/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2004.03.002

microstructures of rock materials have been simulated via statistical tools or random distribution packages [1– 3]. These artificial internal inhomogeneities or microstructures are prevailing in rock mechanics because there had been lack of tools or methods to quantitatively establish the actual rock inhomogeneities. In recent years, digital image processing (DIP), i.e., manipulation of images by computers, has exhibited extensive applications in many disciplines of sciences and engineering. Some recent studies have shown that digitalized image data such as characteristics of cracks on a rock surface and distribution of different minerals in intact rocks can be used to establish the internal structures and inhomogeneities of natural rock materials. In this paper, we attempt to use digital image techniques to acquire the internal structures of different mineral distribution from color rock image and incorporate them into a commercial numerical software

ARTICLE IN PRESS 940

S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

package for mechanical analyses and failure prediction. Numerical results have shown that the proposed digital image-based numerical modeling method can take actual rock mineral distribution into account in the numerical prediction of rock failure process under external loading.

2. Background Rock, as a natural material on the Earth, is an inhomogeneous material in nature because it consists of many individual materials and components including minerals, grains, cement materials, voids and cracks. These individual materials and components usually have different mechanical and physical proprieties and different responses under internal and external loadings. It is well known that constitutive equations governing linear or non-linear relations between stresses and strains at a point are dependent on the properties of materials at that point. Therefore, actual spatial distributions of these individual components in a rock solid can play a dominant role in determining stress and strain distributions in the rock due to loading and further control the failure process that may be induced in the rock. Such basic knowledge and principle have been well recognized and accepted in rock mechanics communities around the world. During the last few decades, many laboratory tests, analytical studies and numerical modeling have been carried out to examine the engineering properties and mechanic responses of rock materials. An examination of relevant literature has shown that a majority of the investigations have been concentrated on macroscopic behavior of rocks. These studies established theoretical and numerical models of rock mechanics on the basis of general principles of continuum mechanics by considering rock materials as homogeneous or piecewise homogeneous materials and using the classical phenomenological approach. In these studies, details of the inhomogeneity and microstructure of rock solids are considered only on a phenomenological basis. In recent years, a number of researchers [1–9] have attempted to carry out numerical analysis of rock responses by taking into account rock material inhomogeneities. In these studies, rock inhomogeneities were incorporated into numerical simulations for examining mechanical behavior, crack initiation, crack propagation and coalescence of rocks under loadings. Some special models such as lattice model, particle model and rock failure process analysis (RFPA) model have been developed. Microstructures of rocks in these models are all generated by using statistical tools. These statistical tools can simulate numerical material inhomogeneity models that are statistically equivalent to those of actual rock materials with known statistical parameters. As

stated by Blair and Cook [3], these statistical methods have emerged as a promising technique for the analysis of heterogeneous rock materials. The so-called lattice model, which was used by Schlangen and Van Mier [10] for concrete studies, was adopted in rock mechanics to simulate the rock inhomogeneity. In this model, lattice consisting of regular or random triangle elements is generated in the rock material domain. Then particles of lumped masses are overlaid on this lattice net and were connected by massless springs along the edges of the triangle elements. The inhomogeneity of rock material can be simulated by assigning random strengths and sizes to the individual members of a lattice. A number of researchers such as Song and Kim [1], Li et al. [11] and Napier and Dede [12] used the lattice models to generate the rock inhomogeneity and examined the fracturing process of intact rock materials. In 1979, Cundall and Strack [13] proposed a classical discrete element method based on particle models for simulation of granular and particulate systems. The granular materials are represented by random package of two-dimensional or three-dimensional arbitrary size spherical hard particles. These particles bond together and move independently of one another and interact only at contacts. The particle model has been proved to be an ideal numerical tool to analyze rock materials [13–15]. Further applications of this model in rock mechanics can be found in a recent conference proceedings edited by Konietzky in 2002 [16]. Recently, Tang and his co-researchers [2,5] developed a finite element program RFPA technique to model the observed evolution of damage and induced seismic activities. Inhomogeneity in this model is introduced by defining an element stiffness/strength distribution via the Weibull statistical function. The Weibull statistical function used one homogeneity index m to describe the degree of inhomogeneity in rocks. A larger value of m implies a more heterogeneous material and vice versa. By adopting the same distribution function, Fang and Harrison [8,9] developed a new technique to simulate the brittle fracture in rocks under compressive loading using a commercial finite difference software— FLAC [17]. The above literature review may have indicated that most current studies on inhomogeneous rock have to establish rock microstructures as the basis of numerical analyses. These microstructures are artificial microstructures because they are generated from statistical tools or random distribution packages numerically. It is further noted that the above statistics-based mechanical analyses have revealed that it is significant to take into account the rock inhomogeneity and microstructure in the examination and modeling of the rock behavior and responses under loading. It seems, however, that there

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

are few publications available in the relevant literature to compare the extent and accuracy of the simulated inhomogeneities and microstructures with actual rock material inhomogeneities. It is believed that the statistical results cannot exactly represent the actual rock inhomogeneities and microstructures. Statistics-based modeling methods have been prevailing in rock mechanics mainly because it was hard to look inside the rock materials and obtain accurate measurements of actual geometrical distributions of minerals and other components in rock solids. In order to have a better representation of the geomaterial inhomogeneities and microstructures, Sepher et al. [18] manually measured the microstructure of an asphalt concrete and then carried out a finite element analysis of an asphalt pavement under traffic loading. Such manual measurements of geomaterial microstructures have limitations since they are time consuming and lack of control in measurement accuracy. There is a need that we should develop new techniques to incorporate actual rock inhomogeneities and microstructures into numerical analyses. For many years in science history, images have been used as a tool to examine material internal structures and flaws because they can truly reflect the geometrical distribution of different materials through different degrees of gray levels or colors. In recent years, numerical techniques (hardware and software) have been developed to transfer image into digital forms. As a result, DIP has become a powerful new technology that can provide a new means for accurate and quantitative measurements of inhomogeneities of rock and other natural made materials. Generally speaking, DIP is the term applied to convert an image into a digital form, and then apply mathematical algorithms to extract significant information from the image. With the rapid development of computer hardware and digital image theories, DIP has been used widely for the morphological features in many fields of sciences and engineering including medical science, aviation, biology and civil engineering [19–31]. In particular, Yue and his co-researchers [24,25] carried out quantitative investigations of the orientations, distributions and shapes of aggregates in asphalt concrete using an innovative DIP procedure. Kwan and his co-researchers [26,27] examined the particle shape characteristics of coarse aggregates and calculated their flakiness and elongations using DIP-based methods. The DIP also has been successfully applied in rock mechanics. For instance, Reid and Harrison [28] used a series of DIP techniques to detect discontinuity geometry from the rock mass surface. Hadjigeorgiou et al. [23] developed a digital face mapping system to acquire the discontinuity network data from rock exposure images. Li [6] and Tham et al. [7] used a double-replica technique to analyze the image of Hong

941

Kong granite and revealed the failure process and cracking behavior of granite samples under uniaxial compression. The present literature review indicates that most applications of DIP in rock mechanics have been mainly concentrated on quantitative analysis of material inhomogeneities and observations of geometrical features on rock surfaces. Few publications are available in incorporating DIP into mechanical modeling. Digital image-based numerical modeling can be considered as a new approach for geometry construction of material internal structures and inhomogeneities. Recently, Yue and his co-researchers [31–33] have launched a research project to develop techniques to examine mechanical responses of rock and other geomaterials by taking into account their actual internal structures and inhomogeneities. Yue et al. [31] developed several DIP techniques, region-growing method and edge detection method to establish the actual microstructures of geomaterials. The actual microstructures are then converted into vector polygon data via geometry transformation algorithms. An automatic mesh-generation technique is then employed to generate finite element meshes. The meshes are incorporated into a linear finite element software for the examination of mechanical responses of asphalt concrete samples in Brazilian tensile tests. Yue et al. [32,33] further applied this method to analyze effect of minerals on the linear mechanical response of granites under loading and to examine effect of inhomogeneous soils on seepage flows, respectively. The main objective of this paper is to present a digital image-based numerical modeling method for the prediction of rock material failures taking into account actual mineral distribution in the rocks. Specifically, DIP techniques are developed for quantifying actual distribution of main minerals in color rock cross-section images. The actual distribution results are used for automatical mesh generation and can be incorporated into conventional numerical software packages for various mechanical analyses and failure predictions. The presentation of this paper can be outlined as follows: (a) New DIP techniques for color images are introduced for acquiring the actual inhomogeneities and microstructures from color image of rock material cross-sections. The main minerals such as feldspar, biotite and quartz and their actual geometrical distributions in granites are identified specifically. (b) Digital image data that contain the rock internal inhomogeneity information are then transformed into a new type of vector data. It is found that the data structures of digital images are exactly the same as those of conventional square element

ARTICLE IN PRESS 942

S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

meshes. Therefore, the vector data can be used automatically as square element meshes or grids in conventional numerical models. (c) The square element meshes are reformatted to incorporate them into existing software packages for mechanical analysis and failure prediction. (d) Finally, the commercial numerical software package FLAC is used to show the applicability of the proposed approach involving finite difference method (FDM) for failure prediction. The conventional Brazilian tensile testing is used to show the prediction of tensile crack initiation and propagation in a circular granite cross-section on which the actual mineral distribution is taken into account in the numerical modeling.

3. Digital image measurement method In the ensuing, fundamental aspects of digital images are introduced first. Then algorithms to extract spatial distribution data from color images are presented. New DIP techniques based on the HSI color space are further developed to measure and identify the multi-individual minerals in granitic rocks. 3.1. Digital image and image data Field cores or laboratory prepared rock samples can be cut into multiple vertical or horizontal plane crosssections with a circular masonry saw. A ruler is placed beneath the section for scaling and calibration. The fresh cross-sections are then photographed with either a digital camera or a conventional camera or a digital scanner. If a conventional camera is used, then a digital scanner has to be used to transform the photograph into digital form and save it in the computer. In some laboratory, if the rock sizes permit, rock samples can be just put on a scanner and their surfaces can be transformed into the digital images directly. By this way, the best quality of rock cross-section images can be obtained. Image files saved in computer are a type of electronic files and cannot be processed for numerical modeling directly. Image data have to be extracted from the electronic file for further analyses. Basically, the digital image consists of a rectangular array of image elements or pixels. Each pixel is the intersection area of a horizontal scanning line with a vertical scanning line. These lines all have an equal width h: For the gray level image, at each pixel, the image brightness is sensed and assigned with an integer value that is named as the gray level. For the mostly used 256 gray images and binary images, their gray levels have the integer interval from 0 to 255 and from 0 to 1, respectively. As a result, the digital image can be expressed as a discrete function

f ðx; yÞ in the i and j Cartesian coordinate system: digital image ¼ f ði; jÞ 2 f ð1; 1Þ 6 f ð2; 1Þ 6 ¼6 4 ^

f ð1; 2Þ f ð2; 2Þ ^

? ?

3 f ð1; MÞ f ð2; MÞ 7 7 7 5 ^

f ðN; 1Þ f ðN; 2Þ ? f ðN; MÞ where i from 1 to N; j from 1 to M: ð1Þ For color images, at each pixel, there are three integer values to represent the red, green and blue, so the image data of color image consists of three discrete functions fk ði; jÞ; where k ¼ 1; 2 or 3. An algorithm is developed to extract the colors of each pixel from the image file and to form the discrete functions fk ði; jÞ: These discrete functions fk ði; jÞ are then used in further DIP investigation. Fig. 1 shows a typical example of a color granite cross-section image and an example of the discrete function fk ði; jÞ for the digital image, where the center of the i- and j-coordinate system is located at the top left corner of the image. A small rectangle area of this image on the top right is isolated and enlarged to illustrate the image data structure clearly in Fig. 1(b). For this small area image in Fig. 1(b), the numbers of the scanning lines along both the iaxis and the j-axis are 10. The red, green and blue image data are shown in Figs. 1(c)–(e), respectively. For the original image, the scanning lines along the i and j directions are 200. The actual size of the image is 50.65 mm  50.65 mm. So, the scanning line is 0.253 mm wide. The following two criteria are used for selection of scanning resolution: (a) the pixel should be small enough to cover all the major minerals. In this application, biotites are the smallest major minerals in granite. Each biotite occupies an area varying between 10.75 and 0.26 mm2. Their mean area is 1.24 mm2. The rock sample diameter is 50.65 mm. The resolution of image is 200, which stands for the each pixel (or grid) of 0.253 mm in width and length. The pixel area is 0.064 mm2. So, each main mineral can be represented by more than four pixels. (b) The computation time in the numerical simulation requires the number of pixels (or grids) to be limited. In this paper, the 200  200 grids are used in FLAC. For a homogeneous case with elasticplastic material, the time for achieving equilibrium state on Pentium 866 PC is about 10 min. It is an acceptable computation time. 3.2. Digital image techniques for color image Granite is a common type of igneous rock materials and has been examined by many researchers in rock mechanics in recent years. Studies have proved that

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

943

200 Scanning lines Length=50.65mm

j : 1...10

200 scanning lines Length=50.65mm

O

j

i Biotite Feldspar Quartz

(a)

(b)

(c)

143 139 129 78 50 66 92 111 124 123

141 141 130 91 50 48 55 95 135 147

145 143 135 113 76 57 55 83 154 170

145 152 147 141 128 108 103 138 168 175

137 148 149 146 142 141 149 151 163 177

136 149 148 149 148 152 152 152 154 169

156 151 149 149 153 153 153 156 153 152

177 155 155 155 157 152 153 153 152 151

176 162 161 161 159 151 143 148 147 148

166 162 159 158 158 147 142 146 149 151

(d)

133 129 120 77 53 68 92 111 124 121

130 131 122 89 52 50 58 96 133 144

133 134 126 108 73 57 57 83 153 168

134 143 139 134 120 105 103 136 166 172

125 137 139 137 133 133 142 145 159 175

125 139 137 139 138 142 143 143 148 165

146 141 138 138 141 142 143 147 145 144

170 144 143 143 145 139 142 142 142 143

168 151 150 149 147 139 131 138 137 140

157 151 147 146 146 135 131 136 139 141

(e)

118 116 107 71 54 66 87 104 116 114

115 116 109 81 55 53 59 89 125 136

117 120 112 98 71 58 59 79 143 162

117 130 125 120 108 93 91 124 156 165

110 123 125 122 119 118 127 131 148 166

109 124 122 124 123 126 127 128 134 154

131 126 122 121 123 124 126 130 129 131

158 128 126 126 126 121 125 126 125 127

155 134 132 131 129 120 114 121 120 124

142 133 128 127 127 116 114 119 122 124

Fig. 1. Digital image and discrete functions: (a) original rock image; (b) rectangle gird (local i and j coordinate system); (c) discrete function f1 ði; jÞ (red color); (d) discrete function f2 ði; jÞ (green color); (e) discrete function f3 ði; jÞ (blue color).

inhomogeneity of granite influenced greatly on their mechanical behaviors, physical properties, crack patterns and failure models. There are three main minerals in granite. They are quartz, feldspar and biotite. DIP techniques are used to identify these three individual materials from digital images and to form the internal structure of intact granite masses. In Yue et al. [31], two different digital image techniques are adopted to acquire the actual inhomogeneity of asphalt concrete. They are the region-growing method and the edge detection method. The asphalt concrete is considered to have the

two different individual materials. They are aggregates and matrices. Since rocks are polyphase materials, different digital techniques are used in this paper in the ensuing. The output scanning rock surface images are color images. Color spaces are three-dimensional arrangements of color sensations. Owing to the eyes’ structure of the human beings, all colors are seen as variable combinations of red (R), green (G) and blue (B). Such RGB space is the most frequently used color space for image processing. Most of color cameras, scanners and

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

944

displays are provided with direct RGB signal input and output. As Eq. (1), the basic image data are three scalar functions based on the RGB system. Besides, there are other kinds of color spaces such as the YUV and the YIQ used in television color spaces or the HSI color space in the visual system. The HSI stands for hue, saturation and intensity (brightness) and is an important color space in the human vision system. A normal person cannot make intuitive estimates of the RGB components. But human beings can easily recognize basic attributes of the HSI system. The character H stands for hue and represents the impression related to the dominant wavelength of the color stimulus. In other word, the hue is the domain color perceived by human beings. The character S stands for saturation and corresponds to relative color purity or amount of white mixed with the color. The character I stands for brightness or lightness and is irrelative to colors. The HSI data can be transformed from RGB data. One simplest transformation can be achieved by the following equation [34]: 2 3 ðR  GÞ þ ðR  BÞ 6 7 H ¼ cos1 4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi5; 2 ðR  GÞ2 þ ðR  BÞðG  BÞ

ð2aÞ

where RaG or RaB: If B > G H ¼ 2p  H S ¼1 I¼

3 minðR; G; BÞ; RþBþG

RþBþG : 3

ð2bÞ ð2cÞ

In the above equations, the R; G and B have been transferred from their original interval [0 255] to [0 1]. The H; S and I have their real values between 0 and 1. Digital image techniques are used to acquire actual distribution of rock minerals on a cross-section. The image segmentation method partitions an image into disjoint regions and can be used for this purpose. Our experiences in the image segmentation have shown that the HSI color model is a good color space for identifying individual materials. The advantages are as follows: (a) three components are independent and can be used in color image segmentation process simultaneously; and (b) the H and S components can be sharply recognized by human eyes. Using the color image in Fig. 1(a), we cut a horizontal line (i ¼ 55) in the image and plot the three components along this line below. The values of the three components H; S and I along the blue line are plotted in Figs. 2(b)–(d), respectively. In Fig. 2(a), along the blue line, there exist two concentrated parts of biotite. The first part of the j-coordinate is from 85 to 99 and another part of the j-coordinate is from 149 to 159. Similarly, there exist one concentrated region of quartz, whose the j-coordinate is from 55 to 85. The rest parts along the

blue line are feldspar. From Figs. 2(b)–(d), it can be observed that the variations of the hue and intensity have some connections with the changes of materials. For the H graph, as we mentioned above, H is a color attribute that describes a pure color. Since the biotite is purifying dark, so its H value is as high as that of outside dark regions. But for the feldspar and quartz, their H values are the same. From the I graph, it can be noticed that variations of I value are consistent with the changes of materials. The non-material areas (out spaces of rock circle) have the lowest I values while the biotite materials have the second lowest values. I values of quartz materials are higher than these of biotite. The feldspar materials have the highest I values. Hence, different materials can be distinguished by using thresholds of I values. A histogram is used to display the distribution of the three-color values in the image and is a function to show, for each color value, the number of pixels in the image that have that value. Fig. 3 gives a plot of the histogram for I values associated with the complete image of the square area in Fig. 1. To distinguish individual materials, multi-thresholds can be found in the histogram to segment the image. In Fig. 3, a trough position can be observed between I values 0.1 and 0.2 of the horizontal axis. The high concentration near I value 0.1 indicates that there are many pixels in the image have that value. These pixels are outside rock circle spaces. So, the granite material pixels can be isolated from the entire photo area by setting the threshold as 0.1647. All the pixels whose I values are higher than 0.1647 will be considered as rock materials or components. The results of this preliminary auto-identification are shown in Fig. 4, where the yellow color materials represent the region outside the granite cross-section. It can be observed that the circular disk is detected well. Next, the second threshold can be set as I value equal to 0.2627. The biotite pixels can be automatically detected by defining their I values between 0.1647 and 0.2627. Moreover, the quartz and the feldspar can be detected by setting the third threshold as 0.5059. The quartz pixels have I values between 0.2627 and 0.5059. The feldspar pixels have I values greater than 0.5059. Figs. 5 and 6 show the results of the second and third steps of automatic detections using I values equal to 0.2627 and 0.5059, respectively. The dark color regions represent the biotite materials and the purple color regions represent the quartz materials. By comparing with the original image in Fig. 1(a), one can observe that the individual minerals have been automatically detected correctly and successfully. It is noted that all above threshold values should be selected carefully and can be adjusted via trial-and-error method for different color images of minerals and rock types. The above DIP techniques can produce good quality microstructures for direct-scanning images of granite

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

945

j

o

i=55

(a)

i 0.7 0.6

H

0.5 0.4 0.3 0.2 0.1

(b)

0 0

20

40

60

0

20

40

60

0

20

40

60

80

100

120

140

160

180

200

100

120

140

160

180

200

100

120

140

160

180

200

0.25 0.2

S

0.15 0.1 0.05

(c)

0 80

0.8 0.7 0.6 0.5

I 0.4 0.3 0.2 0.1

(d)

80

Fig. 2. Variation of the discrete function fk ði; jÞ with the j-coordinate at i ¼ 55: (a) a line cut from the original image; (b) discrete function (H) along the blue line; (c) discrete function (S) along the blue line; (d) discrete function (I) along the blue line.

cross-sections. However, since there are different types of images with different dimensions, for example, scanning electron microscopy (SEM) image or X-ray CT image, etc., one may choose appropriate color spaces and adopt appropriate DIP techniques for different images. For example, edge detection method can be also used for acquiring actual material inhomogeneities in color images. It can be considered as an alternative region segmentation method. The aim of this method is to find the boundaries between different minerals or materials and not to detect the regions of

interested. The digital image techniques used in [31] tried to find the maximum of first derivative of gray level and can be applied for color image detection too. Since the color images are viewed as a two-dimensional threevector field, a typical class color edge detector can be defined as linear combinations of the individual detected result for three-color spaces. Different coefficients in the linear combinations can result in different edge detectors. Detailed information of other color edge detectors such as vector gradient operator and directional operator can be found in relevant textbooks such as

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

946

1000 900 800 700 600 500 400 300 200 100 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 6. Image segmentation result with the third threshold of 0.5059 for I value.

Fig. 3. Histogram of I value in Fig. 1(a).

Fig. 4. Image segmentation result with the first threshold of 0.1647 for I value.

Fig. 5. Image segmentation result with the second threshold of 0.2627 for I value.

the book by Plataniotis and Venetsanopoulos [35] and others. In brief, a careful examination of the individual mineral color characteristics is always a need for quality detection and identification of actual internal structures of rock materials. Due to noise influence or dirty on scanning plate during the image digitalization process, the detection

Fig. 7. Revised microstructure of granite.

result may not be 100% correct. Some digital image software packages such as Painter or Photoshop are used here to do some minor refinements manually by comparing the automatic results with the original color images via human eye view. The degree of the manual detection depends on the image contrasts and the quality of automatic results. The final result of the identified internal distribution of the three minerals in the granite cross-section is shown in Fig. 7. In summary, the number of the total pixels is equal to 40,000. The numbers of the pixels occupied by feldspar, quartz and biotite are found to be 22,267, 8257 and 744, respectively. Since each pixel occupies a plan area of 0.253  0.253 mm2, so the real areas occupied by feldspar, quartz and biotite are 1428.1, 529.5 and 47.7 mm2, respectively.

4. Vector transformation Fig. 7 shows the processed image containing the actual spatial distribution of three main minerals on the granite cross-section. In order to carry out the numerical analysis, the image microstructure must be transformed

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957 j=1:12

o

y=0~3.04mm

h

j

947

y

o

nodes i=1:10

x=0~2.53mm

h

i

x pixel

(a)

(b)

grid

Fig. 8. An example of transformation from pixel to mesh (grid): (a) image and (b) mesh (grid).

into the vector microstructure. We have developed two methods to carry out the transformation. One method is called mesh free vector transformation [33] and the other is called mesh-generation vector transformation [31,32]. 4.1. Mesh free vector transformation In conventional numerical methods such as finite element method (FEM) and FDM, the rock material domains have to be discretized into many connected small regions. Since a digital image consists of a rectangle array of pixels. Each pixel occupies a square area on the rock cross-section. So, each pixel can be considered as either a finite element or a finite difference grid. For each square grid, its four corner coordinates can be calculated by the corresponding positions in the digital image plus the scale distance h: Fig. 8 shows a comparison of the digital image pixels and the vector square element data. The scale distance h can be calculated as follows:

h¼

image consists of 1000 horizontal scanning lines and 1000 vertical scanning lines, then the total number of square element meshes or grids will be equal to 1,000,000. Mechanical computations with such large number of finite elements would require some long computational time and storage space. Such requirements would not be a problem with the fast development of more and more powerful personal computer hardware and software techniques. 4.2. Mesh-generation vector transformation In order to reduce the mesh sizes and to optimize and customize meshes or girds for specific rock mechanics problems, a so-called geometry vectorization transform algorithm has been developed in [31,32]. In this algorithm, the image microstructure data are transformed into vector polygon data. The principles of this transformation are summarized as follows: (1) The actual microstructures of geomaterials are represented by the image containing the boundaries between different materials.

length of rock cross-section in vertical or horizontal direction : total number of pixels on the cross-section image along the selected direction

By using this equation, we can transform the image data into a vector square mesh data easily. The mesh free vector has the following advantages: (1) The process is simply and avoids of complicated and time-consuming mesh-generation steps. (2) All the grids or meshes are square and have exactly the same sizes. (3) The square mesh data can be incorporated into many mechanical software packages for numerical modeling in rock and soil mechanics. The main disadvantage of the mesh free vector data is related to the size of the finite square meshes. If a digital

ð3Þ

(2) The shapes of boundaries can be replaced by polygons if polygons have sufficient sidelines. (3) A geometry transformation algorithm is developed to transform the cluster of close pixels into polygons that have a very small data size. (4) The output vector data can be recorded as the coordinate values of polygon vertices. (5) The automatic mesh-generation algorithms are used for mesh generation. (6) The original image sizes can be very large. Such original images with large sizes can be reduced to small elements or grids for efficient computations in subsequent numerical predictions.

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

948

Original image

Microstructure detection Digital image processing Image segmentation method

Edge detection method

Output revised microstructure

Square mesh vector data

Polygon vector data

Vector data transformation

Fig. 10. Generated grid in FLAC. Mesh generation

Output mesh (grid) data

Incorporated with FEM or FDM software for numerical analyses

Fig. 9. Flow chart of procedures in a digital image modeling method.

The square meshes (grids) or generated meshes (grids) can be incorporated into the existing FEM or FDM for numerical predictions. In the ensuing, the mesh free vector data and the commercial software package FLAC are used to show the digital image-based numerical modeling method for prediction of inhomogeneous rock failure under static loading. Fig. 9 shows a flow chart of the examination procedures for the proposed method.

5. Numerical evaluations FLAC is a FDM-based commercial software package. It stands for the fast Lagrangian analysis of continua. FLAC implements an explicit time marching scheme to solve Newton’s second law to describe material deformation and embodies some basic constitutive models for simulation of highly non-linear, irreversible response representative of geomaterials. FLAC contains a built-in programming language FISH, which enables users to define new variables and functions. This provides an easy way to establish FLAC models for any rock geometries and solve complex problems in rock mechanics. In the present study, a FISH program is developed for establishment of actual material microstructure according to the vector data. Firstly, grids are generated for the whole geometry. According to the original image size, there are 200  200 grids in the whole model. As shown in Fig. 7, there are four different material regions associated with the rectangle grids. The first material region is the space outside the rock cross-section and is represented by

Fig. 11. Inhomogeneity in the FLAC model.

yellow color. The second, third and fourth material regions are those occupied by biotite, quartz and feldspar minerals, respectively. They are represented by dark, pink and blue colors, respectively. The yellow non-rock material region can be represented as the ‘‘null’’ material model in FLAC. This approach has the same affect as to empty the yellow region. For the three mineral regions, the FISH can be used to assign each grid specific material properties according to its mineral type. They were represented by different colors in the image microstructure. Fig. 10 shows the generated circular grids in FLAC. Fig. 11 shows the introduced inhomogeneities in the FLAC grid. In the ensuring, linear and non-linear analyses are presented to show the response of the rock cross-section under the Brazilian split tension testing with or without consideration of the material non-homogeneities. The Brazilian test is a laboratory test to measure the tensile strength of the rock indirectly (Fig. 12a). A pair of vertical compressive loads is applied over the two thin steel strips on the circumference of a circular cylindrical specimen. One of the two thin strips is located at the circumference crown (Point A in Fig. 12a) and the other is located at the circumference bottom (Point B in Fig. 12a). Actual distribution of the load F at the

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957 20 grids

σγ

y

Feldspar

Quartz

We further define the non-dimensional normal stresses s0y and s0r as follows: sy ð6aÞ s0y ¼ ; s0 s0r ¼

σθ x 2α

α=7.5°

(a)

(b)

o Center of the disc

Fig. 12. A sketch diagram illustrating the conventional Brazilian split tensile test: (a) sketch of Brazilian test and (b) details of load applied on square grids.

Point A is shown in Fig. 12b. From Fig. 12b, it can be observed that the load F is applied at the node points of 20 horizontal square elements on the boundary of the sample. Fourteen of the 20 square elements are occupied by quartz and the other six are occupied by feldspar. The actual distribution of the load F at the Point B is similar to those in Fig. 12b. Assuming that the material is homogeneous, isotropic and linear elastic, the two normal stresses along the vertical loading diameter AB can be estimated using the following solution (Hondros, [36]): ( F ½1  ðr=r0 Þ2 sin 2a sy ¼  pr0 ta 1  2ðr=r0 Þ2 cos 2a þ ðr=r0 Þ4 " #) 2 1 1 þ ðr=r0 Þ tan tan a ; ð4aÞ 1  ðr=r0 Þ2 ( F ½1  ðr=r0 Þ2 sin 2a sr ¼ pr0 ta 1  2ðr=r0 Þ2 cos 2a þ ðr=r0 Þ4 " #) 2 1 1 þ ðr=r0 Þ þtan tan a ; ð4bÞ 1  ðr=r0 Þ2 where sy is the normal stress component along the hoop direction, sr the normal stress component along the radial direction, F the applied force, r0 the radius of the circular cylindrical specimen, t the thickness (or height) of the circular cylindrical specimen, a the arc angle over which the F is uniformly applied to and r the radial distance from the center of the circular cylindrical specimen. For r ¼ 0; the sy at the center is given by

 F sin 2a F 1 E ¼ s0 ; sy ¼  ð5Þ pr0 t a pr0 t where s0 ¼ F =pr0 t:

949

sr : s0

ð6bÞ

The Brazilian test is applicable when the primary fracture initiates from the center of the disk. Its tensile strength can be estimated by calculating the horizontal stress sy that corresponds to the applied peak load. Eq. (6) can be used to estimate the tensile strength indirectly from the applied peak load. 5.1. Linear analyses A linear elastic and isotropic constitutive model is used to represent the stress–strain relationship for the granite. The basic parameters for this model are elastic modulus and Poisson’s ratio. FLAC, during the calculations, uses the bulk moduli (K) and shear moduli (G). These two parameters can be calculated by E G¼ ; ð7aÞ 2ð1 þ vÞ K¼

E : 3ð1  2vÞ

ð7bÞ

Two cases are analyzed. One is to assume all the minerals have the same properties, i.e. conventional homogeneous rock material case. The other is to assume the three rock minerals to have different material properties. Using the two cases, effects of the granite internal inhomogeneity on the mechanical response can be examined. 5.1.1. Homogeneous cases All the grids in the FLAC model in Fig. 10 have the same material properties. The modulus E is equal to 50 GPa and the Poisson’s ratio n is equal to 0.2. The analytical results of the non-dimensional normal stresses s0y and s0r from Eq. (4) along the loading diameter AB are plotted in Figs. 13 and 14, respectively. The numerical results from FLAC for the two non-dimensional normal stresses are also plotted in Figs. 13 and 14. It is evident that the numerical results from the present DIP approach agree well with the analytical results. The above simply numerical comparison also revealed an important aspect of solution convergence in using FLAC. In FLAC, the term maximum unbalanced force is used to determine whether or not the computation has reached an equilibrium or steady plastic flow stage. As an explicit code, the solution of FLAC problems requires a number of computational steps to approach a small unbalanced force. When FLAC is used for examining actual problems in geotechnical and mining

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

950

engineering fields with dimensions between 10 and 1000 m, the unbalanced force has a limit of 100 N. When FLAC is used for the present problem of very small sizes (50 mm), the corresponding limit of the maximum unbalanced force has to be changed to a very small value. Fig. 15 shows the variation of the maximum unbalanced force with the calculation steps. The maximum unbalanced force approaches to a small value of 0.6 N after the 5000th calculation steps.

16 σ θ / σ0 (Compression is positive)

14

Analytical solution FLAC result

12 10 8 6 4 2 0 -2 0.0

0.5 1.0 1.5 Normalized coordinate (y/r0)

2.0

Fig. 13. Comparison of the exact solution with the FLAC result for the non-dimensional hoop normal stress sy =s0 along vertical loading diameter.

σr / σ0 (Compression is positive)

5

Analytical solution FLAC result

4 3

5.1.2. Inhomogeneous cases It is known that quartz is harder than feldspar and feldspar is harder than biotite. Based on experience, quartz and feldspar are assumed to have their elastic modulus values equal to 90 and 70 GPa, respectively. The biotite is assumed to have its elastic modulus value equal to 40, 20 or 10 GPa. Figs. 16 and 17 show the variations of the normal stresses s0y and s0r along the vertical loading diameter AB, respectively. The mineral distribution along the diameter AB is also illustrated above the abscissa in Figs. 16 and 17, where the upper parts represent the higher elastic modulus values. From Figs. 16 and 17, the following observations can be made:

2 1 0 0.0

0.5

1.0

1.5

2.0

Normalized coordinate (y/r0) Fig. 14. Comparison of the exact solution with the FLAC result for the non-dimensional radial normal stress sr =s0 along vertical loading diameter.

(1) Some local fluctuation points of normal stress curves are present. Their locations have strong correlations with the spatial distribution of the three minerals along the vertical diameter AB. The curves change greatly at the boundaries of the three minerals. (2) Variations of the two normal stresses are more significant along the diameter AB as the differences in the mineral modulus values increase.

1.60E+03

Maximum unbalanced force

1.40E+03

1.20E+03

1.00E+03

8.00E+02

6.00E+02

4.00E+02

2.00E+02

0.00E+00 0.00E+00 1.00E+03 2.00E+03 3.00E+03 4.00E+03 5.00E+03 6.00E+03 7.00E+03 8.00E+03 9.00E+03 1.00E+04

Calculation steps Fig. 15. Variation of maximum unbalanced force.

ARTICLE IN PRESS

0.5

1.0

1.5

2.0

11.0 7.0

7.0 14.0

Fig. 16. Variations of the non-dimensional hoop normal stress sy =s0 along the vertical loading diameter.

14.0 11.7

Tensile strength (MPa)

Normalized coordinate y/r0

5.2.1. Homogeneous cases The grid in Fig. 10 is further used. The three minerals are assumed to have the same properties. The elastic modulus, the Poisson’s ratio, the cohesion, the friction angle and the tension strength for the homogeneous rock surface are given in Table 1. Three cases are to be

40 30 40.0 25.0 70.0 40.0 Feldspar Biotite

0.2 0.2

30 60 25.0 50.0 0.2 0.2

60 45 50.0 42.3

40.0 90.0

In the ensuing, the classical Mohr–Coulomb plastic and tensile failure criteria are used to model the failure behavior of the three main granite minerals. The corresponding mechanical parameters for this model are the elastic modulus, the Poisson’s ratio, the cohesion, the internal friction angle and the tensile strength. The thickness of the sample t=25.25 mm is used in the computational analysis. Similarly, the above homogeneous and inhomogeneous approaches are also adopted for the failure predictions and comparisons.

Inhomogeneous material (actual distribution of minerals on the granite crosssection)

5.2. Failure predictions

H3 IH

(3) The maximum hoop tensile stress does not occur at the geometrical center of the circular cross-section. As difference in the modulus values increases, the maximum hoop tensile stress can reach the maximum tensile strength at some local locations. Cracks could then be initiated at some local high hoop tensile stress regions instead of the geometrical center in the homogeneous case.

Table 1 Material properties for the three homogeneous cases and the one inhomogeneous case

Fig. 17. Variations of the non-dimensional radial normal stress sr =s0 along the vertical loading diameter.

Homogeneous material

1.5

H1 H2

1.0

Cohesion (MPa)

0.5

Normalized coordinate (y/ro)

Poisson’s ratio

0.0

Material distribution along the loading diameter 2.0

Elastic modulus (GPa)

0

0.2 0.2

1

90.0 74.6

2

Description

p σr / σ0

3

Quartz (100%) Weighted average of case IH () Biotite (100%) Quartz

Case 1 Case 2 Case 3

4

Case no.

5

Friction angle ( )

6

Note: =Properties equal to the weighted averages of the properties of the inhomogeneous cross-section in Case IH by the areas occupied by quartz (26%), feldspar (71%) and biotite (3%).

0.0

951

47.6 43.3

Material distribution along the loading diameter

43.3 28.9

Case 1 Case 2 Case 3

28.9 42.3

14 12 10 8 6 4 2 0 -2 -4

Tensile strength from Mohr–Column failure criterion=c=tanðfÞ (MPa)

σθ / σ0

S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

952

(a)

(b)

(c)

(d)

Stage 1 (Point 1)

Stage 1 (Point 1) Point 4

Stage 1 (Point 1)

24.615

Force (kN)

Force (kN)

29.020 29.014 29.008

Point 3 Stage 2

29.002

(e)

Stage 3 (Point 2)

Point 3

24.597

Stage 2

24.590 24.584

Stage 3 (Point 2) Stage 4

14.499 14.496

Point 4 Stage 2

14.493

Point 3 14.490

Stage 3 (Point 2)

14.487

Stage 4

21.06 21.0 7 21.08 21.09 21.10 21.11

Displacement (µm)

24.603

24.578 24.572

28.996 28.990

Point 4

24.609

Force (kN)

29.026

(f)

21.36 21.40 21.45 21.49 21.55 21.59

Displacement (µm)

(g)

23.99

24.01

Stage 4

24.03

Displacement (µm)

Fig. 18. Failure propagation process in homogeneous rock material: (a) Stage 1; (b) Stage 2; (c) Stage 3; (d) Stage 4; (e) force–displacement curve for H1 case; (f) force–displacement curve for H2 case; (g) force–displacement curve for H3 case.

examined. Case H1 assumes that all the three minerals have the property values as exactly the same as those of quartz. Case H3 assumes that all the three minerals have the property values as exactly the same as those of biotite. Case H2 assumes that all the three minerals have the weighted average values of those of quartz, feldspar and biotite. In addition, the assumed property values of the three minerals are listed under the case IH for the actual inhomogeneous case in this study. From Table 1, it can be observed that the tensile strength values associated with the Mohr–Coulomb plastic failure criterion are much greater than those associated with the tensile failure criterion. So, the expected failure mode will be in tension. Numerical representation of a failure such as cracking can be represented by a collection of failed square elements. When the stress status of a square element violates the failure criteria, then element fails and is marked to distinguish it from the un-failed elements. As loading increases, the number of failed elements increases, and the failed elements may connect to from

elongated regions that correspond to macroscopic fractures. Two common types of failure criteria are the shear failure and the tensile failure for Mohr–Coulomb plastic material model. In the numerical computation, the loading is increased gradually and the failed elements are recorded. Figs. 18(a)–(d) show the process of crack initiation and propagation, where the purple elements have violated the tensile failure criterion. The three cases H1, H2 and H3 have the same failure process pattern as shown in Fig. 18(a)–(d). However, their actual stress and displacement values are different. Figs. 18(e)–(g) show the relations between the applied displacement and the associated applied load F at the loading Point A for the three cases H1, H2 and H3 immediately before the peak and after the tensile failure. In particular, the four stages in Figs. 18(a)–(d) correspond the Points 1, 2, 3 and 4 in Figs. 18(e)–(g), respectively. From Figs. 18(a)–(d), it can be observed that the crack initiates at the center of the granite circular crosssection. The coordinates of the first failed grid (i; j) are

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

equal to (101, 100). The crack propagates almost vertically to the two loading positions simultaneously. It is noted that the failed elements in lower semi-circular cross-section form a vertical seam and the failed elements in the upper semi-circular cross-section also form a vertical seam. The lower vertical seam and the upper vertical seam are relocated by one grid. The lower vertical seam has its horizontal coordinate 100 in the (i; j) coordinate system. It can be considered as the vertical center of the circular cross-section because the grid consists of 200 vertical columns and 200 horizontal rows of the same size square elements. The crack propagation process is symmetry to the central horizontal axis. All the failed square elements violate only the tensile failure criterion. So, the two vertical seams can be considered to form a tensile crack line along the central vertical position. Such predicted failure process can be considered to agree well with the homogeneous material assumption. Furthermore, from Figs. 18(e)–(g), it can be observed that the three cases H1, H2 and H3 have the similar load and displacement curves at the failure and have substantially different values. Further details will be analyzed in the ensuing. 5.2.2. Inhomogeneous cases In the ensuing, the quartz, feldspar and biotite are assumed to have different mechanical properties. Table 1 lists the values of the mechanical properties for the three minerals under the case No. IH. Figs. 19(a)–(h) show the predicted failure behavior for the inhomogeneous rock under the Brazilian split loading condition. It is evident that the failed elements do not appear at the center of the circular cross-section. The coordinates (i; j) of the first failed element are (104, 109). This square element is occupied by biotite. In Fig. 11, it can be observed that two biotite regions are located near the cross-section center. Elements that have maximum tensile stresses in the elastic range are near the disk center. When its tensile stress violates the tensile strength, the element will fail and can be considered as a failed element. The first failed biotite element is among the biotite region that is larger than the other biotite region and is closer to the geometrical center. It would experience a higher tensile stress in the elastic range and reach to the biotite tensile strength earlier. This first failed element would continue to affect adjacent square elements and cause more elements experiencing failure. Figs. 19(a)–(g) show the crack initiation and propagation process in the granite cross-section. All failed elements violated only the tensile failure criteria. Therefore, there represent tensile cracking failure. The failed square elements form a near vertical cracking seam toward to the two loading ends. The numerical results have shown that the tensile stress and strengths control the failure model in the Brazilian test for the stiff

953

granite. The actual spatial distribution of the three main minerals dominates actual crack patterns in the granite cross-section. Fig. 19(h) shows the relationship between the applied displacement and the associated load F at the Point A during immediate before and after the tensile failure. The Stages 1–7 in Figs. 19(a)–(g) are also shown on the curve. 5.2.3. Analysis and discussions Fig. 20 further shows a comparison among the relations between the applied displacement and the associated total force F at the loading Point A (or B) from the initial loading to the post-peak behavior for the three homogeneous cases and the one inhomogeneous case. Table 2 gives a summary of the results of the applied forces at the failure Stages 1–4 and the analytical prediction of the peak tensile force for the three homogeneous cases, where the given tensile strength for each of Case Nos. H1, H2 and H3 are used in Eq. (6). Table 3 gives a comparison between the given tensile strength values and the predicted tensile strength values. From the results presented in Figs. 18–20, Tables 2 and 3, the following important observations can be made: (a) For the four cases, the tensile failure process is a quick process since the post-peak displacements are extremely small (see Figs. 18(e)–(g) and 19(h)). (b) For three homogeneous cases H1, H2 and H3, the numerically predicted applied forces at the split tensile failure are about 2.70%, 4.25% and 2.68% greater than those analytical predicted results, respectively. The numerical predicted tensile strength values are about 2.73%, 4.21% and 2.61% greater than the given values, respectively. Furthermore, the corresponding failure displacements at the loading point are higher for those cases with smaller modulus. These results are within the expectation, which also gives confidence in the numerical computation. (c) For the inhomogeneous case IH, its load and displacement curve is close similar to that of the case H2 before the tensile failure at the peak (Fig. 20). This result is also expected since the case H2 has the weighted average value of the modulus values associated with the case IH. (d) However, the numerically predicted failure loads for the case IH are much less than those homogeneous cases H1, H2 and H3. Since the cases H1 and H2 have higher tensile strength values, this result is expected. (e) But, all the tensile strength values for the three minerals in the case IH are larger or equal to that in the case H3. The predicted peak force for the case

ARTICLE IN PRESS 954

S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

Fig. 19. Failure propagation process in inhomogeneous rock material (purple grids represent tension failed elements, other colors represent different materials): (a) Stage 1; (b) Stage 2; (c) Stage 3; (d) Stage 4; (e) Stage 5; (f) Stage 6; (g) Stage 7; (h) crack process in the force–displacement post peak curve.

IH is about 25% less than the peak force for the case H3. Correspondingly, the numerical predicted average tensile strength for the case IH is about 158%, 103% and 29% less than those of the three minerals quartz, feldspar and biotite, respectively.

This result is unexpected and of great interesting if it can be further confirmed. The above results have indicated that the predicted failures for the four cases are tensile failure, as expected.

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

A majority of the results are within the expectations. The expected low peak force result for the inhomogeneous case IH may be explained by the fact that the stress field in the rock cross-section is highly affected by the local material differences (see, Figs. 14, 16 and 17). Tensile failures can be initiated at local weak points due to material inhomogeneity induced stress concentrations under a much lower applied load. Such initial failure can be quickly propagated to cause the global failure in tension.

30 H1

Force (kN)

25 20 15

H2

10

H1 H2 H3 IH

H3

5 IH

0 0

5

10 15 20 Displacement (um)

25

Fig. 20. Force and displacement response for different cases.

955

6. Summary and conclusions This paper has presented a digital image processingbased numerical modeling method for two-dimensional failure prediction of rocks by taking into account the actual material non-homogeneity. The method is a combination of digital color image processing techniques, vector data transformation and conventional computational mechanics methods. Using granite as an example, the paper has shown that this numerical method can be used to carry out the prediction of rock failure behavior by taking into account the actual mineral distribution in the rock mass with efficiency and accuracy. A finite difference software FLAC is adopted in this paper to show the flexibility of the proposed method. The conventional Mohr–Coulomb plastic model is used for examining the failure behavior of the granite. The numerical results presented in the paper have demonstrated that the material non-homogeneity has significant effects on the tensile stress distribution along the vertical loading diameter for the classical Brazilian split test in determination of an average tensile strength. The numerical results have also shown that the maximum tensile stress cannot occur at the center of Brazilian circular disk if rock inhomogeneity is considered. The crack initiation and propagation can be

Table 2 Results of the total applied failure forces from the present numerical method for the four cases and the conventional analytical method for homogeneous rocks Case No.

H1 H2 H3 IH

Numerical prediction of the failure forces at

Analytical prediction

Difference between A and B

Point 1 (kN)

Point 2 (kN)

Point 3 (kN)

Point 4 (kN)

Average (A) (kN)

Given tensile Predicted failure strength force (B) (kN) (MPa)

Difference (A–B) (kN)

Percentage (A–B)/B100 (%)

29.026 24.620 14.499 10.976

28.996 24.585 14.487 10.926

29.009 24.592 14.489 10.941

29.029 24.610 14.492 10.966

29.015 24.602 14.492 10.952

14.0 11.7 7.0 ?

0.762 1.001 3.378 ?

2.70 4.25 2.68 ?

28.253 23.600 14.114 ?

Note: =from Eq. (5) where r0=0.025325 m, a=5.71 and t=0.025025 m.

Table 3 Results of the tensile strength prediction from the present numerical method Case no.

H1 H2 H3 IH

Given tensile strength (B) (MPa)

Numerical predicted results of tensile strengths at failure

Difference between A and B

Point 1 (MPa) Point 2 (MPa) Point 3 (MPa) Point 4 (MPa) Average (A) (MPa)

Difference (A–B) (MPa)

Percentage (A–B)/ B100%

0.3825 0.4925 0.1825 8.5725 5.5725 1.5725

2.73 4.21 2.61 157.95 102.67 28.97

14.0 14.39 11.7 12.20 7.0 7.19 (Quartz) 14.0 (Feldspar)11.0 5.44 (Biotite) 7.0

14.37 12.18 7.18

14.38 12.19 7.18

14.39 12.20 7.18

14.3825 12.1925 7.1825

5.42

5.42

5.43

5.4275

ARTICLE IN PRESS 956

S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957

recorded at each step of the numerical calculations. It has shown that the weakest mineral biotite controls the crack initiation and failure pattern for a granite crosssection subject to the Brazilian test loading. The results have further shown that a majority of the numerically predicted failures compare well with the analytical results and are within the expectations. The numerically predicted failure load and tensile strength values for the inhomogeneous case are much lower than expectation. If it can be confirmed, this finding can be used as an additional important justification to take into account actual material inhomogeneities in failure prediction of granitic rocks under loading. The method presented in this paper can offer an efficient approach to extend the conventional and powerful numerical methods to mechanical analysis of actual non-homogeneous materials subject to loading. Results of such mechanical analyses can be used to examine the applicability and limitations associated with the conventional macro-mechanical approaches for rock mechanics. Furthermore, although this paper deals with mainly on two-dimensional problems in rock mechanics, it is believed that the approach and techniques presented in the paper can be further applied to the examination of three-dimensional problems in rock mechanics via a natural extension of the square elements in twodimension to the regular cubic elements in threedimension.

Acknowledgements The authors would like to thank financial supports from the Research Grants Council of Hong Kong SAR Government. One of the authors Mr. S. Chen would also like to thank the University of Hong Kong for postgraduate scholarship.

References [1] Song J, Kim K. Micromechanical modeling of the dynamic fracture process during rock blasting. Int J Rock Mech Min Sci 1996;33:387–94. [2] Tang CA, Liu H, Lee PKK, Tsui Y, Tham LG. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—Part I: effect of heterogeneity. Int J Rock Mech Min Sci 2000;37:555–69. [3] Blair SC, Cook NGW. Analysis of compressive fracture in rock using statistical techniques: Part I. A non-linear rule-based model. Int J Rock Mech Min Sci 1998;35:837–48. [4] Cundall PA. Formulation of a three-dimensional distinct element model—Part I: a scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int J Rock Mech Min Sci 1988;25:107–16. [5] Tang CA, Tham LG, Lee PKK, Tsui Y, Liu H. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—Part II: constraint, slenderness and size effect. Int J Rock Mech Min Sci 2000;37:571–83.

[6] Li L. Microscopic study and numerical simulation of the failure process of granite. PhD thesis, Department of Civil Engineering, University of Hong Kong, Hong Kong, 2001. p. 175. [7] Tham LG, Li L, Tsui Y, Lee PKK. A replica method for observing microcracks on rock surfaces. Int J Rock Mech Min Sci 2003;40:785–94. [8] Fang Z, Harrison JP. Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks. Int J Rock Mech Min Sci 2002;39:443–57. [9] Fang Z, Harrison JP. Application of a local degradation model to the analysis of brittle fracture of laboratory scale rock specimens under triaxial conditions. Int J Rock Mech Min Sci 2002;39: 459–76. [10] Schlangen E, Van Mier JGM. Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cem Concr Compos 1992;14:105–18. [11] Li LH, Bai YL, Xia MF, Ke FJ, Yin XC. Damage localizations as a possible precursor of earthquake rupture. Pure Appl Geophys 2000;157:1929–43. [12] Napier JAL, Dede T. A comparison between random mesh schemes and explicit growth rules for rock fracture simulation. Int J Rock Mech Min Sci 1997;34:356. [13] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979;29:47–65. [14] Hazzard JF, Young RP. Simulating acoustic emissions in bondedparticle models of rock. Int J Rock Mech Min Sci 2000;37: 867–72. [15] Fakhimi A, Carvalho F, Ishida T, Labuz JF. Simulation of failure around a circular opening in rock. Int J Rock Mech Min Sci 2002;39:507–15. [16] Konietzky H. Numerical modeling in micromechanics via particle methods. Proceedings of the First International PFC Symposium, Gelsenkirchen, Germany, 6–8 November 2002, Lisse, 2003. p. 321. [17] Itasca. Fast Lagrangian analysis of continua, version 3.3. Minnesota: Itasca Consulting Group, Inc.; 1995. [18] Sepher K, Svec OJ, Yue ZQ, Hussein HM. Finite element modelling of asphalt concrete microstructure. In: Aliabadi MH, Carpinteri A, Kalisky S, Cartwright DJ, editors. Localized damage III: computer-aided assessment and control. Southampton, Boston: Computational Mechanics Publications; 1994. p. 225–32. [19] Frost JD, Wright JR. Proceedings of the EF/NSF Conference on Digital Image Processing: Techniques and Applications in Civil Engineering. Hawaii: A.S.C.E.; 1993. p. 301. [20] Lee H, Chou E. Survey of image processing applications in civil engineering. Proceedings of the EF/NSF Conference on Digital Image Processing: Techniques and Applications in Civil Engineering, 1993. p. 203–10. [21] Chermant JL. Why automatic image analysis? An introduction to this issue. Cem Concr Compos 2001;23:127–31. [22] Chermant JL, Chermant L, Coster M, Dequiedt AS, Redon C. Some fields of applications of automatic image analysis in civil engineering. Cem Concr Compos 2001;23:157–69. [23] Hadjigeorgiou J, Lemy F, Cote P, Maldague X. An evaluation of image analysis algorithms for constructing discontinuity trace maps. Rock Mech Rock Eng 2003;36(2):163–79. [24] Yue ZQ, Bekking W, Morin, I. Application of digital image processing to quantitative study of asphalt concrete microstructure. Transportation Research Record 1492, Transportation Research Board, National Research Council, Washington, DC, 1995. p. 53–60. [25] Yue ZQ, Morin I. Digital image processing for aggregate orientation in asphalt concrete mixtures. Can J Civ Eng 1996; 23:479–89.

ARTICLE IN PRESS S. Chen et al. / International Journal of Rock Mechanics & Mining Sciences 41 (2004) 939–957 [26] Kwan AKH, Mora CF, Chan HC. Particle shape analysis of coarse aggregate using digital image processing. Cem Concr Res 1999;29:1403–10. [27] Mora CF, Kwan AKH, Chan HC. Particle size distribution analysis of coarse aggregate using digital image processing. Cem Concr Res 1998;28:921–32. [28] Reid TR, Harrison JP. A semi-automated methodology for discontinuity trace detection in digital images of rock mass exposures. Int J Rock Mech Min Sci 2000;37:1073–89. [29] Hollister SJ, Kikuchi N. Homogenization theory and digital imaging: a basis for studying the mechanics and design principles of bone tissue. Biotechnol Bioeng 1994;43:586–96. [30] Terada K, Kyoya T, Kazama M, Lee K, Oyang L. Image-based modeling and analysis of microstructures for two-scale problems in geomechanics. Int J Numer Anal Methods Geomech 2002; 26:273–97. [31] Yue ZQ, Chen S, Tham LG. Finite element modeling of geomaterials using digital image processing. Comput Geotech 2003;30:375–97.

957

[32] Yue ZQ, Chen S, Tham LG. Digital image processing based finite element method for rock mechanics. In: Lin YM, Tang CA, Feng XT, Wang SH, editors. New development in rock mechanics and rock engineering. The Proceedings of the Second International Conference, PR China; 2002. p. 609–15 [Supplement]. [33] Yue ZQ, Chen S, Tham LG. Seepage analysis in inhomogeneous geomaterials using digital image processing based finite element method. Proceedings of the 12th Panamerican Conference for Soil Mechanics and Geotechnical Engineering and the 39th US Rock Mechanics Symposium, Soil and Rock America 2003. Boston, USA: MIT; June 2003. p. 1297–302. [34] Haralick RM, Shapiro LG. Glossary of computer vision terms. Pattern Recognition 1991;24:69–93. [35] Plataniotis KN, Venetsanopoulos AN. Color image processing and applications. Berlin: Springer; 2002. p. 355. [36] Hondros G. The evaluation of Poisson’s ratio and the modulus of materials of a low tensile resistance by the Brazilian (indirect tensile) test with a particular reference to concrete. Aust J Appl Sci 1959;10:243–68.