Analyzing fracture properties of the 3D reconstructed model of porous rocks

Engineering Fracture Mechanics xxx (2017) xxx–xxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Analyzing fracture properties of the 3D reconstructed model of porous rocks Xiao-Ping Zhou ⇑, Nan Xiao School of Civil Engineering, Chongqing University, Chongqing 400045, PR China

a r t i c l e

i n f o

Article history: Received 2 May 2017 Received in revised form 22 October 2017 Accepted 23 October 2017 Available online xxxx Keywords: The 3D reconstructed model The topological and mechanical properties The two-point correlation function Lineal-path function Fractal function

a b s t r a c t The non-negligible fissures, cracks or pore structures, which are contained in rocks, have significant influences on their mechanical properties. In this paper, a novel method is developed to reconstruct a three-dimensional multiple phase rock model from the limited morphological information of Chongqing sandstone. Firstly, X-ray micro-CT experimental tests are applied to Chongqing sandstone to obtain the raw CT images and the mechanical parameters. Then, the two-point correlation function, lineal-path function and fractal function are combined to generate the reconstructed model. Next, the reconstructed model subjected to compressive loads is simulated by Finite Element-Smoothed Particle Hydrodynamics (FE-SPH) coupling method with a brittle model formulation. Finally, the initiation, propagation and coalescence of cracks are investigated during the numerical simulation. The topological and mechanical properties of the reconstructed model are carefully analyzed and compared with that of the reference model. The comparison of the numerical modeling results and the experimental data shows a good consistency. This indicates that the proposed method has the ability to predict macroscopic behaviors of rocks from limited experimental information. Ó 2017 Published by Elsevier Ltd.

1. Introduction With the development of the science and technology, the world is recognized in macroscopic and microcosmic scale. The advanced numerical reconstruction methodology is the essential mean to access the various types of information about the complex microstructure in the multiple phase materials such as rocks and is also capable of predicting the behavior of the multiple phase materials to help us better recognize this world. The formed irregular micropore structures in rocks significantly affect the macro-mechanical properties of porous rocks, such as porosity, permeability, elastic modulus, Poisson’s ratio, compressive and tensile strength. Among them, the topological properties, such as porosity and permeability, are directly influenced by the micropore structures, and the mechanical properties, such as elastic modulus, Poisson’s ratio, compressive and tensile strength, are indirectly related to the microstructures. Thus, the better understanding of the distribution and characteristics of microstructure in rocks is crucial to the successful performance of rock engineering. Rocks are composed of solid and pore phases, namely, two-phase materials. The pore phase usually coexists with the solid phase in the form of joints, fissures and cracks. Some experimental methods can be used to obtain a sequence of twodimensional images from the multiple phase materials such as serial sectioning [1–3], confocal laser scanning microscopy [4], FIB-SEM [5], micro-CT [6,7] and Nano-transmission X-ray microscopy [8]. Compared with other methods, the ⇑ Corresponding author. E-mail address: [email protected] (X.-P. Zhou). https://doi.org/10.1016/j.engfracmech.2017.10.021 0013-7944/Ó 2017 Published by Elsevier Ltd.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

2

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Nomenclature C1; C2 e0

centers of the maximal balls difference between the reference control functions and the calculated control functions for the two-point probability function e1 difference between the reference control functions and the calculated control functions for the lineal-path function in the pore phase difference between the reference control functions and the calculated control functions for the lineal-path funce2 tion in the solid phase e3 difference between the reference control functions and the calculated control functions for the fractal function in the pore phase difference between the reference control functions and the calculated control functions for the fractal function in e4 the solid phase b target mean energy Er g positive constant T control parameter of the system temperature T0 the initial temperature u distance between two randomly selected point v a constant parameter controlling the annealing speed wr standard deviation of the energy distribution x1 ; x2 two randomly selected points in the system C general control function DE ¼ Ej  Ej1 constant porosity k the temperature control parameter rI the direct tensile components of stress in a tension softening model rII the shear component of stress in a shear softening/retention model n accepting probability

micro-CT method is a nondestructive method to obtain the information of the composition, structure, and properties of the nontransparent rock core samples. Then, based on the generated image data, the geometrical model which contains the configuration of the pore structures in rocks can be accurately reconstructed by a series of mathematical procedures, such as Marching cubes algorithm, Delaunay refinement algorithm[9–16]. Although the reconstructed model is able to partially characterize porous rocks in the numerical simulation for prediction purpose, the size of the reconstructed model is too small and limited to the resolution of the CT scanner. In order to analyze porous rocks under different engineering condition and in different engineering zones, countless CT observations may be needed. Due to the high cost of various rock samples preparation, several stochastic methods have been developed to reconstruct multi-phases and heterogeneous model for predicting the behaviors of rocks and giving helpful advice to engineers, such as geological process-based reconstruction method [17,18], Gaussian random fields reconstruction method [19–21], multiple-point statistics reconstruction method[22–24] and simulated annealing reconstruction method[25–28]. Among them, the simulated annealing reconstruction method appears to be the most attracting method for its capability to introduce any number and arbitrary type of control function. With this capability, the reconstructed rock model has the possibility to characterize the natural rocks which contain various and complex pore structures. Moreover, with the appropriate introduction of the correct control function, the reconstructed rock model is able to predict the natural rock performance in various engineering project under different geological conditions. Especially, cracks are often initiated from the micro defects in the porous materials such as porous rocks. The micro defects, such as micropore structures, are the key factor to the fracture mechanism of porous rocks. With the appropriate introduction of the comprehensive consideration of the environment and the geology, the reconstructed rock model is able to predict the engineering disasters caused by the fracture of rocks, such as the rockfall and landslide caused by the crack propagation and coalescence, the mud-rock-flow induced by the coupling of the seepage and the expansion of the pore structure, and the rock burst in the deep stratum caused by the fast crack propagation with high elastic potential energy. Torquato [26] formulated a procedure to reconstruct the disordered heterogeneous materials based on simulated annealing algorithm. The next year, Torquato [27,28] extended the previous methodology developed for dispersion and formulated a methodology with the lineal path function and the two-point correlation functions to reconstruct 1D, 2D and 3D general random heterogeneous media from limited morphological information. Torquato [27,28] also employed lower order statistical function to reconstruct microstructures with stochastic optimization. However, he found that the lower order statistical function is not capable of fully characterizing the two-phase heterogeneous materials, and more than one solution is suitable for a specific lower order correlation function. Furthermore, the higher-order correlation functions can be used to illustrate the contribution of shape and geometry effects [29]. Later, Sheehan and Torquato [30] optimized the specified statistical corPlease cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

3

