Rock Mech. Rock Engng. (2007) 40 (2), 145–168
DOI 10.1007/s00603-006-0105-8
Printed in The Netherlands
Digital Image Based Approach for Three-Dimensional Mechanical Analysis of Heterogeneous Rocks
By
S. Chen, Z. Q. Yue, and L. G. Tham
Department of Civil Engineering, The University of Hong Kong,
Hong Kong, P.R. China
Received June 30, 2005; accepted May 22, 2006
Published online September 4, 2006 # Springer-Verlag 2006
Summary
This paper presents a digital image based approach for three-dimensional (3-D) numerical simulation and failure analysis of rocks by taking into account the actual 3-D heterogeneity. Digital image techniques are adopted to extract two-dimensional (2-D) material heterogeneity from material surface images. The 2-D image mesostructures are further extrapolated to 3-D cuboid mesostructures by assuming the material surface as a representation of the inner material heterogeneity within a very small depth. The iterative milling and scanning system is set up to generate the 3-D rock mesostructures. A Hong Kong granite specimen is used as an example to demonstrate the procedure of 3-D mesostructure establishment. The mechanical responses and failure process under the conventional Brazilian tensile test condition are examined through numerical analyses. The stress distribution, crack propagation process and failure model of heterogeneous material cases are simulated with a finite difference software. The numerical results indicate that material heterogeneity plays an important role in determining the failure behavior of rocks under external loading.
Keywords: Heterogeneous materials, rock, digital image processing, numerical simulation, finite difference method, and failure prediction.
1. Introduction
Rock is a natural material and is widely encountered in engineering and scientific activities. It is essential for engineers to accurately predict the mechanical behavior of rock under external loading. However, even in laboratory testing environment, rock samples with the same dimensions and loading conditions can have different failure patterns. It is believed that the heterogeneity of rocks is one of the crucial parameters that controls rock failure models and crack patterns. Importance of such meso-material heterogeneity has been well recognized by researchers in rock mechanics.
During the last decades, some researchers (Cundall and Strack, 1979; Tang et al., 2000; Fang and Harrison, 2002a, b) have attempted to examine the material properties by taking the heterogeneity into consideration. Due to measurement difficulty, statistical models were widely used to analyze the heterogeneity of rock. Such models utilized statistical tools or random generators to simulate rock mesostructures via some theoretical material heterogeneity parameters. The simulated rock mesostructures were considered to be statistically equivalent to that of the actual rocks.
Using such statistical models, various numerical analyses based on finite element method (Tang et al., 2000) or finite difference method (Fang and Harrison, 2002a, b) have been developed to examine mechanical behavior, crack initiation, crack propagation and crack coalescence of rocks under loading. It is therefore necessary to compare the results of the simulated mesostructures with those of the actual rock mesostructures. The accuracy of the statistics-based models can be enhanced if actual rock heterogeneous properties can be taken into account in the numerical modeling. Direct measurements of actual mesostructures become important since they can give an actual exploration of the rock material properties and behavior.
In recent years, digital image processing techniques have become a powerful tool for measuring quantitative material data and establishing material spatial distribution (Frost and Wright, 1993; Frost and McNeil, 1998). In particular, Yue and his colleagues (Yue et al., 2002, 2003a, b; Chen et al., 2004a, b) have used digital image techniques and developed numerical modeling methods for mechanical analysis of heterogeneous geomaterials. The methods rely on digital image processing techniques to establish actual spatial distribution of different particles and components in geomaterials. The digitalized actual mesostructures of geomaterials can then be integrated with conventional numerical methods for mechanical analysis and failure prediction.
Literature review has shown that a majority of the research developments has focused on the understanding of various problems in heterogeneous geomaterials under 2-D conditions. Mechanical problems in geomaterials are generally three dimensional in nature. It is logical to extend the investigations to cover the behavior of heterogeneous geomaterials under 3-D conditions. A brief account of the relevant studies is summarized below.
In conventional 3-D finite element analyses, geomaterials are usually treated as homogeneous or piece-wisely homogeneous. Hexahedral or triangular elements are often used. Details of the geomaterial heterogeneity at meso-level are largely ignored.
The statistical tools or random generators were sometimes used to develop artificial 3-D distribution of different particles and components in geomaterials. For example, the three-dimensional particle flow code (PFC3D) has been developed for heterogeneous granular material analysis based on the bonded-particle method. In this method, the spherical particles of arbitrary diameters are used. The heterogeneity is generated using size distribution functions with random size particles. Applications of this method can be seen in a conference proceeding edited by Konietzky (2003). More recently, Liang et al. (2004) presented a finite element simulation of 3-D failure process in heterogeneous rocks. A statistical Weibull distribution function was used to simulate the 3-D distribution of rock particles.
For actual 3-D model development, non-destructive CT scanning techniques have been applied in rock analysis. For example, Kawakata et al. (2002) used X-ray CT techniques to investigate fault formation process of crustal rocks in triaxial compression tests. Ichikawa et al. (2001) examined the viscoelastic behavior of granitic rock by using an ultra-bright synchrotron radiation-CT method. These CT methods examined the voids or cracks inside the geomaterials. Presently, it is costly to use these techniques to establish a detailed spatial distribution of different particles and components in geomaterials because a large number of geomaterial sections have to be scanned.
In this paper, the authors propose a digital image-based approach to incorporate the actual 3-D mesostructures into numerical analysis. The method adopts several digital image techniques to extract the actual 3-D mesostructures. Firstly, digital images of geomaterial cross-section surfaces are captured. Each digital image contains the square pixels of the same size. Each 2-D square pixel can be assumed to represent a homogeneous square material on its pixel area. Secondly, this 2-D material heterogeneity from surface image can be considered as a representation of 3-D mesostructures from a very thin thickness in the depth direction. As a result, each pixel square area can be extended to form a 3-D cuboid. An iterative milling and scanning technique can be engaged to establish a series of 3-D cuboids to represent the actual 3-D geomaterial. Lastly, the digital cuboid representation of the granite sample is input into conventional numerical software package for further mechanical analysis. This paper utilizes the finite difference method for numerical analysis. The approach and numerical analysis are illustrated by a cylindrical granite sample. Details of the approach and analysis are given below.
-
Procedure for Establishing the 3-D Mesostructures
In this section, we will present the digital image based approach for establishing the digital representation of the 3-D spatial distribution of actual heterogeneous geomaterials. The process of 3-D modeling includes digital image processing, vector data transformation and iterative milling and scanning.
2.1 Digital Image and Plane Functions
A granitic rock core is used for illustration. A digital scanner is used to capture the rock cross-section surface and transform its colors into a digital image format. The rock sample is put on the scanner plate for a direct transformation. Hence, the high quality rock cross-section images can be obtained. An example of the rock core crosssection image is shown in Fig. 1.
Basically, the digital image is 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. Hence, each pixel occupies a square area. For the digital image in Fig. 1, the number of scanning lines along each of the horizontal (i) and the vertical (j) directions is equal to 60. The diameter of this rock cross-section area is 59.86mm. Therefore, the scanning line is 0.998mm wide and each pixel has a plane area of 0.996mm2.
For a gray 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, their gray levels have the integer intervals from 0 to 255. For binary images, their gray levels
Fig. 1. An illustration of rock surface image
are either 0 or 1. As a result, the digital image can be expressed as a function of piecewise continuity fðx;yÞ on plane image area in the x and y Cartesian coordinate system.
For a color image, at each pixel, three integer values are used to represent its red, green and blue colors, respectively. Hence, the color image data consists of three functions of piecewise continuity fkðx;yÞ over the plane image area, where k ¼ 1;2 or 3. The three plane functions fkðx;yÞ are the basic data in the further investigation.
It is noted that the functions fðx;yÞ or fkðx;yÞ are continuous and constants for the x and y coordinates within each square pixel and can be either continuous or discontinuous between any two connected square pixels. The continuity of the functions is stepped because it is caused by the gray level changes or the color level changes. Therefore, we can call that the functions fðx;yÞ or fkðx;yÞ are piecewise continuous over the entire image plane.
2.2 Digital Representation of Actual 2-D Granite Mesostructures
To identify the objects of interest from digital image is one of the major research areas in realm of computer science. Through different gray levels or color intensities, digital images can depict and preserve exactly the details of individual materials and their relative positions as well as spatial distribution in material surface. Digital image processing techniques such as the region segmentation method and the edge detection method are commonly used methods to identify individual materials from digital images. The region segmentation method partitions the image into regions of similar attribute. The most basic attribute for segmentation is image amplitude-brightness for gray image or color components for color image. In contrast, the edge detection method separates different materials by finding material boundaries via high gradient magnitudes among image pixels. At the different material boundaries, the first-order derivatives of the image data usually have a positive maximum or a negative minimum value. A threshold can distinguish the boundary points from the whole image pixels. Some specific methods, such as Prewitt method (Prewitt, 1970), Roberts method (Roberts, 1965), Canny method (Canny, 1986) have been developed.
In Yue et al. (2003a), digital image techniques are adopted to acquire the actual heterogeneity of asphalt concrete. The asphalt concrete is considered to have the two different individual materials. They are aggregates and matrices. Since granite is a heterogeneous solid material and consists of mainly three minerals: quartz, feldspar and biotite, Chen et al. (2004a, b) adopted a multi-threshold region segmentation method to extract the actual mesostructures from images of granite. This method tries to distinguish different minerals from the variation in their image data. The variation can be either the intensity value based on the hue, saturation and intensity (HSI) space for color images (Plataniotis and Venetsanopoulos, 2002) or gray values for gray level image. Figure 2a shows a histogram of I index of HSI color space of image in Fig. 1. Our experience in image segmentation has noted that variations of the I value are consistent with the changes of materials. The biotite materials have the lowest I values and the quartz materials have the second lowest values while the feldspar materials have the highest I values. Hence, different materials can be distinguished by using different thresholds of the I values. A careful selection of threshold on I of 0.2039 and
Fig. 2. HSI-based multi-region segmentation method for color image. (a) A histogram of I index of HSIbased color image in Fig. 1; (b) result of multi-region segmentation method
0.4352 will produce an image segmentation result in Fig. 2b. Compared with the original image, we can observe that three different minerals have been well detected.
The image mesostructures are the digital representation of the inner mineral distribution. They consist of only three different colors. Each color represents one type of three minerals. Because it occupies a square area of image, each pixel can be automatically considered as either a finite element or a finite difference grid. Moreover, each pixel has its own x and y coordinate data. A simple linear transformation can be performed to convert the image pixel into square vector data. The 2-D mesostructures can then be easily incorporated with conventional finite element or finite difference software for mechanical analysis. More details can be found in papers (Yue et al., 2003; Chen et al., 2004a, b).
2.3 Milling and Scanning System
The above method can result in a digital representation of actual granite mesostructures. The digital representation is for the rock cross-section surface. It is a digital representation of the granite heterogeneity on the straight plane. In order to have a 3- D mesostructure of the entire granite solid, an iterative milling and scanning technique is developed and used.
As discussed above, each pixel on the digital image can be considered as a square mesh or a square grid of constant material properties. Such square pixels can be used to develop cuboid elements for the 3-D mechanical analysis. The following assumptions and approach are used to extend a square pixel element into a cuboid element.
The surface image mesostructures can be considered as a representation of the inner material distribution within a small depth l. Since the rock body can be considered as continuous solid, the difference in this representation can be small as long as the depth l is very small.
As a result, the 2-D square element can be extrapolated to form a 3-D cuboid element with a uniform thickness l in the depth direction (z axis). Each cuboid element will have the same material properties as its parent square pixel element.
The surface rock layer of thickness l can then be milled with a rock core grinder machine. A new rock surface will be exposed and digitized with a digital scanner.
Similarly, the new surface image mesostructures after milling can be used to develop the second layer of 3-D cuboid elements.
The above process can be repeated for the third layer, the fourth layer...of the 3- D cuboid elements.
As a result, the entire granite specimen can be represented by a series of layers of the cuboid elements, where each cuboid element is homogeneous and has the properties of an assigned mineral.
Figure 3 shows an example of the transformation process from a 2-D plane image mesh to 3-D layer of cuboid mesostructures. The thickness l should be made equal to the scanning line thickness h. As a result, the cuboid element can be a cube, which will make the cuboid digital approximation have the similar accuracy level along the three orthogonal directions (x, y, z) as shown in Fig. 3 c.
An automatic rock core grinder machine is used to remove the thin rock layer and to form fresh and flat rock cross-section surface. The rock sample is fixed firmly on a
Fig. 3. An example of transformation from two-dimensional surface image to three-dimensional cuboid element. (a) Two-dimensional surface image; (b) two-dimensional square element; (c) three-dimensional cuboid element
steel core jig. This core jig ensures that the rock core surface is placed parallel to the grinder machine surface. Then a thin layer of rock in uniform thickness l is removed by milling. A new rock cross-sectional surface is formed and parallel to the grinder surface. The procedure to establish a digital representation of the actual 3-D rock mesostructures can be summarized as follows:
Step 1: A rock specimen is put on a scanner plate to transform the rock surface into a digital image.
Step 2: Digital image processing techniques are used to extract the 2-D rock actual mesostructures from the plane image.
Step 3: The actual plane image mesostructures are transformed as square vector mesostructures.
Step 4: The 2-D square mesostructures are extrapolated to form a one layer of cuboid mesostructures by elongating the zero thickness 2-D square mesostructure to a thickness l along the depth (z) direction.
Step 5: The rock specimen surface is removed again by the determined uniform thickness l via a rock core grinder machine. The removed thickness is measured precisely by the digital gauge and scale.
Step 6: Steps (1) to (5) are repeated until the internal mesostructures of the entire rock core has been measured and digitized.
2.4 Digital Representation of 3-D Rock Mesostructures
In this section, a cylindrical rock sample is used as an example to develop a digital representation for its 3-D mesostructures. In order to build accurate 3-D mesostructures, the standard practice for preparing rock core specimens and determining dimensional and shape tolerances stated in ASTM-D4543 is adopted for the tolerance checks on the straightness of the elements on the cylindrical surface, the flatness of the end bearing surface, and the perpendicularity of the end surfaces with the axis of the core.
Moreover, Fig. 4 shows the digital image acquisition system as well as its 3- D coordinate axis for building correct 3-D mesostructures. Two fixed points are marked on the scanner plate. At the corresponding positions, two vertical lines along the depth direction are drawn on the rock cylindrical surface. Every time when the designed thickness l is removed and the rock sample is placed on the scanner plate, the two lines have to match with two points separately for exact location calibration. The approach ensures that the digital images of the rock cross-section surfaces at different heights can match exactly in horizontal directions. By this way, several 2-D mesostructures from different layers can form precise 3-D mesostructures.
An ideal image resolution is important for the 3-D modeling. A high image resolution will have very small pixels and consequently result in too many cuboid elements.
Fig. 4. A schematic illustration for building correct 3-D actual mesostructures. (a) Front view of scanned rock sample; (b) left view of scanned rock sample; (c) plan view of scanned rock sample; (d) rock surface
image acquisition system
Fig. 5. The generated actual rock 3-D heterogeneity
On the other hand, a low resolution may have a higher degree of ignorance in the fine details of rock material heterogeneity and may build a 3-D model without enough information of material spatial distribution.
In this study, the number of pixels along each of the x and y coordinates is determined to be 60. The diameter of rock cylinder is 59.86mm. The length and width of each cuboid element are both 0.998mm. The mechanical apparatus of the grinder machine can produce a precise removed depth l of 0.991mm. This depth l is close to the width of cuboid.
For illustration purpose, a total of 12 layers of the 3-D mesostructures have been built. A digital representation of 60 60 12 cuboids is established to represent the rock internal structure. Thus, this model has 43200 cuboid elements. The entire modeled rock sample is a disc of a diameter 59.86 mm and a height 11.89 mm. By inheriting the same 3-D coordinate system in Fig. 4, the digital representation of the entire 3-D mesostructures is shown in Fig. 5, where 12 digital images for 12 rock cross-section surfaces are used. In Fig. 5, three colors are used to represents the three main minerals Feldspar, Quartz and Biotite. Each cuboid element represents one of the three minerals and has its assigned constant material properties.
3. Numerical Evaluations
The above cuboid element mesh can be incorporated into a conventional finite element or finite difference software package. In this simulation, 3-D finite difference code-FLAC3D (Itasca, 1997) is used. FLAC3D contains a programming language FISH for user to define new variables and functions. Following the rules of FISH language, the FISH program is a text code which can automatically generate the geometry and material distribution as well as loading conditions of a very complicated 3-D model.
According to the 3-D mesostructures obtained from digital images and milling and scanning system, there are 12 layers along the z direction. At each layer, the number of elements along x and y direction is 60. The material distribution of each layer is determined by the image mesostructures. The type of color of each cuboid represents its material type. Hence, the coordinate of each cuboid in the 3-D grid with the following expression of model and its material type can be determined precisely. A FISH program can be developed to automatically incorporate the actual 3-D mesostructures into the numerical software.
Firstly, a 3-D grid system of cuboid elements 606012 is generated. Secondly, a material type number is assigned to this element according to the color number of each cuboid element in each layer. Thirdly, the 3-D grid is adjusted to the actual size. Finally, the material type (i.e. either linear-elastic or Mohr-Coulomb material) and material properties (i.e., elastic modulus, Poisson's ratio, tensile strength, friction angle and cohesion) are assigned to each cuboid element.
The conventional Brazilian tensile test is analyzed numerically in this paper. This test is commonly used to estimate the rock tensile strength indirectly. A disc or cylindrical sample is loaded diametrically between the platens of a testing machine. The vertical compressive force is applied to induce a tensile failure where the disc or cylindrical sample is split along the vertical loading diameter. The boundary conditions of the 3-D Brazilian test model are shown in Fig. 6.
Several numerical analyses are performed below to examine the effects of material heterogeneity. The results of homogeneous and heterogeneous simulations are compared. Also, the elastic-plastic studies are carried out to investigate the failure process and to estimate rock tensile strength.
Fig. 6. Boundary conditions of the 3-D Brazilian test model. (a) Plan view of 3-D Brazilian test model;
(b) side view of 3-D Brazilian test model
3.1 Linear-Elastic Analysis
In the homogeneous cases, all cuboids in the 3-D model are assumed to have the same material properties. Elastic modulus of 70GPa and Poisson's ratio 0.3 are assigned to the three rock minerals.
Figures 7 and 8 show the plane view of contour lines of the numerical results for the maximum and minimum principal stresses in the homogeneous elastic rock disc, respectively. The contour lines are symmetric with respect to the loading plane.
In the heterogeneous case, the cuboids are assumed to have different material properties according to the minerals. Quartz, feldspar and biotite are assumed to have elastic modulus values of 90, 70 and 40GPa, respectively. Their Poisson's ratio values are equal to 0.3. These values for the mineral elastic constants are determined based on the references (Li, 2001; Vutukuri et al., 1974). Other relevant data are presented in Table 1.
Fig. 7. Plane view of contour lines of minimum principal stress in the 3-D Brazilian homogeneous model (interval 30 kPa )
Fig. 8. Plane view of contour lines of maximum principal stress in the 3-D Brazilian homogeneous model
(interval 0.1 MPa )
Fig. 9. Plane view of contour lines of minimum principal stress in the 3-D Brazilian heterogeneous model (interval 30 kPa )
Fig. 10. Plane view of contour lines of maximum principal stress in the 3-D Brazilian heterogeneous model (interval 0.1 MPa )
Figures 9 and 10 show the plane view of contour lines of the numerical values for the maximum and minimum principal stresses in the heterogeneous elastic rock disc, respectively. Compared to the contours in Figs. 7 and 8, it is evident that the material heterogeneity alters the distribution of stresses and changes the contour lines from symmetric to irregular pattern.
Figures 11 and 12 present the numerical results for the non-dimensional normal and radial stress along the loading plane of the heterogeneous elastic rock disc, respectively. The normal stresses in one disc layer (i.e., the 6th layer) are presented for simplicity. The elastic modulus values of the cuboids along the loading plane for the disc layer 6 are also plotted at the lower portion of each of Figs. 11 and 12. Along the line in the graph, the highest points represent quartz materials, the middle points represent feldspar materials, and the lowest points represent biotite materials.
Fig. 11. Comparison of the non-dimensional hoop normal stress along loading diameter in the 6th layer
Fig. 12. Comparison of the non-dimensional radial normal stress along loading diameter in the 6th layer
The stress along the loading diameter of the 6th layer varies from place to place due to the different mineral distribution. Stresses always change greatly at the boundaries of different minerals. The maximum tensile stress points sometimes are not at the center of the loading axis due to the influence of material heterogeneity. It means that as the loading increases, the tensile failure might not occur at the center of the disc. This demonstrates the characteristic of the failure model of the heterogeneous case.
3.2 Linear-Plastic Analyses
The classical Mohr-Coulomb plastic model with tension cutoff is adopted in FLAC3 D to simulate the failure behavior of the three main granite minerals. For the MohrCoulomb model, a shear yield function and a non-associated shear flow rule are used. The corresponding mechanical parameters for this model are elastic modulus, Poisson's ratio, cohesion, internal friction angle and tensile strength. Similarly, the same homogeneous and heterogeneous approaches are adopted for the failure predictions and comparisons.
In FLAC3D, when the status of a cuboid element violates the failure criteria, the element will fall into plastic stage and will be marked as failed elements. In this numerical simulation, the failure can be represented by a collection of failed square elements and these failed elements may connect to form elongated regions that correspond to macroscopic fractures. A FISH program is developed to record the plastic elements in FLAC3D model while the loading is increased. According to the numbering sequence of cuboid elements, the exact position of the failed elements can be determined.
3.2.1 Homogeneous Cases
All minerals are assumed to have the same material properties in FLAC3D model. Based on experience and data of other references (Li, 2001; Vutukuri et al., 1974), the mechanical properties for the homogeneous rock body are given in Table 1. Three cases are studied:
In case H1, all the three minerals are assumed to have the same property values as those of quartz.
In case H2, all the three minerals are assumed to have the weighted average values of those of quartz, feldspar and biotite.
In case H3, all three minerals are assumed to have the same property values as those of biotite.
In addition, the assumed property values of the three minerals are listed under the case HE for the actual heterogeneous case.
In order to determine the collapse of the disc, the histories of the x displacement of the point A and point B are monitored. The results are shown in Fig. 13. The points A and B are on the boundary of the disc at the coordinates: A (59.86 , 29.93, 11.88 mm) and B (0.0, 29.93, 11.88 mm). The two points can form a radius which is perpendicular to the loading diameter. It can clearly observe that there is a turning position in the curves. After that point, the x displacements at the two monitored points increase rapidly. Also, the length of the induced crack exceeds half of the diameter. It indicates that the disc begins to rupture and the system is collapsing. The measured force at this stage is defined as the failure force.
Figure 14a–c show the process of crack initiation and propagation, where the dark grids indicate elements that have violated the tensile failure criterion. The three homogeneous cases have the same failure patterns, but the stress magnitudes are
Fig. 13. Histories of x displacement of monitored points
different. Crack propagation process in three different layers (the 1st, 6th, and 12th layers) is shown for a clear presentation in Fig. 14.
As shown in Fig. 14a, the crack first initiates in the 1st and the 12th layer, but no crack occurs in the 6th layer. The exact coordinates of the failed cuboids on the 1st and 12th layers are (30, 30, 1), (30, 31, 1), (30, 31, 1), (31, 31, 1), (30, 30, 12), (30, 31, 12) , (30, 31, 12), (31, 31, 12). The number of cuboids along each of the x and y coordinates is 60. It means that the crack initiates at the center of the loading diameter of the two surface layers. Then the crack propagates along the loading diameter and towards the two ending positions simultaneously on the x–y plane. In the mean time, the crack zones also occur along the z direction inside the rock layers. Failed cuboids appear at the center of the 2nd and 11th layer, and finally appear at the center of the 6th layer. At each layer, the failed cuboids will propagate vertically towards the loading positions. They also follow the same crack propagation pattern on the x–y plane as shown in the 1st and 12th layer in Fig. 14. This crack propagation process is symmetric to both the central horizontal x axis and the central z axis.
Finally, these failed cuboids will form a 12-layer-thick and 2-grid-wide line at the center of the rock cylinder. All the failed square cuboids violate the tensile failure criterion only. Thus the vertical seam can be considered as a tensile crack line along the central vertical position. The tensile strength of the rock is then calculated via the measured applied boundary force as shown in Table 2. It is noticed that the predicted tensile strengths for the three homogeneous cases agree well with the pre-given assumptions.
3.2.2 Heterogeneous Cases
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 heterogeneous case no. HE.
Fig. 14. Failure propagation process in homogeneous rock cylinder (dark grids represent the tension failed cuboids). (a) Stage 1 (measured force is equal to 92% failure force); (b) Stage 2 (measured force is equal to 95% failure force); (c) Stage 3 (measured force is equal to 100% failure force)
Table 2. Results of the total applied failure forces from the present numerical method for the homogeneous cases and the conventional analytical method for homogeneous rocks
Case Numerical prediction Analytical prediction Difference between A and B no. of the failure forces
(A) (kN) Given tensile Predicted failure Difference Percentage
strength (MPa) force (B) (kN) (AB) (kN) (AB)=B100 %
| H1 | 16.26 | 14.0 | 16.24 | 0.02 | 0.1 |
| H2 | 13.59 | 11.7 | 13.53 | 0.06 | 0.5 |
| H3 | 8.13 | 7.0 | 8.12 | 0.01 | 0.1 |
Fig. 15. Failure propagation process in heterogeneous rock cylinder (white grids represent the tension failed cuboids and other colors represent material heterogeneity). (a) Stage 1 (measured force is equal to
80% failure force); (b) Stage 2 (measured force is equal to 88% failure force); (c) Stage 3 (measured force is equal to 93% failure force); (d) Stage 4 (measured force is equal to 100% failure force)
Figure 15a–d show the simulated crack initiation process for the heterogeneous rock specimen under the Brazilian loading condition, where the white grids indicate tensile failed cuboids and other colors represent different minerals. Crack propagation process in the 1st, 6th, and 12th layer is also shown. The failure load is determined using the approach discussed in the above three homogeneous cases.
First, several failed cuboids appear in the disc. Unlike the homogeneous cases, these failed elements are offset from the center of the disc. As the loading increases, more failed cuboids appear at other locations. These failed cuboids are very close to the loading axis. All failed cuboids violate only the tensile failure criteria. The accumulation of failed elements eventually forms a tensile crack line insider the disc body. The cracks begin to appear in the disc when the measured force is about 80% of the failure load. In contrast, in the homogenous cases, the applied force inducing the initial tensile cracks is about 92% of the corresponding failure load. It is believed that material heterogeneity generates some high stress concentrations and causes the tensile failure to occur earlier.
In Fig. 15, the initial cracks appear in the 1st layer at the cuboids Nos. (29, 43, 1) , (31, 43, 1) and (31, 40, 1). These cuboids are either quartz or feldspar that surround a large biotite. This biotite region is located near the y axis and is the biggest biotite region adjacent to the center. In the 6th and 12th layers, cracks also usually appear at the feldspar or quartz cubes. These failed cubes are also close to the biotite cubes. This numerical result shows that the tensile stress and the tensile strengths control the failure model in the Brazilian test for the stiff granite. The actual spatial distribution of the three main minerals dominates actual crack patterns in the granite.
For the purpose of comparison between the crack patterns under 2-D and 3- D conditions, the 1st layer mesostructures are further used to examine the failure pattern under a 2-D plane stress condition.
The 2-D plane stress problem is calculated again with the FLAC2D. The 1st layer heterogeneity is incorporated into this modeling. The material properties of the three minerals are the same as those in modeling with FLAC3D. The 2-D plane stress results are shown in Fig. 16. As shown in Fig. 16a, the cracks initiate at the elements (29, 41) and (29, 40). The two elements are slightly different from the failed cubes in the 3- D
Fig. 16. Crack pattern in 2-D model with the 1st layer material heterogeneity (white line represent the crack and other colors represent material heterogeneity). (a) Crack initiation points; (b) crack propagation path
simulations shown in Fig. 15a. The crack then propagates downward and upward toward the two loading ends as shown in Fig. 16 b.
The results shown in Figs. 15 and 16 clearly indicate that the 2-D crack pattern is different from the 3-D crack pattern. Furthermore, the 3-D heterogeneity has a higher influence on the stress distribution and crack propagation process than the 2-D heterogeneity.
3.2.3 Analysis and Discussions
Figure 17 shows a comparison among the relations between the applied displacement and the associated total force F at the loading position from the initial loading to the post-peak behavior for the three homogeneous cases and the heterogeneous case. Table 2 gives a summary of the results of the applied forces and the analytical prediction of the peak tensile force for the three homogeneous cases. Table 3 gives a comparison between the given tensile strength values and the predicted tensile strength values.
Fig. 17. Force and displacement responses for different cases
Table 3. Results of the tensile strength prediction from the present numerical method
| (MPa) | A (MPa) | Difference (AB) (MPa) | Percentage (AB)=B100 % | |
| H1 | 14.0 | 14.02 | 0.02 | 0.1 |
| H2 | 11.7 | 11.72 | 0.02 | 0.1 |
| H3 | 7.0 | 7.01 | 0.01 | 0.1 |
| HE | (Quartz) 14.0 | 8.49 | 5.51 | 39.4 |
| (Feldspar) 11.0 | 2.51 | 22.8 | ||
| (Biotite) 7.0 | 1.49 | 21.3 |
Case no. Given tensile Numerical predicted results Difference between A and B strength (B) of tensile strengths at failure
From the results presented in Fig. 17, Tables 2 and 3, the following important observations can be made:
-
For three homogeneous cases H1, H2 and H3, the numerically predicted peak loads at the split tensile failure are almost the same as analytical predicted results. The numerical predicted tensile strength values are also consistent with the given values. Furthermore, the corresponding displacements at the loading point are higher for the cases with smaller modulus values. These results are within the expectation, which also is an indirect verification of the numerical computation.
-
For the heterogeneous case HE, its load and displacement curve is similar to that of the case H2 before the tensile failure at the peak (Fig. 17). This result is also expected since case H2 has the weighted average value of the modulus values associated with the case HE.
-
However, the numerically predicted failure load for the case HE is much less than that of the homogeneous cases H2. Although case HE and case H2 have the same overall tensile strengths, the calculated tensile strength for case HE decreases significantly. The tensile strength of the heterogeneous case is about 27.4% less than the tensile strength of case H2, which indicates that the material heterogeneity greatly influences the mechanical strength of granitic rocks.
3.3 Applications of the Proposed Method
The present literature review has indicated that a majority of current investigations in geomechanics could be classified into two categories, namely (1) studies based on the phenomenological approach and (2) studies based on the statistical approach.
The phenomenological approach, which considers a geomaterial to be an ideal homogeneous or piecewise homogeneous material, is a mainstream method in civil and geotechnical engineering. It has been widely adopted in theoretical studies, laboratory testing and numerical modeling. Although phenomenological-based studies have made a great contribution towards the development of present geomechanics, this simplifying assumption induced some problems. For example, in current laboratory system, the geomaterial specimens are assumed to be homogeneous although most specimens are heterogeneous materials which consist of different individual materials and components. The test data, such as the stress-strain curves, gained from physical laboratory tests, actually are overall responses of the specimens under external loading. An overall response for each specimen is an averaging result of the responses of the specimen at numerous different positions. Hence, the material heterogeneity is largely ignored and the phenomenological approach-based computational results are not precise and sometimes cannot predict the mechanical behaviour of practical engineering work.
Studies based on the statistical approach have generated new insights into the mechanical behavior of heterogeneous geomaterials. It has been shown that the analysis based on statistical data can take the material heterogeneity into consideration. However, although achievements were made in these studies, it seems that there are few publications available in the relevant literature to discuss the extent and accuracy of how the statistical approach-based studies represent the actual mechanical behavior of geomaterials.
Based on the above discussion, the possible applications of the proposed digital image based numerical approach for geomechanics are:
-
This approach takes the actual material heterogeneity into the mechanical consideration and can be a better method to study the complex material properties of geomaterials.
-
By introducing the actual material heterogeneity and mesostructures into numerical models, digital image based numerical approach can employ well-established constitutive models to examine the mechanical performance of materials. This can establish a direct connection between the mechanical properties and information of material spatial distribution.
-
It is possible to accurately predict the mechanical behavior of geomaterial under external loadings. For example, studies based on actual heterogeneity could explain the variation of material properties of specimens in laboratory tests.
-
Comparing the data of actual mesostructures-based numerical tests with data from physical tests, back analysis can be conducted. It is possible to provide a coupling method to determine the constitutive models and material parameters of individual materials.
-
It is possible to develop a coupling method to incorporate the structures of geomaterials at different levels (micro-, meso- and macro-) and examine the mechanical properties of geomaterials at all scales.
4. Conclusions
This paper has presented a digital image based approach for examining rock properties by taking the 3-D heterogeneity into consideration. An iterative milling and scanning system and several digital image techniques are adopted to generate the 3-D actual rock mesostructures. Using a granite cylinder as an example, the paper has presented the procedure to establish the actual 3-D mesostructures.
The finite difference code FLAC3D is adopted to simulate Brazilian tensile tests. The stress analysis, crack pattern and failure prediction are performed. The numerical results presented have demonstrated that the material heterogeneity has significant effects on the stress distribution, crack initiation and maximum failure force under loading. The numerical results have also shown that the maximum tensile stress does not occur at the centre of a circular disk if rock heterogeneity is considered.
The crack initiation and propagation can be recorded at each step of the numerical calculations. It is shown that the three different mineral distributions considered in the paper control the crack initiation and failure pattern for granite subject to the Brazilian test. The results have further shown that a majority of the numerically predicted failures compare well with the analytical results and are within expectations.
The method presented in this paper is an efficient approach to extend the conventional and powerful numerical methods to mechanical analysis of actual heterogeneous materials subject to loading. Although this paper only uses a rock cylinder as an example, it is believed that this method can be incorporated with conventional macromechanical approaches for all heterogeneous geomaterials such as asphalt concrete, cement concrete and soils.
Acknowledgements
The authors would like to thank the Research Grants Council of Hong Kong SAR Government for the financial support. One of the authors Mr. S. Chen would also like to thank the University of Hong Kong for granting him a postgraduate scholarship.
References
Canny, J. (1986): A computational approach to edge detection. IEEE Trans. Pattern analysis and machine intelligence 8, 679–714.
Chen, S., Yue, Z. Q., Tham, L. G., Lee, P. K. K. (2004a): Modeling of the indirect tensile test for inhomogeneous granite using a digital image-based numerical method. Int. J. Rock Mech. Min. Sci. 41, 447–447 (SINOROCK Paper No. 2B 01, CD ROM).
Chen, S., Yue, Z. Q., Tham, L. G. (2004b): Digital image-based numerical modeling method for prediction of inhomogeneous rock failure. Int. J. Rock Mech. Min. Sci. 41, 939–957.
Cundall, P. A., Strack, O. D. L. (1979): A discrete numerical model for granular assemblies. Geotechnique 29, 47–65.
Fang, Z., Harrison, J. P. (2002a): Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks. Int. J. Rock Mech. Min. Sci. 39, 443–457.
Fang, Z., Harrison, J. P. (2002b): 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. 39, 459–476.
Frost, J. D., McNeil, S. (1998): Proc., Second International Conference on Imaging Technologies: Techniques and Applications in Civil Engineering. A.S.C.E. Switzerland, 354.
Frost, J. D., Wright, J. R. (1993): Proc., EF=NSF Conference on Digital Image Processing: Techniques and Applications in Civil Engineering. A.S.C.E. Hawaii, 301.
Ichikawa, Y., Kawamura, K., Uesugi, K., Seo, Y. S., Fujii, N. (2001): Micro- and macrobehavoir of granitic rock: observations and viscoelastic homogenization analysis. Comput. Methods Appl. Mech. Engng. 191, 47–72.
Itasca (1997): FLAC3D Fast Lagrangian analysis of continua in 3 Dimensions, version 2.0. Itasca Consulting Group Inc, Minnesota.
Kawakata, H., Cho, A., Yanagidani, T., Shimada, M. (2002): The observations of faulting in westerly granite under triaxial compression by x-ray CT Scan. Int. J. Rock Mech. Min. Sci. 34 , 375.
Konietzky, H. (2003): Numerical modeling in micromechanics via particle methods. In: Proc., First International PFC Symposium, Gelsenkirchen, Germany, November 2002 , 321 pp.
Li, L. (2001): Microscopic study and numerical simulation of the failure process of granite. Ph.D. thesis, Department of Civil Engineering, University of Hong Kong, Hong Kong.
Liang, Z. Z., Tang, C. A., Li, H. X., Zhang, Y. B. (2004): Numerical simulation of the 3- D failure process in heterogeneous rocks. In: Hudson, J. A., Feng, X. T. (eds.). Proc., ISRM SINOROCK 2004 China.
Plataniotis, K. N., Venetsanopoulos, A. N. (2002): Color image processing and applications. Springer, Berlin Heidelberg New York, pp 355.
Prewitt, J. M. S. (1970): Object enhancement and extraction. In: Lipkin, B. S., Rosenfeld, A. (eds.) Picture processing and psychopictorics. Academic Press, New York.
S. Chen et al.: Digital Image Based Approach for 3-D Analysis
Roberts, L. G. (1965): Machine perception of three-dimensional solids. In: Tippett, J. T. et al. (eds.), Optical and electrooptical information processing. MIT Press, Cambridge, MA, 159–197.
Tang, C. A., Liu, H., Lee, P. K. K., Tsui, Y., Tham, L. G. (2000): Numerical studies of the influence of microstructure on rock failure in uniaxial compression–part I: effect of heterogeneity. Int. J. Rock Mech. Min. Sci. 37, 555–569.
Vutukuri, V. S., Lama, R. D., Saluja, S. S. (1974): Handbook on mechanical properties of rocks, vol. II. Trans Tech, Clausthal.
Yue, Z. Q., Chen, S., Tham, L. G. (2002): Digital image processing based finite element method for rock mechanics. In: Li, Y., Tang, C., Feng, X., Wang, S. (eds.) Proc., 2nd International Conference on New Developments in Rock Mechanics and Rock Engineering, Supplement. October 10–12, 2002, Shenyang, China. Rinton Press, Princeton, USA, 609–615.
Yue, Z. Q., Chen, S., Tham, L. G. (2003a): Finite element modeling of geomaterials using digital image processing. Comput. Geotech. 30, 375–397.
Yue, Z. Q., Chen, S., Tham, L. G. (2003b): Seepage analysis in inhomogeneous geomaterials using digital image processing based finite element method. In: Proc., 12th Panamerican Conference for Soil Mechanics and Geotechnical Engineering and the 39th US Rock Mechanics Symposium, Soil and Rock America 2003. MIT, Boston, USA, June, 2003 , 1297–1302.
Authors' address: Dr. Z. Q. Yue, Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, P.R. China; e-mail: yueqzq@hkucc.hku.hk
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