relation functions to eliminate the artificial anisotropic effects on a system possessing significant short-range order. Kröner [31,32] and Beran [33] also proposed a theoretical model with the N-point statistical functions to illustrate the effects of the volume fraction and the number of the phase in the multiple-phase model. Manwart and Hilfer [34] proposed a novel method to reconstruct random media based on Monte Carlo methods. Jiao and Stillinger [35,36] introduced a novel heterogeneous-media reconstruction algorithm called ‘‘lattice-point”, and proposed a theoretical method to reconstruct the heterogeneous materials based on two-point correlation functions. Davis et al. [37] characterized the microstructure of the two-phase media and generated statistically similar numerical reconstructed model based on simulated annealing algorithm. Ju et al. [38] presented a novel method for reconstructing the microstructure of sandstone. In this paper, an exhaustive numerical model is established by using the proposed reconstruction methodology. Then, the reconstructed model is simulated by using Finite element-Smoothed Particle Hydrodynamics (FE-SPH) coupling method with the Maximum Stress Criteria (Rankine criterion), which has the ability to predict the failure of brittle rocks. To verify the proposed reconstruction method, Chongqing sandstone samples are used. The SOMATOM scope is employed to generate X-ray computed tomograph (CT) images from Chongqing sandstone samples. A series of triaxial compression experiments are conducted on Chongqing sandstone samples by the multi-function triaxial rock testing system. The comparison results show that the proposed reconstruction method is able to safely predict the macroscopic behaviors of porous rocks by limited experimental morphing information. 2. X-ray micro-CT observations and experimental testing procedure 2.1. X-ray micro-CT observation X-ray CT has been widely applied to investigate the influence of internal microstructures on the mechanical behaviors of porous rocks. The X-ray CT scanner employed in this paper is named SOMATOM Scope made by SIEMENS, as shown in Fig. 1 (a). The basic working principles of SOMATOM Scope can be described in Fig. 1(b). During the working process of the CT device, the rock sample is fixed on the middle of the cardboard boxes on the scanning platform to sufficiently sustain the X-ray beam. The location of the rock sample is firstly calibrated by the infrared calibration instrument before the X-ray Micro-CT observations. Then the X-ray source continues emitting X-ray beam which can penetrate the sample placed on the sample stages. The emitted X-ray is attenuated to a certain degree according to the composite of the sample, and the attenuated X-ray is received by the receiver as raw information for imaging and reconstruction. The values of the pixel in the image data are directly proportional to the attenuation of X-ray which follows the famous Maxwell’s equations. During the CT scanning process, the detector and receiver are rotated by the scanning ring of the CT machine to obtain raw data or image data from various directions. The detector and receiver achieve extremely fast decay and afterglow times due to the adoption of Ultra Fast Ceramics (UFC) materials which have the high X-ray absorption combined with extremely effective conversion of X-ray energy into visible light reducing noise in the images. The whole process is monitored by Syngo Scope one console, a workplace solution combining powerful hardware with a complete clinical application suite. The X-ray tube of this system has a focal spot size 0.8  0.5 mm. High-contrast resolution is 15.6 lp/cm under 2%MTF. The X-ray beam penetrating the rock sample is measured by 24 detector rows. Using Siemens’s proprietary Ultra-Fast Ceramics detector, SOMATOM Scope provides excellent detailed images. The electron current is 81 mA, the accelerating voltage is 130 kV and the scanning time is 1500 ms. The reconstruction matrix of the Syngo Scope one console is 512  512. Three rock samples employed in this paper are respectively marked by C2, C3, C4. The rock samples are scanned by SOMATOM scope to obtain the CT images, as shown in Fig. 2(a) and (b). Before the novel reconstruction, the raw CT images should be pre-processed. In this paper, a novel denoising algorithm [39] is employed to denoise the CT images, as shown in Fig. 3(a). Then, the Isodata

Fig. 1. CT scanner and working principle of the X-ray CT observation: (a) SOMATOM scope CT scanner. (b) Working principle of the X-ray CT observation.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

4

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 2. Chongqing Sandstone sample: (a) the raw CT images of samples C2 and C3. (b) The raw CT images of sample C4.

(a)

(b)

(c)

Fig. 3. The CT images pre-processing: (a) the raw CT images, (b) the CT images after denoising algorithm, (c) the CT images after threshold algorithm.

threshold algorithm [40] is employed to segment the CT images, as shown in Fig. 3(b). Finally, the images suitable for numerical reconstruction are obtained, as shown in Fig. 3(c).

2.2. Experimental testing procedure Chongqing sandstone is a typical porous rock. We choose Chongqing sandstone for numerical reconstruction and simulation. First, the mechanical properties of Chongqing sandstone are defined by experiments. Then, the mechanical properties of the reconstruction or reference model are determined by the numerical procedure. Finally, the mechanical properties obtained from the numerical simulation are compared with those obtained from experiments. An uniaxial compression test on Chongqing sandstone sample C1 is conducted to obtain its uniaxial compressive strength (UCS) by the hydraulic controlled uniaxial loading system. The diameter and height of the cylindrical rock sample are respectively 50 mm and 100 mm. The UCS of Chongqing sandstone sample C1 is 103.684 MPa. Then, the triaxial compression experiments are conducted on Chongqing sandstone samples C2, C3 and C4 by the multifunction triaxial rock testing system 600–50 HT Plus. Compared with other triaxial rock testing system, this multi-function triaxial rock testing system adopts the adaptive confining pressure compensation technology developed by Université Lille 1, so that the influence of the piston movement on the confining pressure can be eliminated to ensure the stability of the confining pressure acting on the rock sample. During experiments, the maximum axial load is 1000 kN, and the maximum confining pressure is 60 MPa. In order to estimate the mechanical properties of Chongqing sandstone, namely, cohesion c and friction angleFðuÞ, at least three different kinds of loading stage are required to obtain the representative Mohr circles and the failure envelope of Chongqing sandstone. Then, uniaxial compressive strength (UCS) of Chongqing sandstone can be calculated by using the linear Mohr-Coulomb relationship. The full procedure of the triaxial compression test can be found in the International Society for Rock Mechanics (ISRM) suggested test method [41]. In this paper, three different confining pressures of 10, 20 and 30 MPa are applied in triaxial compression tests. The failure of Chongqing sandstone samples C2, C3 and C4 after tests are shown in Fig. 4(a)–(f). Young’s modulus of Chongqing sandstone sample can be calculated by the slope of different stages, Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

5

Fig. 4. The failure pattern of Chongqing sandstone sample after triaxial test: (a) View of sample C1. (b) View of sample C2. (c) View of sample C3 (d) View of sample C4.

and the Poisson’s ratio can be estimated by the ratio of radial strain to axial strain, as shown in Fig. 5. The cohesion u and the internal friction angle F s ðuÞ can be estimated by the average of the intercept and slope of the tangent line between the Mohr circles. The experimental results are listed in Table 1. The estimated properties of Chongqing sandstone are shown in Fig. 6 and Table 2. 3. Stochastic reconstruction 3.1. Control function Chongqing sandstone can be regarded as a two-phase medium. Therefore, naturally, after the Isodata threshold algorithm is used, CT images of Chongqing sandstone can be regarded as a two-phase digitized system, and the system is definitely defined for the solid phase and void phase by the indicator function as follows [27]:



vðxÞ ¼

1; x 2 void phase 0; x 2 solid phase

ð1Þ

where x is the position vector. In the 2D two-phase digitized system, x can take the discrete valuex ¼ ði; jÞ, i ¼ 0; . . . ; l1  1, j ¼ 0; . . . ; l2  1, which can represent the position of the voxel, and the size l1  l2 represents the resolution of the images. Similarly, in the 3D two-phase digitized system, the discrete values of the position vector x ¼ ði; j; kÞ are assigned as i ¼ 0; . . . ; l1  1, j ¼ 0; . . . ; l2  1, k ¼ 0; . . . ; l3  1, which can represent the position of the voxel, and the resolution of the 3D digitized system is defined by l1  l2  l3 . Then, the volume fraction of the void phase or porosity can be described as [27]:

u ¼ hvðxÞi

ð2Þ

Fig. 5. Stress-strain curve of Chongqing sandstone sample: (a) Axial stress-strain curve of Chongqing sandstone sample. (b) Axial stress versus radial strain curve of Chongqing sandstone sample.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

6

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Table 1 Triaxial compressive experimental results of Chongqing sandstone. Samples

Confining pressure (MPa)

Density (kN/m3)

Young’s modulus (GPa)

Poisson’s ratio

C2 C3 C4 Average

10 20 30

2376.13 2375.85 2451.44 2401.14

17.663 19.543 20.719 19.308

0.174 0.183 0.179 0.179

Fig. 6. The linear Mohr-Coulomb relationship.

where the angular brackets denote an ensemble average. In fact, Eq. (2) also defines the probability to be the void for one point in the digitized system, namely, the one-point probability function for the void phase in the two-phase system. The two-point probability function for the void phase can be written as [27]:

Sðx1 ; x2 Þ ¼ hvðx1 Þvðx2 Þi

ð3Þ

where x1 and x2 are two randomly selected points in the system. When the whole system is statistically heterogeneous and isotropic, the two-point probability function is only related to the distance u ¼ jx1  x2 j between two randomly selected point, naturally, Sðx1 ; x2 Þ ¼ SðuÞ. When u ¼ 0, SðuÞ ¼ u. When u ! 1, SðuÞ ¼ u2 . Under the same condition, SðuÞ  u2 ¼ Ss ðuÞ  u2s , in which Ss ðuÞ ¼ SðuÞ  2u þ 1 is the two-point probability function for the solid phase, us is the volume fraction of the solid phase. Therefore, in the two-phase system, the void phase and solid phase is not able to be distinguished by the two-point probability function. According to the definition of the one-point and two-point probability functions, we can know that this type of function provides the basic information about the statistical distribution of the void space in the system. Lineal-path function Lðx1 ; x2 Þ is another descriptive function for reconstructing the microstructure in stochastic reconstruction. The lineal-path function is defined as the probability of finding a line segment directed from x1 to x2 which entirely lays in the void phase. In a statistically heterogeneous and isotropic system, the lineal-path function is only affected by the distance u between two randomly selected points. Unlike SðuÞ, the lineal-path function can distinguish between the void phase and the solid phase. LðiÞ is defined as the probability of finding a line segment of length ði; i þ 1 which entirely lies in void phase. Considering that the void phase with a chord of length l is presented in a one-dimensional system, LðiÞ can be expressed as [27]:

 LðiÞ ¼

ðl  iÞ=N; when 0 6 i 6 l 0;

ð5Þ

otherwise;

where N is the voxel size of the system. For a system with more chords or with a higher dimension, its principle is the same as that in a system with fewer chords or with a lower dimension. For the solid phase, the lineal path function Ls ðuÞ can be defined by the same function for void phase. It should be noted that Ls ðuÞ–LðuÞ þ a, in which a is a constant. Compared with Table 2 The properties of Chongqing sandstone by triaxial compressive test. Cohesion (MPa)

Internal friction angle (°)

UCS (MPa)

35.6136

27

116.98

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

7

the two-point probability function, the lineal-path function provides the basic information about clustering or connecting characteristics of the microstructure. Besides this two functions mentioned above, we can introduce an arbitrary number of functions characterized the microstructure to the reconstruction process. The probable limitation of those functions on the proposed reconstruction method is the increase of computational resource. Aiming at completely and exhaustively reconstructing the rocks, considering the discontinuous and broken characteristics of rocks, the following fractal function is introduced to the reconstruction process [38]:

hvðx; y; zÞvðx þ 1; y; zÞ    vðx þ d; y; zÞ vðx; y þ 1; zÞvðx þ 1; y þ 1; zÞ    vðx þ d; y þ 1; zÞ FðuÞ ¼ FðxÞ ¼ . ..

ð6Þ

vðx; y þ d; zÞvðx þ 1; y þ d; zÞ    vðx þ d; y þ d; zÞi where d is the side lengths of the squares, vðxÞ is the probability of the point x in void phase, the product of vðxÞ determines whether or not all the points are in the white phase in one square with side d, the angular brackets denote the global average of the product of vðxÞ over the whole reconstruction system, and the value of FðuÞ or Fðu0 Þ denotes the percentage of the covering boxes with d completely filled by void voxels. The fractal function is intended to characterize the discontinuous and broken characteristics of microstructures during the reconstruction process. The corresponding function F s ðu0 Þ for the solid phase can be defined as the same as that for the void phase. For the evaluation of SðuÞ, LðuÞ and Fðu0 Þ, it may be helpful to assign an integer value to SðuÞ or u0 . 3.2. Random number generator During the reconstruction process, the system needs to be updated for thousands of times to obtain the optimal solution. Naturally, thousands of random numbers are needed. At present, the generated random number is pseudo-random one. The true random numbers are generated from physical phenomena. Generally, there are three types of pseudo-random number generator in C++ programming environment, namely, rand function, Mersenne Twister and CryptGenRandom. The rand function is a conventional random number generator included in C++ standard library based on linear congruential method. Although the numbers generated by rand function in a long period can be treated as random numbers, it does not meet the requirements of independence, homogenization, and randomness for reconstructing the porous rock model. The Mersenne twister is by far the most widely used general-purpose pseudorandom number generator, and has the advantages of long-period and k-distributed, which avoids the stratification in reconstructing porous rock model. But the Mersenne twister places a heavy load on the memory caches with the performance cost, and it is somehow slow under today’s standard. The CryptGenRandom is a cryptographically secure pseudo-random number generator function, which is included in Microsoft CryptoAPI and used in native TLS support and code signing. Although the CryptGenRandom is not still a real random number generator, its randomness and unpredictability are good enough to meet the requirements for reconstructing porous rocks. Moreover, the CryptGenRandom also avoids the stratification in reconstructing porous rock model when it is compared with rand function. Therefore, the CryptGenRandom is adopted as the default random number generator to generate the coordinate value of the exchanged pixel points. 3.3. The hybrid reconstruction algorithm and FE model generation The core ideology of the hybrid reconstruction algorithm originates from simulated annealing algorithm, which originates from the solid annealing phenomenon. During the solid annealing process, the solid is fully heated, the particles in the solid become disorder and internal energy increases, then it is slowly being cooled, the particles in the solid tend to be more orderly and the solid reaches equilibrium state at every temperature, finally the solid reaches ground state at normal atmospheric temperature and the internal energy decreases to the minimum. The modified simulated annealing algorithm employed in the reconstruction process can be described in Fig. 7. Simulated annealing (SA) algorithm is firstly proposed for combinatorial optimization problems [42]. In the proposed reconstruction algorithm, a high-energy state of void and solid voxels is slowly transformed into a low-energy configuration. The general objective function or control function of the whole digitized system can be defined as the sum of the square of the difference between an arbitrary number of reference control functions and the calculated control function determined by experiments. There are five control functions employed during our reconstruction process, namely, SðuÞ, LðuÞ, Ls ðuÞ, Fðu0 Þ and F s ðu0 Þ. The sum of the square of the difference between the reference control functions and the calculated control functions can be defined as follows:

e0 ¼

Ui n X X 2 ðe S i ðuÞ  Si ðuÞÞ

ð7Þ

i¼1 u¼0

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

8

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 7. The modified simulated annealing algorithm.

e1 ¼

Ui n X X 2 ðe L i ðuÞ  Li ðuÞÞ

ð8Þ

i¼1 u¼0

e2 ¼

Ui n X X 2 ðe L si ðuÞ  Lsi ðuÞÞ

ð9Þ

i¼1 u¼0

e3 ¼

Ui n X X 2 ðe F i ðu0 Þ  F i ðu0 ÞÞ

ð10Þ

i¼1 u¼0

e4 ¼

Ui n X X 2 ðe F si ðu0 Þ  F si ðu0 ÞÞ

ð11Þ

i¼1 u¼0

where the tilde e Sðu0 Þ denotes the reference control functions, U i is the maximum distance for each reference control functions along any direction i, n is the number of the directions with a uniform value in the control functions, which means that all the control functions can only be matched along the principal or the diagonal directions under periodic boundary condition, E can be obtained from the combination of the sums. Considering the importance of various control functions during the reconstruction process, a modified general control function C is defined as below:

C ¼ e0 þ e1 þ e3 þ

ðe2 þ e4 Þg g þ e0 þ e1 þ e2

ð12Þ

The specific weight of e2 and e4 is controlled by the positive constant g, which is only linear with time. During the annealing process, the minimization and maximization of the modified general control function, namely, u and T, lead to the formation of general control function by Eqs. (7), (8), (10) and Eqs. (7)–(11), respectively. At the beginning of the proposed reconstruction algorithm, void or pore and solid voxels are randomly generated to form the initial state of the reconstructed model with the constant porosity DE ¼ Ej  Ej1 . Then, a new state is randomly generated by selecting two voxels from different phase to exchange their phase function values. In the classical simulated annealing algorithm, an annealing step is finished by a pair of voxels from different phase, and the formation of the object blocks needs a slowly iterative annealing step with n, which consumes a lot of time. To accelerate the reconstruction process, a Voxel block exchange algorithm (VBE) is developed to form the preliminary reconstruction pore structure in this paper. Based on the two-point probability function, multiple voxels in one phase are selected from the whole system to exchange synchronously their phase function values for that of the same number of voxels in the other phase, as shown in Fig. 8. The exchanges change the corresponding control functions and result in the energy difference DE ¼ Ej  Ej1 . Such procedure to generate a new state is also referred as a Markov chain. The voxel exchange or the generated new state is accepted with probability n according to the Metropolis criterion, and the probability n can be written as [42]: Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

9

Fig. 8. The comparison between the original model and the preliminary model: (a) The random generated sandstone model. (b) The sandstone model after VBE algorithm.

( ðjÞ

nðDE Þ ¼

if DEðjÞ 6 0

1; eDE

ðjÞ

=T

; if DEðjÞ > 0

ð13Þ

where T is the control parameter of the system temperature. If the new state is rejected, the phase function values of those selected voxels remain the same, and the whole system keeps unchanged. The temperature drop is controlled by a threshold k, which means that the whole system is not cooled down until k new states are generated. The temperature control parameter k should be high enough for the whole system to turn into an equilibrium state. A cooling schedule is a series of process control parameters which force the system to return a satisfying optimum solution in a limited time frame, and is employed in the global optimization of simulated annealing algorithm. Generally, the cooling schedule is constituted of four parameters, namely, the initial temperature T 0 , the deduction rate function of the temperature f ðTÞ, the final temperature T e , and the length of the Markov chain Lk . In this paper, the ‘‘heat capacity” method [43] is introduced to the reconstruction process, and the target mean energy b E r can be calculated as[43]:

Er1  b E r ¼ v  xr1 Pk

ð14Þ

EðjÞ

where Er ¼ j¼1 is the calculated mean energy of the rth Markov chain, v is a constant parameter and controls the annealk ing speed, wr is the standard deviation of the energy distribution. Under the condition of b E 0 ¼ E0 , the temperature deduction can be estimated as below: 2

T T r  T r1 ¼ ð b Er  b E r1 Þ r1 2

xr1

ð15Þ

With the decrease of the temperature in the whole two-phase system, the corresponding energy decreases, and the success probability to exchange the phase function values of the randomly selected two voxels from different phase also decreases. That is to say, the characteristics of the reconstructed system become more like that of the target system. This results in the drop of the evolution rate of the system, which means that a lot of computational resources are consumed by the failure of exchanging the phase function values of the two voxels. Naturally, the generation of a new state becomes more and more difficult. In fact, the efficient generation of a new state needs to ensure that the randomly selected voxels are located in a different phase. In order to accelerate the evolution rate of the system and to generate new states more efficiently, two improved system evolution algorithms are proposed as follows: (1) Algorithm I When the failure time of the exchange reaches a threshold, an edge detection algorithm is employed to detect the interfaces of two different phases in the system. A starting voxel is chosen near the interface in the system. At the starting voxel, an artificial straight line propagates along arbitrary directions until it runs through the first interface; on the other side of the interface, the other point is randomly selected to exchange on the straight line before the straight line runs through the other interface or runs through the system. If there is no interface to be run through along the arbitrary directions, the straight line propagates along the opposite directions. Then, the same procedure is repeated. (2) Algorithm II Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

10

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

When the failure time of the exchange reaches a threshold, an edge detection algorithm is employed to detect the interfaces of two different phases in the system. A pair of starting points is randomly chosen near the interfaces in the system. For each pair of the point, an artificial straight line propagates from one point to the other point. If the artificial straight lines run through even pairs of interfaces, the selected voxels are exchanged with each other. If the steps or the assumptions in those algorithms fail (e.g. In algorithm I, there is no interface to be run through along the arbitrary or opposite directions, or in algorithm II, the artificial straight line runs through odd pairs of interfaces or no interface), the same procedure is repeated until there is a pair of selected voxels suitable for exchange. The corresponding phase function values are restricted to a certain range. It should be noted that the employment of those two algorithms is only effective under a proper threshold of failure times. The applications of improved system evolution algorithm I and II are plotted in Fig. 10(g) and (h), the final results seem the same except that algorithm I is faster and consumes much more computational resources than algorithm II. A complete reconstruction procedure proposed in this paper is described in Fig. 9, and is illustrated in Fig. 10. In general, the proposed reconstruction process can be divided into four stages, namely, the preliminary stage, the initial stage, the intermediate stage and the final stage. The values of the two-point probability function and the lineal-path function for reference and reconstructed model are plotted in Fig. 11(a) and (b). The values of the fractal dimension of the pore structures for the reference model and the reconstructed model are plotted in Fig. 11(c). The box count and box size of the reference and the reconstructed model are listed in Table 3. It is found from Fig. 11 that the values of the two-point probability function and the lineal-path function for the reconstructed model are very close to those for the reference model. Moreover, the values of the fractal dimension of the pore structures for the reconstructed model are also very close to those for the reference model. The similarities in function values indicate that the reconstructed model is almost equal to the reference model in statistics. Finally, the reconstructed 3D porous sandstone model and the reference system are developed by a novel geometrical reconstruction algorithm [39] for visualization and numerical analysis, as shown in Fig. 12. 4. Network extraction and topological comparison The network is the topological description of the pore space in the reconstructed model and the reference model, and provides a detailed statement about the characteristics of pore space in rocks or in the 3D reconstructed models. The network extracted from a 3D model includes pores which represent the large void space, and throats which define the narrow void space connecting the large void space [45]. There are two major algorithms to extract the network, namely, the medial axis based algorithm and the maximal ball algorithm. In this paper, the maximal ball algorithm is applied to extract the network, and this algorithm is modified by introducing the watershed algorithm. The modified maximal ball algorithms can be briefly described as follows: Step 1: A voxelized image from CT scanning is converted into a binary one by employing the OSTU threshold algorithm. Then, the watershed algorithm is employed to segment the void space into isolating small pore space, and the corresponding inscribed sphere is generated at each void voxel. Step 2: With the removal of the inclusions, the maximal balls are constructed and sorted and clustered into pore-throat chains using the following principle:

distðC 1 ; C 2 Þ < RRIGHT1 þ RRIGHT2

C1; C2 2 S

ð16Þ

where C 1 ; C 2 are the centers of the maximal balls, RRIGHT is the square of the Euclidean distance from the center of the maximal ball Cðxa ; ya ; za Þ to the nearest voxel Db ðxb ; yb ; zb Þ, and it can be written as:

Fig. 9. The proposed reconstruction algorithm.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

11

Fig. 10. The evolution process of the reconstructed model by the proposed method: (a) the raw CT images. (b) the Reference model after denoising and threshold algorithm. (c) The random generated sandstone model. (d) The preliminary model by employing VBE algorithm. (e) The initial state of the reconstructing model. (f) The intermediate state of the reconstructing model without the basic porous characteristics. (g) The intermediate state of the reconstructing model with the basic porous characteristics. (h) The final state of the reconstructed model by the improved system updating algorithm.

2

R2RIGHT ¼ dist ðC; Db Þ ¼ ðxb  xa Þ2 þ ðyb  ya Þ2 þ ðzb  za Þ2 ; C 2 S; Db 2 Sb

ð17Þ

where S and Sb represent the void and solid space in the reconstructed system, respectively. Step 3: Along the clustered pore-throat chains, the void space or pore space in the voxelized image is segmented into pores and throats, the corresponding topological parameters are calculated, such as size, volume, length and shape factor. Step 4: The pore network and the corresponding calculated parameters are exported, and the corresponding results are visualized. In fact, to some extent, the modified network extraction algorithm employed in this paper is similar to the maximal ball algorithm proposed by Hu and blunt [44], except that we employ the watershed algorithm to preliminarily segment the void space in step 1. This modified procedure greatly improves the efficiency of the whole network extraction algorithm. For more details of the maximal balls algorithm and its revised version, please refer to the works [44]. Although the distribution of Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

12

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 11. Comparison of the control function values between the reference and reconstructed models: (a) Two-point probability function value. (b) Linealpath function value (c) Fractal dimension value.

Table 3 Box count and box size of reference model and reconstructed model. Type

C2

C3

C4

C6

C8

C12

C16

C32

C64

D

Reference Reconstruction

23778 23669

11477 10909

6974 6825

3587 3351

2291 2182

1278 1169

858 807

306 300

98 92

1.5553 1.564

Fig. 12. 3D visualization of the reference and reconstructed model: (a) Visualization of reference model. (b) Visualization of reconstructed model.

Fig. 13. Network extraction: (a) Network of the reference model. (b) Network of the reconstruction model.

pore and throat in reconstructed model is different from that in the reference model in Fig. 13, the topological properties of the reconstructed model is almost the same as that of the reference model, which proves that the reconstructed model is equal to the reference model in statistics from the corresponding calculated parameters in Table 4.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

13

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 4 Calculated topological parameters of the reference model and reconstructed model.

Porosity Number of voxels Number of balls Number of throats Number of pores Invalid voxels

Reference model

Reconstructed model

18.5081% 11258200 8949362 261968 574629 378

18.3345% 11566040 9194067 244658 562837 224

5. Numerical simulation and verification of the reconstructed model 5.1. Theoretical basis for FE-SPH coupling method and failure criterion FE method is a conventional and efficient method to solve the partial differential equations with special boundary conditions in mechanics and mathematics. However, FE method has difficulties in dealing with the fracture propagation or large deformation under impact because its formulae are based on grids. Smoothed particle hydrodynamics method (SPH) is a mesh-free method, and nodes or elements are not needed in this method, instead, we just need to present an SPH model by a series of points called particles or pseudo-particles. The kernel concept of SPH method is to reproduce the physical field variables in an integral representation. SPH method obtains its adaptability in the preliminary stage of approximation for field variables, which has no effect on the construction of SPH formulae with the arbitrary particle distribution, so SPH method is naturally able to handle the problem of fracture propagation or large deformation under impact easily. Considering the complicated interactions in the simulation of the reconstructed porous rock model, a novel FE-SPH coupling method is introduced to investigate the stress-strain response of the reconstructed model under uniaxial compressive load. FE method is employed to calculate the stress and strain in the zone undergoing small deformation, while the SPH method is applied to simulate large deformation zone under compression. The FE grids or meshes within the large deformation zone are removed based on the threshold of conversion to particles and then are replaced by SPH particles located in the center of the deleted element. The generated particles inherit all the physical and mechanical properties of the deleted element including distance, quality, density, velocity, Poisson’s ratio and constitutive model. The calculation and conversion mode are plotted in Fig. 14. In this paper, the Maximum Stress Criterion (Rankine criterion), is used to predict the failure of the brittle rock materials. This criterion indicates that the failures appear when the tensile stress reaches the maximum tensile strength of rocks. In this criterion, cracks are considered as irrecoverable throughout the whole simulation, and the interactions among cracks are neglected. By introducing the condition of consistency, Rankine criterion can be written as

Fig. 14. The basic procedure of coupling FE-SPH method: (a) The initial establishment of the FE model. (b) Generation of SPH particles in the element. (c) The generated SPH particle inherits the mechanical properties of the original FE meshes. (d) Deletion of the original finite element. (e) Final stage of the SPH model.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

14

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

  C ¼ C t; rI;II ¼ 0

ð18Þ

where the superscript I; II respectively indicate that r is the direct tensile components of stress in a tension softening model (Mode I fracture) and rII is the shear component of stress in a shear softening/retention model (Mode II fracture), t represents r in the local coordinate system. In fact, the cracking of brittle materials is more complex than the classical plastic yielding of elastic-plastic material. In the cracking of brittle materials, the cracks may be in an active & opening status or a closing & reopening status, while the elastic-plastic material is only in a single plastic yielding status. For Mode I fracture and an active & opening status, Eq. (18) can be further developed as I

    C nn ¼ C nn tnn ; rIt ¼ tnn  rIt eck nn ¼ 0

ð19Þ

where rIt ðeck nn Þ is the tension softening evolution for an active opening crack, t nn represents r in the local coordinate system in a tension softening model. For a closing & reopening status, Eq. (18) can be obtained as open C nn ¼ C nn ðt nn ; rIc Þ ¼ t nn  rIc ðeck nn Þjenn ¼ 0

where

ð20Þ

open is the crack closing & reopening evolution, and depends on the maximum crack opening strain, rIc ðeck nn Þjenn

eopen nn

¼ maxov erhistory ðeck nn Þ. The cracking condition for Mode I fracture is plotted in Fig. 15(a). It should be noted that, although the cracking condition is stated for all preliminary existing cracks in general, only the components of C, which refer to the preliminary existing cracks, are considered in the computations. For Mode II fracture, Eq. (18) can be rewritten as

    ck ck C nt ¼ C nt tnt; rIIs ¼ t nt  rIIs g ck nt ; enn ; ett ¼ 0

ð21Þ

II ck ck ck s ðg nt ; enn ; ett Þ

where r is the shear evolution that linearly depends on the shear strain and the crack opening strain. The cracking condition for Mode II fracture is plotted in Fig. 15(b). 5.2. Numerical simulation and comparison Considering the huge number of small and inferior accompanying finite element meshes in the generation of finite element mesh, the combination of vertex welding algorithm with Temporal Branch-on-Need Octree (T-BON) algorithm [39] is introduced to remove the duplication without any loss of detail. In this paper, the Maximum Stress criterion is used to calculate the strength of Chongqing sandstone. The physical and mechanical properties obtained from experiments should not be directly assigned to the reconstructed model because it leads to the second deduction in properties. Considering this, based on Chongqing sandstone experiments, the finite element method (FEM) is employed to calculate the proper parameters in the reference model. Then the reconstructed model is simulated under different confining pressure with the proper parameters. The calculated parameters used in the numerical simulation are listed in Tables 1 and 2. The reconstructed model after the vertex welding and the T-BON algorithm is shown in Fig. 16(a). Chongqing sandstone is considered as a two-phase material. Although the pore phase may contain some lowstrength bonding material, the pore structure in the reconstructed model is set to its ultimate state as void space, and the properties of the pore are set as null rather than a low value. At the same time, considering the influence of structure evolution on the global stiffness, the flexibility matrix is calculated based on the displacement of the elements. Considering the uneven surface, the reconstructed model and reference model cannot be loaded and restricted conventionally. Therefore, we use rigid plates to load the displacement. Displacements and rotations in the lower plate are fixed along all direction. Displacements in the upper plate are fixed along the direction of X- and Z-axis, and rotations in the upper plate are fixed along all direction, as shown in Fig. 16(b). The confining pressure is also applied to the FE model by the rigid plates. The upper plate

Fig. 15. The cracking condition: (a) The cracking condition for Mode I fracture. (b) The cracking condition for Mode II fracture.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

15

Fig. 16. The reconstructed FE model and its constraint condition: (a) Visualization of the reconstructed FE model. (b) Constraint condition of the reconstructed FE model.

is loaded by displacement along the direction of Y-axis, as shown in Fig. 17. The constraint condition of the reconstructed model is the same as that of the reference model. The maximum principal stress contour under confining pressure of 10 MPa is plotted in Fig. 18. It is found from Fig. 18 that at the microscopic scale, the distribution of stress is uneven and stress concentration appears in the local discontinuous area due to the pore structure. This may have a significant influence on the crack growth in rocks. Although the distributions of stress are different between those two models, the numerical physical and mechanical responses are almost the same in those two models. The fractures in rock materials are largely caused by native micro-defects, such as native micropore structures. Therefore, the numerical simulation results concretized by the ABAQUS code have an important meaning on the study of the rock fracture mechanics under simulating annealing reconstruction. Fig. 19 presents the crack propagation process during the compressive test, and it also indicates that the whole model experiences the crack initiation stage, the crack propagation stage, and the crack coalescence stage which finally leads to the failure of the model. Those figures are plotted from the cross-section of the numerical results in Fig. 18. As mentioned above, the FEM elements are converted into SPH particles with a critical degree of deformation when the rock is fractured to form cracks, so the location and area of crack propagation can be determined by the FEM elements which are converted into the SPH particles. In this paper, in order to highlight the differences for different crack propagation states in the numerical results, the SPH particles are enlarged to indicate the crack propagation path, and the overlaid SPH particles are hidden. At the beginning of the test, with the loading of the rigid plate on the rock sample, a number of tensile and compressive failure elements are observed on the tip edge of the cracks during the crack initiation stage, as shown in Fig. 19(b). Hertzian cracks are initiated at the top and bottom of the model, and some thin cracked zones form underneath the loading rigid plate. The whole rock model remains relatively intact owing to the high hydrostatic pressure [46–48]. With the increase of the loading displacement, the Hertzian cracks grow into micro-cracks in the cracked zone during the crack propagation stage, as shown in Fig. 19(c) and (d), and the thin cracked zones also expand horizontally and are connected with each other to form a larger thin cracked zone. After the thin cracked zone expands to a certain degree, the expansion and propagation of the lateral cracks slow down, while the median cracks start to propagate along the loading direction.

Fig. 17. Comparison of the displacement contour between reconstructed model and the reference model: (a) Displacement contour of the reference model. (b) Displacement contour of the reconstructed model.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

16

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 18. The comparison of the maximum principle stress contour between reconstructed model and the reference model: (a) The maximum principle stress contour of the reference model. (b) The maximum principle stress contour of the reconstructed model.

Fig. 19. Crack propagation process: (a) The rock sample before compressive test. (b) the crack initiation stage, (c) and (d) the crack propagation stage, (e) and (f) the crack coalescence stage.

With a further increase of the loading displacement, as shown in Fig. 19(e) and (f), both the thin cracked zone and the lateral cracks remain slowly changing during the crack coalescence stage, while the median cracks keep on propagating along the loading direction. As a result, the median cracks are connected with each other, which finally leads to the failure of the whole model.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

17

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

In the compressive test, the stress field undergoes a dynamic expansion process, as shown in Fig. 19(b)–(f). Moreover, the failure elements are more inclined to be initially observed near the pore structures and the cracks are more inclined to propagate along the pore structures in the numerical model. This indicates that the pore structures have a very important influence on the initiation and propagation of the cracks. As shown in Fig. 20, the axial stress-strain responses of the reconstructed model are very close to those of the reference model and experimental data under confining pressure of 20 MPa. Although the reference model is accurately reconstructed based on the CT images data by a series of mathematical procedures, and can also be used in the numerical simulation to analyze the basic characteristic of porous rocks, the size of the reference model is limited to the scanning area (mm) and the features of porous rocks are only partially reflected. Therefore, the reconstructed model is developed. Compared with the reference model, the reconstructed model is able to more comprehensively characterize the porous rock in various pore structures configurations and to predict the mechanical properties of porous rocks in a larger size based on the statistical control function. Moreover, as shown in Fig. 21(a), the numerical results of the axial stress-strain of the reconstructed model are in good agreement with experimental data under different confining pressure. The numerical results of peak strength of Chongqing sandstone are also in good agreement with experimental data, as shown in Fig. 21(b). Table 5 shows the numerical param-

Fig. 20. Comparison of the numerical results of the reconstructed and reference model with experimental data under confining pressure of 20 MPa.

Fig. 21. Comparison of the numerical and experimental results under different confining pressure (a) Numerical results of the axial stress-strain curves of the reconstructed model under different confining pressure (b) Numerical and Experimental results of peak strength of Chongqing sandstone sample.

Table 5 The numerical parameters calculated by the reconstructed model. Density (kg/m3)

Young’s modulus (GPa)

Poisson’s ratio

UCS (MPa)

3271.49

24.015

0.2664

116.98

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

18

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

eters calculated by the reconstructed model. Compared with Tables 1 and 2, it is found that the numerical parameters obtained by the reconstructed model are close to experimental data. All the results imply that the reconstructed model has the ability to safely predict mechanical properties of porous sandstone under the different configuration of pore structures with limited experimental information. Therefore, the proposed method provides the possibility to predict the macroscopic mechanical behaviors of porous rocks. 6. Conclusions In this paper, a novel 3D reconstruction method is proposed to generate the 3D reconstruction model of porous rocks based on two-point correlation function, lineal-path function, and fractal function. Firstly, the 3D porous rock model is reconstructed based on the 3D reference model generated by CT scanning of Chongqing sandstone sample. The reconstruction model and reference model are compared in topological parameters and mechanical properties to validate the proposed method. The comparison results showed that those two models are equal in topology and mechanical properties. It is demonstrated that the proposed method can exactly describe the complexity of the pore structures. Then, the reconstructed model is simulated by Finite Element-Smoothed Particle Hydrodynamics (FE-SPH) coupling method under compressive loads. The material behaviors of the reconstructed model are simulated by using a brittle material model formulation. The numerical results show that at the microscopic scale, the reconstructed porous rock model goes through crack initiation stage, crack propagation stage and crack coalescence stage under compressive loads. Finally, the comparison of the numerical results and experimental data shows a good agreement, which indicates that the proposed reconstruction method has a good ability to predict the macroscopic behaviors with limited experimental data. Acknowledgements The work is supported by the National Natural Science Foundation of China (Nos. 51325903, 51279218), project 973 (Grant no. 2014CB046903) and Natural Science Foundation Project of CQ CSTC (No. CSTC, cstc2013kjrc-ljrccj0001). References [1] Lymberopoulos DP, Payatakes AC. Derivation of topological, geometrical, and correlational properties of porous media from pore-chart analysis of serial section data. J Colloid Interface Sci 1992;150(1):61–80. [2] Vogel HJ, Roth K. Quantitative morphology and network representation of soil pore structure[J]. Adv Water Resour 2001;24(3–4):233–42. [3] Tomutsa L, Silin D. Nanoscale pore imaging and pore scale fluid flow modeling in chalk. Office of Scientific & Technical Information Technical Reports; 2005. [4] Fredrich JT, Menéndez B. Imaging the pore structure of geomaterials. Science 1995;268(5208):276. [5] Zhou S, Yan G, Xue H, et al. 2D and 3D nanopore characterization of gas shale in Longmaxi formation based on FIB-SEM. Mar Pet Geol 2016;73:174–80. [6] Dunsmuir JH, Ferguson SR, D’Amico KL, et al. X-Ray microtomography: a new tool for the characterization of porous media. Spe Technical Conference and Exhibition; 1991. [7] Coenen J, Tchouparova E, Jing X. Measurement parameters and resolution aspects of micro X-ray tomography for advanced core analysis. In: Proceedings of International Symposium of the Society of Core Analysts, Abu Dhabi, UAE; 2004. [8] Yu W, Jie P, Wang L, et al. Characterization of typical 3D pore networks of Jiulaodong formation shale using nano-transmission X-ray microscopy. Fuel 2016;170:84–91. [9] Lorensen WE, Cline HE. Marching cubes: a high-resolution 3D surface construction algorithm. ACM Siggraph Comput Graph 1987;21(4):163–9. [10] Vignoles GL, Donias M, Mulat C, et al. Simplified marching cubes: an efficient discretization scheme for simulations of deposition/ablation in complex media. Comput Mater Sci 2011;50(3):893–902. [11] Delibasis KS, Matsopoulos GK, Mouravliansky NA, et al. A novel and efficient implementation of the marching cubes algorithm. Comput Med Imaging Graph 2001;25(4):343–52. [12] Rajon DA, Bolch WE. Marching cube algorithm: review and trilinear interpolation adaptation for image-based dosimetric models. Comput Med Imaging Graph 2003;27:411–35. [13] Chew LP. Guaranteed-quality triangular meshes. USA: Cornell University; 1989. [14] Danlov AA. Unstructured tetrahedral mesh generation technology. Comput Math Math Phys 2010;50(1):139–56. [15] Cheng SW, Poon SH. Three-dimensional Delaunay mesh generation. Discr Comput Geom 2006;36:419–56. [16] Du Q, Wang DS. Recent progress in robust and quality Delaunay mesh generation. J Comput Appl Math 2006;195:8–23. [17] Bakke S, Ren PE. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. SPE J 1997;2(2):136–49. [18] Øren PE, Bakke S. Process based reconstruction of sandstones and prediction of transport properties. Transp Porous Media 2002;46(2–3):311–43. [19] Joshi MY. A class of stochastic models for porous media. University of Kansas. Chemical and Petroleum Engineering; 1974. [20] Adler PM, Jacquin CG, Quiblier JA. Flow in simulated porous media. Int J Multiph Flow 1990;16(4):691–712. [21] Levitz P. Off-lattice reconstruction of porous media: critical evaluation, geometrical confinement and molecular transport. Adv Coll Interface Sci 1998; s76–77(2):71–106. [22] Strebelle S. Conditional simulation of complex geological structures using multiple-point statistics. Math Geol 2002;34(1):1–21. [23] Arpat GB, Caers J. Conditional Simulation with Patterns. Math Geol 2007;39(2):177–203. [24] Mariethoz G, Caers J. Multiple-point geostatistics: stochastic modeling with training images; 2014. [25] Hazlett RD. Simulation of capillary-dominated displacements in microtomographic images of reservoir rocks. Transp Porous Media 1995;20(1– 2):21–35. [26] Rintoul MD, Torquato S. Reconstruction of the structure of dispersions. J Colloid Interface Sci 1997;186(2):467–76. [27] Yeong CLY, Torquato S. Reconstructing random media. Phys Rev E 1998;57:495–506. [28] Yeong CLY, Torquato S. Reconstructing random media. II. Three-dimensional media from two-dimensional cuts. Phys Rev E 1998;58:224–33. [29] Torquato S. Random heterogeneous materials: microstructure and macroscopic properties. New York: Springer; 2002. [30] Sheehan N, Torquato S. Generating microstructures with specified correlation functions. J Appl Phys 2001;89(1):53–60. [31] Ekkehart Kröner. Statistical continuum mechanics. Vienna: Springer; 1971. [32] Kröner E. Bounds for effective elastic moduli of disordered materials. J Mech Phys Solids 1977;25(2):137–55. [33] Beran M. Statistical continuum theories. Phys Today 1969;22(7). 94-94.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021

X.-P. Zhou, N. Xiao / Engineering Fracture Mechanics xxx (2017) xxx–xxx

19

[34] Manwart C, Hilfer R. Reconstruction of random media using Monte Carlo methods. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 1999;59(5 Pt B):5596–9. [35] Jiao Y, Stillinger FH, Torquato S. Modeling heterogeneous materials via two-point correlation functions: basic principles. Phys Rev E: Stat, Nonlin, Soft Matter Phys 2007;76(3 Pt 1):031110. [36] Jiao Y, Stillinger FH, Torquato S. Modeling heterogeneous materials via two-point correlation functions. II. Algorithmic details and applications. Phys Rev E: Stat, Nonlin, Soft Matter Phys 2008;77(1):031135. [37] Davis MA, Walsh SD, Saar MO. Statistically reconstructing continuous isotropic and anisotropic two-phase media while preserving macroscopic material properties. Phys Rev E: Stat, Nonlin, Soft Matter Phys 2011;83(2 Pt 2):026706. [38] Ju Y, Zheng J, Epstein M, et al. 3D numerical reconstruction of well-connected porous structure of rock using fractal algorithms. Comput Methods Appl Mech Engng 2014;279(9):212–26. [39] Zhou XP, Xiao N. A novel 3d geometrical reconstruction model for porous rocks. Engng Geol 2017;288(13):371–84. [40] Ridler TW, Calvard S. Picture thresholding using an iterative selection method. IEEE Trans Syst Man Cybern 2007;8(8):630–2. [41] Franklin JA. Suggested methods for determining the strength of rock materials in triaxial compression: revised version. Int J Rock Mech Min Sci Geomech Abstr 1978;20(6):285–90. [42] Aarts E, Korst J. Simulated annealing and Boltzmann machines: a stochastic approach to combinatorial optimization and neural computing. Chichester: Wiley; 1989. [43] Frost R, Heineman P. Simulated annealing: a heuristic for parallel stochastic optimization. In: Arabnia HR, editor, Proceedings of the International Conference on Parallel and Distributed Processing Techniques and Applications, PDPTA 1997; June 30–July 3, 1997, Las Vegas, Nevada, USA, pp. 1595– 1564. CSREA Press. [44] Dong H, Blunt MJ. Pore-network extraction from micro-computerized-tomography images.[J]. Phys Rev E Stat Nonlinear Soft Matter Phys 2009;80 (2):036307. [45] Wilhelms J, Van Gelder A. Octrees for faster isosurface generation. The Workshop on Volume Visualization. ACM 1990; 57–62. [46] Cook NGW, Hood M, Tsai F. Observations of crack growth in hard rock loaded by an indenter. Int J Rock Mech Min Sci Geomech Abstr 1984;21 (2):97–107. [47] Chiaia B. Fracture mechanisms induced in a brittle material by a hard cutting indenter. Int J Solids Struct 2001;38(44–45):7747–68. [48] Liu HY, Kou SQ, Lindqvist PA, et al. Numerical simulation of the rock fragmentation process induced by indenters. Int J Rock Mech Min Sci 2002;39 (4):491–505.

Please cite this article in press as: Zhou X-P, Xiao N. Analyzing fracture properties of the 3D reconstructed model of porous rocks. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.021