A novel 3D geometrical reconstruction model for porous rocks

Engineering Geology 228 (2017) 371–384

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

A novel 3D geometrical reconstruction model for porous rocks Xiao-Ping Zhou a b

a,b,⁎

, Nan Xiao

MARK

a

School of Civil Engineering, Chongqing University, Chongqing 400045, China State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400045, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Pore-matrix interface extraction The surface reconstruction, 3D numerical analysis Delaunay refinement algorithm Porous media

In this paper, a novel geometrical reconstruction model is proposed for numerical simulation of pore structure in rocks. The spatial configuration and geometrical information of the microstructure are firstly gathered from Chongqing sandstone by X-ray micro-CT observations, and experimental tests on Chongqing sandstone are applied to obtain the mechanical parameters. Then, on the basis of the gathered geometrical information, the surface geometrical model is established by using the modified Marching cubes algorithm. Next, the established surface geometrical model is extended to 3D tetrahedral model by using Delaunay refinement algorithm. Finally, the 3D tetrahedral model is transformed into finite element (FE) numerical model. Moreover, the relationship between the porosity, fractal dimension of pore structure and mechanical properties of Chongqing sandstone is also investigated. The comparison of the numerical results and the experimental data shows a good consistency.

1. Introduction Rock naturally contains a great many discrete, multiple scales and irregular-in-shaped pores (Cueto et al., 2009), the microstructure in rocks has not only a direct relation to the porosity and permeability of rocks, but also has an indirect relation to Young's modulus, Poisson's ratio, and uniaxial compression strength, etc. And the schedule or evaluation of the engineering project depends on the measurement of the mechanical properties of rocks. However, in the conventional numerical modeling of rocks, the microstructures in rocks are often eliminated or simplified. This simplified measure more or less results in a simple investigation in the mechanical behaviors of rocks. Thus, a more accurate geometrical model, which can reflect the true features of rocks and predict its mechanical behaviors of rocks more accurately, is required. Various methods are employed to investigate the microstructures in rocks, such as serial sectioning method (Schlueter and Zimmerman, 1997), confocal laser scanning microscopy method (Uday et al., 2013), FIB-SEM method (Zhou et al., 2016), micro-CT method (Sarkar and Siddiqua, 2016) and nano-transmission X-ray microscopy method (Yu et al., 2016). Compared with other direct imaging methods, the micro-CT or X-ray CT method is a non-destructive method and becomes an effective method in engineering geology. Sarkar and Siddiqua (2016) used X-ray computed tomography (X-ray CT) to study the light backfill under the distilled water and two other pore fluid conditions, and explored the effect of the pore fluid on the microstructure, pore geometry and hydraulic conductivity of a barrier material. Zhou et al. (2015) introduced a mathematical procedure using ⁎

spherical harmonics to characterize and reconstruct the particle micromorphology in three dimensions based on micro X-ray CT data. Lanzón et al. (2014) performed X-ray spectrometry and thermogravimetric analysis on the chemical and mineralogical characterization of the calcarenitic stone quarried from ancient time in Cartagena, Spain. Yun et al. (2013) presented an assessment of anisotropy in rock by systematic and planar clustering of three-dimensional X-ray attenuation values and associated statistical analysis on the basis that X-ray CT numbers directly represent the internal density of the mineral fabric in the rock mass. Based on XRCT data, Christe et al. (2011) tested a new index for evaluating the degree of alteration of highly tectonized carbonate rock samples and compared the new index with visual-based classification systems. Then, the 3D geometry model can be reconstructed to provide a fairly accurate geometrical description of the microstructures in the rock materials by a series of mathematical procedures, such as Marching cubes algorithm and Delaunay algorithm (Lorensen, 1987; Bowyer, 1981). However, due to numerous small void space inside rocks and the random distribution of the microstructures, there are still some differences between the natural rocks and the reconstructed model by those method mentioned above, which make the prediction results unsatisfied. Therefore, some models were proposed to establish the relationship between microstructure and porous material properties, such as the elastic spider network model (Gent and Thomas, 1959), cube structure model (Gent and Thomas, 1963) and GibsonAshby model (Gibson and Ashby, 1997). Those models are applied to simulate the internal microstructures of porous materials based on mechanical consideration or morphological consideration. However,

Corresponding author at: School of Civil Engineering, Chongqing University, Chongqing 400045, China. E-mail address: [email protected] (X.-P. Zhou).

http://dx.doi.org/10.1016/j.enggeo.2017.08.021 Received 2 May 2017; Received in revised form 18 August 2017; Accepted 21 August 2017 Available online 01 September 2017 0013-7952/ © 2017 Elsevier B.V. All rights reserved.

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

C2, C3, C4. The rock samples are scanned by SOMATOM scope to obtain the CT images, as shown in Fig. 2(a) and (b).

those models are still somehow simplified, and therefore it results in the calculation error when the thickness of the wall of the pore is uneven or the surface of the microstructure is rough. In this paper, the SOMATOM scope is employed to generate X-ray computed tomograph (CT) images from Chongqing sandstone samples. Then, a series of triaxial compression experiments on Chongqing sandstone samples are conducted to obtain the mechanical properties by the Multi-function triaxial rock testing system. The mathematical reconstruction procedures are modified and a novel reconstruction model is proposed. The relationship between the microstructure and properties of porous rocks is investigated based on the proposed reconstruction algorithm. We tried to improve and guarantee the efficiency and accuracy of the reconstructed model by a series of mathematical derivations and deductions. First, the image data obtained from CT scanning is processed by a denoising algorithm based on a nonlinear total variation (TV) algorithm. Then, the surface model is reconstructed by a modified Marching cubes algorithm based on the processed image data. Next, the 3D geometrical model is reconstructed by using the improved Delaunay refinement algorithm and employing the surface model as the input data. Finally, the 3D reconstructed model is simulated under compressive loads, and the numerical results are compared with the experimental data. The relationship between the microstructure and properties of porous rocks is also investigated.

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 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 experiment results are shown in Fig. 4(a). 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 rock samples C2, C3 and C4 by the multi-function triaxial rock testing system 600–50 HT Plus, as shown in Fig. 3. 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 angle c ∈ IR, at least three different kinds of loading stages 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. In this paper, three different confining pressures of 10 MPa, 20 MPa and 30 MPa are respectively applied to Chongqing sandstone sample C2, C3 and C4. The failure of Chongqing sandstone sample C2, C3 and C4 after tests are shown in Fig. 4(b)–(d). Young's modulus of Chongqing sandstone can be calculated by the slope of elastic stage, and Poisson's ratio can be estimated by the ratio of radial strain to axial strain, as shown in Fig. 5. The cohesion c and the internal friction angle c ∈ IR 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 given in Fig. 6 and Table 2.

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 summarized in Fig. 1(b). During the working process of the CT device, the X-ray source continues emitting X-ray beam which can penetrate the sample placed on the sample stages, then 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. During this process, the detector and receiver are rotated to obtain raw data or image data from various directions. 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

3. The basic principles of the reconstruction method and preliminary process In the proposed reconstruction model, the gathered data are mainly derived from the images generated by the CT scanning. Furthermore, 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.

372

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 2. Chongqing sandstone sample. (a) The raw CT images of samples C2 and C3. (b) The raw CT images of sample C4. Fig. 3. The multi-function triaxial rock testing system. (a) The operating console of the multi-function triaxial rock testing system. (b) The loading platform of the multi-function triaxial rock testing system.

important and necessary preliminary procedure. Compared with other denoising algorithms, the denoising algorithm based on nonlinear total variation (TV) (Jia and Zhao, 2010) has a better performance on keeping the details of visual effect like edge or corner details. In this paper, the shock filter is combined with anisotropic diffusion in developing an adaptive total variation model for image denoising, then the edge detection operator is employed to distinguish the edge from the smooth area to obtain the adaptive parameters which control the diffusion level of the total variation model, and the sensitivity of the adaptive parameters to the noise is decreased so as to apply the morphology principle of image processing to the denoising model. Moreover, both of the image smoothing denoising and preservation of edge or corner details are taken into consideration while the piecewise constant effect is prevented. The general form of TV algorithm is shown below (Grasmair, 2009):

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

f=u+n the gathered data is the pixel on the grayscale images from the CT scanning, and the pixel is the smallest control element of the images. The values of the pixel in the Grey Scale images are directly proportional to the attenuation of the electromagnetic wave, which varies from 0 (black) to 255 (white). Therefore, we can set up the well-determined geometry model of Chongqing sandstone by the volumetric data from CT scanning. However, during the formation and transmission of the image signal, the noise disturbance may lower the quality of the image, and leads to the imprecisely modeling of the microstructure in rocks. Before employing the reconstruction algorithm, image denoising is an

(1)

where u represents the original image, f represents the degraded image, n represents the Gaussian white noise with zero means and δ is the standard deviation. The energy functional can be written as follows:

TV Fu (u) = arg min TV(u) + u

μ 2

∫ Ω |f (x ) − u (x )| 2 dx

(2)

where μ > 0 represents regularization parameter, TV(u) represents the total variation of the restored images u, Ω represents the image space and x ∈ Ω. The total variation TV(u) in Eq. (2) is anisotropic, namely, TV (u) = |∇xu| + |∇yu|, and the method noise corresponding to Eq. (2) can be rewritten as follows: 373

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

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.

with the anisotropic diffusion (Alvarez and Mazorra, 1994) in the following form:

Table 1 Test results of Chongqing sandstone sample. Sample no.

C2 C3 C4 Average

Confining stress (MPa)

Density (kN/ m3)

Young's modulus (GPa)

Poisson's ratio

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

∂f = −sign(Gσ ⊗ fNN ) sign(Gσ ⊗ fN )|∇f | + c1 fTT ∂t

(5)

where ⊗ represents convolution operator, Gσ represents a Gaussian function in which the standard deviation is σ , c1 represents a positive constant, fNN and fTT represent the gradient of the evolution curve and the second derivative along the tangent direction, respectively. Four edge detection operators dθ (Almeida and Almeida, 2010) are introduced to calculate the adaptive parameter ξf. When θ = 0, O1 ⎤ 1 ⎡ d 0 = 12 •⎢O2 M O2⎥, M = ⎡ 1 2 2 1 ⎤, O1 = ⎡ 0 0 0 0 0 0 ⎤, O2 = ⎡ 0 ⎤, ⎣ − 1 − 2 − 2 − 1⎦ ⎣0 0 0 0 0 0⎦ ⎣0⎦ ⎢ ⎣ O1 ⎥ ⎦ the rest of the edge detection operators can be obtained by rotating the angle θ ∈ Θ(Θ = {450, 900, 1350}). Based on the definition of the edge detection operators, the adaptive parameter ξf can be written as follows:

ξf =

1 1+

f )2 ∑θ ∈ Θ (dθ ⊗ 

∈ (0, 1),

Θ = {0°, 45°, 90°, 135°} (6)

and the method noise corresponding to Eq. (2) can be written as:

ξf u (x ) − Fξ , μ (u)(x ) = − curv(Fξ , μ (u))(x ) μ

(7)

 f )2

→ ∞, ξf → 0. Under In the edge area of the images, ∑θ ∈ Θ (dθ ⊗ the condition of μ > 0, u(x) − Fξ , μ(u)(x)→0, u(x) − Fμ(u)(x)→ − ∞. Therefore, when the images are denoised, the traditional total variation model may blur the edges, but the adaptive total variation model can preserve the edge details. The Bregman algorithm (Goldstein and Osher, 2009) is introduced to quickly solve Eq. (4) and to avoid the piecewise constant effect. For κ > 0 and c1 ∈ IR, the operator cut is defined as

Fig. 6. The linear Mohr-Coulomb relationship.

Table 2 Estimated properties of Chongqing sandstone sample. Cohesion (MPa)

Internal friction angle (°)

UCS (MPa)

35.6136

27

116.98

1 c 1 cut ⎛c1 , ⎞ ≔ c1 − 1 •max ⎛|c1| − , 0⎞ |c1 | κ ⎠ ⎝ κ⎠ ⎝

1 u (x) − TV Fμ (u)(x ) = − curv(TV Fμ (u))(x ) μ

(3)

where curv(TV Fμ(u))(x) represents the curvature of all the level set. To denoise images,μ should be diminished, which leads to the loss of edge or corner details. In order to preserve the edge or corner details while the images are denoised, an adaptive parameter named ξf is introduced, and the adaptive total variation model can be obtained as:

μ Fξ , μ (u) = arg min ξf TV(u) + u 2

∫Ω

|f (x ) − u(x )|2 dx

(8)

C = ξf,λ > 0 is the Assuming u Auxiliary variable, Eq. (4) can be equivalent to the fast iteration algorithm in the following form: 1

1 bxk ≔ cut ⎛C ∇x uk + bxk − 1 , ⎞ λ⎠ ⎝

(9)

1 byk ≔ cut ⎛C ∇y uk + byk − 1 , ⎞ λ⎠ ⎝

(10)

λ T k (∇x bx + ∇Ty byk ) μ

(11)

uk + 1 ≔ f −

(4)

= f , bx = by0 = 0 (k = 1, 2, …), 0

Jia and Zhao (2010) proved that for k = 0 , 1 , …, there is lim uk = u∗ when 0 < λ  μ < 1 8. Because  μ = μ C , C ∈ (0, 1], we have

To reduce the sensitivity of ξf to noise, the shock filter is combined

k →∞

374

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

(a)

(b)

(c)

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

0<λ  μ < 1 8 ⇒ 0 < C λ μ < 1 8. For the arbitrary C ∈ (0, 1], 0 < λ μ < 1 8 ⇒ 0 < C•λ μ < 1 8 ⇒ lim uk = u∗ can be determined, k →∞

namely, the optimal solution of Eq. (4). The denoising algorithm based on nonlinear total variation (TV) is employed to denoise the CT images, as shown in Fig. 7(a). Then, the Isodata threshold algorithm (Ridler and Calvard, 1978) is employed to segment the CT images shown in Fig. 7(b) to obtain the images suitable for numerical reconstruction depicted in Fig. 7(c).

4. Surface model based on modified marching cubes algorithm Fig. 9. Ambiguity on the volume.

After the image data are processed preliminarily by using TV denoising algorithm and segmented by using the Isodata threshold algorithm, the redundant disturbance is eliminated, and the geometrical model of rocks is able to be accurately reconstructed based on the processed image data. In this paper, the surface model is adopted as the intermediate state from the processed image data to the geometrical model suitable for investigation and simulation of mechanical properties of rocks. There are several surface reconstruction methods to generate a surface model, such as Marching cubes algorithm (Lorensen and Cline, 1987), Cuberille algorithm (Herman and Liu, 2010) and Dividing Cubes algorithm (Cline and Lorensen, 1988), etc. After the Classic Marching Cubes algorithm was firstly proposed by Lorensen, some defects were discovered until now. The main defect of the MC algorithm is topological ambiguity, which leads to the appearance of holes in the constructed isosurface. For example, when the vertices on the surface of the volume cell are arranged in Fig. 8, there may be an alternative choice for the connection modes among the dots with the pixel of intermediate values, then a hole is generated. This is indicated as Type A “hole problem” (Zhou et al., 1994), which is the ambiguity of MC algorithm on the surface. When the vertices on the volume cell are arranged in Fig. 9, there may also be another alternative choice for the connection modes among the dots with the pixel of intermediate values, which is the ambiguity of MC algorithm on the volume. Although this mode does not result in Type A “hole problem”, the mesh generated by this mode may be an error. The other weakness of

MC algorithm is that a huge number of small and inferior triangle meshes are generated and used to construct isosurface so that the construction speed of isosurface is slow when the MC algorithm is used to process a large number of data. There are many methods, which are developed to improve MC algorithm, based on either correction of the ambiguity (Lewiner et al., 2003) or reduction on the huge number of small and inferior triangle meshes (Vignoles et al., 2011), but it is not still up to expectations for the easiness in realization. To solve the Type A “hole problem” and the problem of the ambiguity on volume, a trilinear interpolate T is introduced (Nielson, 2003) as:

T (x , y, z ) = a2 xyz + b2 xy + c2 yz + d2 zx + e2 x + f2 y + g2 z + h2

(12)

The saddles points are limited to each face of the cube, for example, for the face x = 0, T(0, y, z) can be written as

T (0, y, z ) = c2 yz + f2 y + g2 z + h2

(13)

Assuming ∂ T/∂ z = ∂T/∂z = 0, we can find the saddle point at (0, − f2/c2, − g2/c2) with a saddle value of h2 − f2g2/c2. The body saddle point is the solution to the system of Eqs. (14), (15), and (16).

∂T ∂x = a2 yz + b2 y + d2 z + e2 = 0

(14)

∂T ∂y = a2 xz + b2 x + c2 z + f2 = 0

(15)

Fig. 8. Ambiguity on the surface.

375

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 10. Removal of the redundant vertices.

∂T ∂z = a2 xy + c2 y + d2 x + g2 = 0

(16)

interpolation at the saddle point. Similar to the saddle point of the bilinear interpolation in two dimension, the saddle point B(xs,ys,zs) of the trilinear interpolation is the transition point that decides which pair of vertices are disconnected in three dimensions, and the ambiguity on the volume can be eliminated by comparing the threshold with the value of the interpolant at the saddle point B(xs, ys, zs). It is worth noting that conditions 3, 4, 6, 7, 10, 12 and 13 can be respectively disintegrated into 2, 2, 3, 5, 3, 4 and 7 kinds of new conditions by this method. Therefore, the classic 15 kinds of conditions can be extended to 33 kinds of conditions, and the corresponding edgeTable and triTable should be relatively modified. To solve another weakness of MC algorithm, namely, the slow processing speed and the huge number of small and inferior triangle meshes, the combination of vertex welding algorithm with Temporal Branch-on-Need Octree (T-BON) algorithm is introduced (Wilhelms and Van Gelder, 1990). Although the triangular meshes are mapped by the index set referred to the vertices in the program, they originate from the triangular group found by the coordinate of the vertices instead of the index set. It's worth noting that the vertices of the triangle are independent of each other in the data structure, which means that the same vertex can be duplicated in two or more triangular meshes, so the vertices need to be accessed without redundancy to extract the correct triangle meshes, as shown in Fig. 10. The vertex welding algorithm can be summarized in the following steps:

If a2 = 0, the solution is given by:

x = (e2 c2 − d2 f2 − b2 g2) (2b2 c2)

(17)

y = (−e2 c2 + d2 f2 − g2 b2) (2b2 c2)

(18)

z = (−e2 c2 − d2 f2 + b2 g2) (2c2 d2)

(19)

If a2 ≠ 0, the solution is written as:

x=

(a2 f2 − b2 c2 )(c2 d2 − a2 g2 ) ⎤1 2 ⎫ 1 ⎧ −c2 ± ⎡ ⎥ ⎢ a2 ⎨ (a2 e2 − b2 d2 ) ⎦ ⎬ ⎣ ⎭ ⎩

(20)

y = −(d2 x + g2) (a2 x + c2)

(21)

z = −(b2 x + f2 ) (a2 x ) + c2

(22)

Substituting Eqs. (20), (21) and (22) into Eq. (2), the body saddle value can be determined. Considering the trilinear interpolate T within a particular volume cell, for a given threshold t, linear interpolation gives us the intersection points of the edges with the isosurface T(x, y, z) = t, which can be known as edge points. The threshold “t” is the value selected to segment the CT images. There are quite a few methods to determine the threshold “t” to segment the CT images, such as OSTU threshold algorithm (Cui and Jia, 2013), Minimum error threshold algorithm (Kittler and Illingworth, 1986), Isodata threshold algorithm (Ridler and Calvard, 1978), etc. In this paper, by comparing the segmented results of each threshold algorithm, the Isodata threshold algorithm is chosen as the threshold algorithm. The porosity of the segmented results by Isodata threshold algorithm is also the closest to the porosity of the actual porous rocks. If the edge points are connected to form triangles, the triangular approximation to the isosurface is generated. For a particular threshold value t, if two vertices which contain the same data (white or black) can be connected by a path in the cell or on its boundary, these two vertices can be considered as “connectable”. Furthermore, if the path consists of points on the edges of the volume cell, the two vertices are considered as “edge connectable”; if the path consists of points on the faces of the cell, the two vertices are considered as “face connectable”; if the path consists of points inside the body of the volume cell, the two vertices are considered as “body connectable”. The vertices in the volume cell which contain no ambiguity are edge-connectable, the vertices in the volume cell which contain ambiguity on the surface are face-connectable, and the vertices in the volume cell which contain ambiguity on the volume are body-connectable. Therefore, the ambiguity of the marching cubes can be eliminated based on the connection status of the vertices. The method can be best understood in two dimensions. Supposing that we want to construct the contour B(x,y) = t, in which t is the threshold, there are two possible topologies for a case with a saddle point at (xs, ys). As t is reduced, the contour becomes a horseshoe shaped curve that disconnects the black vertices. If t is reduced to the value less than B(xs, ys), the contour disconnects the white vertices. The saddle point B(xs, ys) is the transition point that decides which pair of vertices are disconnected. In other words, if t > B(xs, ys), the contours disconnect the vertices with valuable data which are more than t; If t < B(xs, ys), the contours disconnect the vertices with valuable data which are less than t. The topology for MC algorithm can be correctly determined by comparing the threshold with the value of the

(1) Merging vertices. The vertex welding algorithm begins by rounding vertices to a certain tolerance. (2) Welding vertices. The vertices are sorted and originate from a C ++ template library for CUDA and the duplicates are erased. Actually, after the vertex algorithm is finished in the program, the new triangle mesh group is introduced to a 3D vector for the further merger of coplanar triangles and removal of duplicate triangles by the T-BON Octree algorithm. For example, there are four volume cells arranged side by side, as shown in Fig. 11(a), and all of them follow the condition 8, so the isosurfaces in four volume cells can be fitted by eight triangle meshes, but the isosurface can also be fitted by another two larger triangle meshes, as shown in Fig. 11(b). Compared with the fitting mode in Fig. 11(a), the fitting mode in Fig. 11(b) seems to generate absolutely the same isosurface with fewer triangle meshes. Therefore, the T-BON Octree algorithm is naturally introduced to merge the coplanar triangles and remove the duplicate triangle so as to form the optimized isosurface. It's worth noting that the T-BON Octree algorithm is employed before the triangles are generated in the software. Therefore, the number of the triangles is directly reduced when the meshes are constructed. The parent nodes of Octree have eight child nodes, and mostly the child nodes can be further divided into eight Octants. If the maximum resolution is reached at one node, it is a leave of the tree. If the leaves of a node have exactly the same attributes, it is possible to reduce them to a node, this is called condensation, which is very important for the reduction of triangle meshes in this paper. If the empty volume cell is detected in the parent nodes by checking its value in this algorithm, all 376

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 11. The comparison between classic fitting mode and improved fitting mode.

of the volume cells below the parent nodes are empty volume cell, thus the corresponding child nodes do not need to be traversed so as to accelerate the calculation. Moreover, the branch-on-need strategy normally produces a tree shape quite different from the even sub division strategy in full Octree algorithm. The even-sub division strategy partitions the volume into equal parts. In rare cases, it may be necessary to visit more nodes in a T-BON Octree than in an even sub division Octree to process a particular section of the volume. Based on the modified Marching cubes algorithm, the CT image data of Chongqing sandstone is employed to generate the 3D surface model. High-resolution X-ray micro-computed tomography is applied to gather the CT image data from Chongqing sandstone. Some slices of the rock volumetric data are shown in Fig. 12: By employing the modified Marching cubes algorithm, the surface model is reconstructed, and the results are plotted in Fig. 13. 5. Finite element model based on improved Delaunay refinement algorithm

Fig. 13. The reconstructed surface model of Chongqing sandstone.

Based on the accurate definition in the mathematics of the relationship between the voxels in the CT image data and the vertex in the surface model by modified Marching Cubes algorithm, the surface model of Chongqing sandstone is exactly reconstructed. But the reconstructed surface model only satisfies the visualization of the interface in rocks. It is still unable to provide the satisfactory numerical model for investigation. Therefore, the improved constrained Delaunay triangulation algorithm is adopted to develop the surface model into the Finite element model to further investigate the relationship between the microstructures and the mechanical properties of rocks. The Delaunay triangulation algorithm is an incremental algorithm based on empty circumsphere criterion (Delaunay, 1934). It defines that the circumsphere of an element in 3D space does not contain any other points except for the four vertices of its own. By this algorithm, the point sets in 3D space are connected to each other to form the tetrahedral element sets. The input domain of the Delaunay refinement algorithm is the distinct vertices set (ITV) generated from the surface model reconstructed by Marching cubes algorithm. A modified constrained Delaunay triangulation algorithm is employed to compute the Delaunay tetrahedron set (DTS) to obtain the FE tetrahedral model (FTM). The surface model is defined as the convex domain that holds all points to be inserted. The modified constrained Delaunay triangulation algorithm is

described as follows and illustrated in Fig. 14: Step 1: The surface model is defined as the convex domain that holds all points to be inserted. Step 2: The constrained information of the boundary edge set and triangle set in the surface model is extracted. Step 3: One vertex from the initial tetrahedron vertices set (ITV) is inserted into the convex domain. The elements that violate the empty circumscribed sphere criterion are deleted to form an empty cavity, and the tetrahedral elements are generated by using the inserted vertex and the surface triangles of the cavity. Step 4: The generated tetrahedral elements are saved and Step 3 is repeated to insert all the distinct vertices into the convex domain to obtain the initial tetrahedron set (ITS). Step 5: The self-adaptive algorithm is employed to generate interior point set (IPS) and the IPS is inserted into ITS to form the Delaunay tetrahedron set (). Step 6: The boundary edges and the boundary triangles are recovered. Step 7: The tetrahedral elements outside the domain are eliminated. Step 8: The quality improving algorithm is employed to improve the quality of DTS. Step 9: The sliver elements in DTS are eliminated to obtain the final FE tetrahedral model.

Fig. 12. The processed CT images of Chongqing sandstone.

377

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 14. The modified constrained Delaunay triangulation algorithm.

Fig. 16. The pipe and the recovery of the edges. (a) The definition of a pipe. (b) The convex region formed by ACDE and CDEF. (c) The concave region formed by ACDE and CDEF.

inferior elements are generated due to the high irregularity of the microstructure in the natural rocks. Moreover, the integrity of boundary edges and triangles is not able to be guaranteed. Therefore, three types of modification procedure are employed to remove the bad quality elements and to recover the disappeared edges and triangle to ensure boundary integrity. The first procedure is to recover the disappeared edges and triangles. In the insertion of the vertices and the generation of the tetrahedral elements, some edges are disappeared because they intersect some tetrahedral elements. An element set called pipe is naturally generated. In the pipe, all elements are intersected by the disappeared edges, as shown in Fig. 16(a). The missing edge is AB. Pipe length is defined as the number of the elements in the pipe. When the pipe length equals 2, the disappeared edges are recovered by dividing two adjacent tetrahedra ACDE and CDEF in a convex polyhedron into three tetrahedra ACDF, ACEF and CDEF, as shown in Fig. 16(b). When the pipe length exceeds 2, there are two methods to decrease the pipe length. For the convex polyhedron constituted by two tetrahedral elements, the first two tetrahedral elements can be divided into three elements to decrease the pipe length by 1, as shown in Fig. 16(b). For the concave polyhedron constituted by two tetrahedral elements, the first two tetrahedral elements ACDE and CDEF can be divided into six tetrahedra AGCD,AGCE,AGDE, FGCD, FGDE and FGCE by the intersection point G between edge AB and triangle CDE, as shown in Fig. 16(c). Then the tetrahedral element GCEF is inserted into the pipe to replace ACDE and CDEF so as to decrease the pipe length by 1. Based on the three type of recovery procedure, every disappeared edge in the FE tetrahedral model can be recovered. When the triangles intersect the tetrahedral elements in the insertion of the vertices, the triangles are disappeared in the FE tetrahedral model. In this paper, there are four cases of the intersection between the disappeared triangle and tetrahedral elements. Considering the irregularity of the FE tetrahedral model generated by the CT images, the algorithm employed in this paper extends the method described in reference (Guan et al., 2006). The four cases are sorted by the number of the edges of the tetrahedral elements intersected by the disappeared triangles. When one edge of the tetrahedral elements ABCD is

Fig. 15. The self-adaptive algorithm.

In step 5, some self-adaptive points are added into the FE tetrahedral model to generate interior tetrahedral elements. At the same time, the added point should be avoided more than four points cosphericity. The self-adaptive algorithm can be summarized as follows and illustrated in Fig. 15: Step 1: The diametric ball of the triangle elements in the surface model is calculated to define the protection area set (PAS). Step 2: The grid-based method is employed to generate the initial interior point set (IIPS). Step 3: The random disarrangement algorithm is performed on IIPS to avoid more than four points co-sphericity. Step 4: IIPS is examined to eliminate the point located in PAS to obtain the interior point set (IPS). Step 5: IPS is inserted into the initial tetrahedron set ITS to obtain the Delaunay tetrahedron set (DTS). After the self-adaptive algorithm, a lot of bad quality, small and 378

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 17. Constraint triangle intersects edges of tetrahedron. (a) One edge intersection. (b) (c) Two edges intersection. (d) (e) (f) (g) Three edges intersection. (h) Four edges intersection.

in Fig. 18. But, if v encroaches any edge, v is rejected to be added and is delivered to step 3. Step 3. If one edge of the bad quality tetrahedron elements encroaches, then v is added at the midpoint of the edge. The third procedure is to eliminate the sliver elements. After the self-adaptive algorithm, the Delaunay tetrahedral model DTS almost satisfies the numerical simulation requirement. But there still exist some useless sliver elements which affect the convergence and efficiency of the numerical computation. In this paper, we use a quality parameter 3r/R to indicate the quality of the tetrahedral elements, where r and R are respectively the radius of the inscribed sphere and circumscribed sphere of the tetrahedral elements. By setting a threshold ξ = 0.01, the sliver elements are detected and eliminated if 3r/R < ξ. Based on the surface model, the FE tetrahedral model is reconstructed by using the modified Delaunay algorithm, as shown in Fig. 19.

intersected by the disappeared triangle CDE, CDE can be recovered by dividing ABCD into AECD and BCDE, as shown in Fig. 17(a). When two edges of the tetrahedral elements ABCD are intersected by the disappeared triangles BEF in Fig. 17(b) and (c), BEF can be recovered by dividing ABCD into ABEF, BCEF (BEDF) and BCDF (BCDE). When three edges of the tetrahedral elements ABCD are intersected by the disappeared triangles EFG in Fig. 17(d), (e), (f) and (g), EFG can be recovered by dividing ABCD into AEFG, CEFG, BCEG and BCDG, or by adding a point at the barycenter of EFGBCD to divide ABCD into AEFG, HBEF, HBCF, HBDE, HDEG, HCFG, HCDG, HEFG and HBCD. When four edges of the tetrahedral elements ABCD are intersected by part of the disappeared triangles EFGH in Fig. 17(h), EFGH can be recovered by dividing ABCD into ADEF, DEFG, DFGH, BEFH, BCFH and CFGH. The second procedure is to improve the bad quality elements. The radius-edge ratio ρ(π) = r/l is defined as the tetrahedron quality measure parameter, where r is the radius of the circumscribed ball of the tetrahedron element, l is the shortest edge of the tetrahedron element. When ρ (π ) ≥ 2 2 , the tetrahedron element is considered as bad quality one. An alternative vertex v is added into DTS to remove the bad quality tetrahedron element after the local mesh edges are reconnected to form the new tetrahedron elements. An edge or a triangle encroaches when at least one vertex is contained in its diametric ball. The vertex adding algorithm is described as follows: Step 1. If the radius-edge ratio of a tetrahedron element from DTS is greater than 2 2 , then v is added at the circumcenter of the tetrahedron element. While, if v encroaches any triangle, v is rejected to be added and is delivered to step 2. Step 2. If one triangle of the bad quality tetrahedron elements encroaches, then v is added at the circumcenter of the triangle, as shown

6. Numerical simulation and experimental verification Considering the huge number of small and inferior accompanying finite element meshes in the generation of finite element mesh (Areias and Rabczuk, 2017; Zhuang et al., 2016; Zhu et al., 2016), the combination of vertex welding algorithm with Temporal Branch-on-Need Octree (T-BON) algorithm is introduced to remove the duplication without any loss of detail. In this paper, Mohr-coulomb criterion is used to calculate the strength and plastic strain of Chongqing sandstone. The physical and mechanical properties of Chongqing sandstone obtained from experiment should not be directly assigned to the reconstructed model

Fig. 18. The vertex adding algorithm. (a) The vertex is added at the circumcenter of the bad quality triangle. (b) The local mesh edges are reconnected to form the new triangle.

Fig. 19. The reconstructed FE tetrahedral model.

379

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

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

Fig. 21. The displacement and maximum principle stress contour of the reconstructed model. (a) The displacement contour of the reconstructed model. (b) The maximum principle stress contour of the reconstructed model.

reconstructed model under different confining pressure are plotted in Fig. 22(a). The numerical uniaxial compressive strength (UCS) of the reconstructed model is plotted in Fig. 22(b). Compared with experimental data, the numerical uniaxial compressive strength and axial stress-strain curves are in good agreement with experimental data, which implies that the reconstructed model has the ability to safely predict mechanical properties of the porous sandstone. In order to prove that the numerical results are independent of the discretization, three reconstructed models with the same internal structure and different meshes are prepared for uniaxial compressive strength tests. The corresponding reconstructed models with different meshes are shown in Fig. 23. The total number of nodes and elements for different reconstructed models are listed in Table 4. The different load displacement curves for different meshes of the reconstructed model are plotted in Fig. 24. It is found from Fig. 24 that the load displacement curves for different meshes of the reconstructed model are almost the same. It indicates that the numerical results are independent of the discretization. Porosity has a significant influence on the mechanical properties of rocks. In order to investigate the influence of the pore structures on the mechanical properties of rocks, a series of reconstructed models with different porosity and fractal dimension are obtained by changing the threshold, as shown in Fig. 25. Then, a series of uniaxial compressive strength (UCS) tests are conducted on the reconstructed models. Firstly, by changing the threshold, ten different uniaxial compressive strength tests are conducted on the reconstructed models. The constraint conditions are plotted in Fig. 20(b). The corresponding results are listed in Table 5. The numerical results show that the rock strength decreases with increasing the porosity of rocks. When the porosity is very low, the rock can

because this leads to the second deduction in properties. Considering this, based on the experimental data of Chongqing sandstone, the numerical method is employed to calculate the proper parameters. Then the reconstructed model is simulated under different confining stress with the proper parameters. The calculated parameters used in the numerical simulation are listed in Tables 1 and 2. The reconstructed model is shown in Fig. 20(a). At the same time, considering the influence of microstructure evolution on the global stiffness, the flexibility matrix is calculated by the displacement of the elements. Considering the uneven surface, the reconstructed and reference model cannot be loaded and restricted conventionally. Therefore, rigid plates are applied 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 are fixed along all direction, as shown in Fig. 20(b). The confining pressure is also applied to the FE model by the rigid plates. The upper plate is loaded by displacement along the direction of Y-axis, as shown in Fig. 21(a). The maximum principal stress contour under confining pressure of 10 MPa is plotted in Fig. 21(b).The parameters obtained by the reconstructed model are listed in Table 3. Compared with Tables 1 and 2, the parameters obtained by the reconstructed model are close to the experimental data. The numerical axial stress-strain curves of the Table 3 The numerical parameters calculated by the reconstructed model. Density (kg/m3)

Young's modulus (GPa)

Poisson's ratio

UCS (MPa)

Cohesion (MPa)

Internal friction angle (°)

3171.49

24.015

0.2664

116.98

34.7952

27.4

380

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 22. Comparison of 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 the peak strength of Chongqing sandstone sample.

be considered approximately as intact. There are quite a few studies on the relationship between the strength and the porosity of porous materials. Several general models have been proposed for porous materials. Ryshkewitch (1953) proposed the relation between the strength and the porosity of Al2O3 and ZrO2 as follows:

σ = σ0 e−k1 p

Table 4 The number of nodes and elements for the reconstructed models.

(23)

Model number

Total number of nodes

Total number of elements

Model 1 Model 2 Model 3

89,067 170,406 392,519

70,345 137,183 323,445

where σ is the strength of porous materials, σ0 is the strength of porefree materials, k1 is the empirical constant, p is the porosity. Schiller (1971) obtained the relation between the strength and the porosity of sulfate plasters as follows:

p σ = n ln ⎛⎜ 0 ⎞⎟ ⎝ p⎠

(24)

where n is the empirical constant. Balshin (1949) suggested the relation between the strength and the porosity of the metal ceramics as follows:

σ = σ0 (1 − p)b1

(25)

where b1 is the empirical constant. Hasselman (2010) proposed the relation between the strength and porosity of different refractory materials as follows:

σ = σ0 − c0 p

(26) Fig. 24. The axial stress-strain curves for different reconstructed models.

where c0 is the empirical constant. Zheng et al. (1992) obtained the relation between the strength and the porosity based on the brittle fracture theory proposed by Griffith (1920) as follows:

1 2

⎡ p − p⎞ 2 3)⎤ KIc = KIco ⎢ ⎜⎛ c ⎟ ·(1 − p ⎥ pc ⎠ ⎝ ⎣ ⎦

m 2

p − p⎞ σ = α·⎜⎛ c ⎟ ⎝ pc ⎠

·KIc ,

KIc =

2γE

(28)

where KIco is the fracture toughness of pore-free materials. It is assumed that σ0 = α · KIco is the strength of pore-free materials, Eq. (28) can be rewritten as:

(27)

where α is a coefficient concerning stress state, pc is the percolation porosity at failure threshold, m is the ratio of calculated average distance to the nearest pore, γ is the fracture surface energy. The porosity dependence of the fracture toughness (Wagh et al., 1993) is written as

1.85

⎡ p − p⎞ σ = σ0 ⎢ ⎜⎛ c ⎟ pc ⎠ ⎣⎝

1 2

⎤ ·(1 − p2 3 ) ⎥ ⎦

Fig. 23. The reconstructed models with different meshes. (a) Model 1. (b) Model 2. (c) Model 3.

381

(29)

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 25. 2D CT image data with different segment threshold. (a) Image data with porosity of 1.23% and fractal dimension of 0.9743. (b) Image data with porosity of 24.29% and fractal dimension of 1.5572. (c) Image data with porosity of 49.52% and fractal dimension of 1.7677.

The porosity-UCS curves obtained by the reconstructed model and the theoretical curves of the previous models are plotted in Fig. 26. The correlation coefficient between the numerical results and the theoretical curves is listed in Table 6. It is observed from Fig. 26 that the numerical results are in good agreement with the theoretical curves, and it's found from the correlation coefficient that Zheng's model is the best to illustrate the relationship between the porosity and the strength of the reconstructed model. In addition, Zheng's model was also proved to be the best to illustrate the relation between the porosity and the strength of porous materials by experiments (Chen et al., 2013). The Young's modulus of porous materials, such as porous sandstone, is also deeply related to the porosity. Numerous theoretical models were also proposed to predict the relationship between Young's modulus and porosity. Particularly, reliable relations for predicting the relative Young's modulus of porous materials include the Maxwell-type model (Torquato, 2002), the self-consistent model (Budiansky and O′Connell, 1976), the differential scheme model (Hashin, 1988) and the exponential model (Pabst and Gregorová, 2004), etc. The Maxwell-type model is (Torquato, 2002)

Er =

1−p , 1+p

Er =

E E0

Fig. 26. Numerical and theoretical porosity vs UCS.

Table 6 The correlation coefficient of the porosity-UCS theoretical model.

(30)

where E0 is Young's modulus of pore-free materials, E is Young's modulus of porous materials, Er is relative Young's modulus. The self-consistent model is expressed as (Budiansky and O′Connell, 1976)

Er = 1 − 2p

(31)

(32)

The exponential model is (Pabst and Gregorová, 2004)

−2p ⎞ ⎟ Er = exp ⎛⎜ 1 ⎝ − p⎠

Correlation coefficient

Zheng's model Ryshkewitch's model Schiller's model Balshin's model Hasselman's model

0.995 0.947 0.934 0.957 0.951

good agreement with the theoretical results. Moreover, it is found from the correlation coefficient that the numerical results are closest to the exponential model proposed by Pabst et al., 2013). In addition, the exponential model was also verified to be the best to illustrate the relationship between Young's modulus and the porosity by experiments (Pabst et al., 2013). With the development of rock mechanics, it is proved that the pore size distribution has a significant influence on the strength of porous rocks. The pore size distribution is taken into consideration in the strength model in different forms. However, these models are not still satisfied to illustrate the relationship between pore structure and strength, the reason is that the parameters used to represent the pore size distribution are not accurate enough (Tian and Bu, 2014). Fractal dimension is a very effective parameter to characterize the pore size

The differential scheme model is written as (Hashin, 1988)

Er = (1 − p)2

Theoretical model

(33)

The numerical results and the theoretical curves are depicted in Fig. 27, and the corresponding correlation coefficient is given in Table 7. In general, Young's modulus increases with decreasing the porosity. The numerical relation between Young's modulus and the porosity is in Table 5 The numerical results for different reconstructed model. Porosity (%) UCS (MPa) Young's modulus (GPa)

0.59 199.46 37.18

1.23 190.22 36.70

2.58 187.32 35.68

4.89 175.21 33.94

8.76 150.86 30.05

382

14.93 139.46 27.48

18.33 110.98 23.02

24.29 104.80 20.80

36.03 50.76 11.20

49.52 17.47 5.29

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

Fig. 27. Numerical and theoretical Young's modulus vs UCS.

Fig. 28. The numerical and theoretical fractal dimension vs UCS.

Table 7 The correlation coefficient of Young's modulus vs UCS.

Table 9 The correlation coefficient of fractal dimension vs UCS.

Theoretical model

Correlation coefficient

Theoretical model

Correlation coefficient

Exponential model Maxwell-type model Self-consistent model Differential scheme model

0.992 0.863 0.921 0.918

Extended classical Freundlich model Linear fitting model Rational fitting model Exponential fitting model

0.989 0.863 0.794 0.802

distribution and the complexity of the pore structure in porous rocks. By changing the threshold, a series of reconstructed models with different fractal dimension are obtained. Several regular fitting models are employed tentatively to investigate the relationship between the uniaxial compressive strength and fractal dimension, such as Linear fitting model, Rational fitting model, Exponential fitting model and Extended classical Freundlich model, etc. Linear fitting model is

7. Discussions and conclusions In this paper, the spatial configuration and geometrical information of the microstructures are firstly gathered from Chongqing sandstone by X-ray micro-CT observations, and experimental tests on Chongqing sandstone are conducted to obtain the mechanical parameters. Then, the volumetric CT image data are preprocessed by adaptive total variation algorithm to eliminate the noise, and the proper threshold algorithm is employed to distinguish rock matrix from the pore. The classic Marching cubes algorithm is improved by introducing the vertices welding algorithm and Octree algorithm. Based on the processed CT image data, the surface geometrical model is established by using the modified Marching cubes algorithm. Next, a modified constrained Delaunay triangulation algorithm is employed to generate the FE tetrahedron model based on the input domain from the surface geometrical model. In the modified constrained Delaunay triangulation algorithm, a self-adaptive algorithm is introduced to generate interior tetrahedral elements, the disappeared elements are recovered, and the bad quality elements are eliminated. Finally, the 3D tetrahedral model is transformed into FE model for numerical analysis. The relation among the porosity, fractal dimension of pore structure and mechanical properties of porous sandstone is investigated. The numerical results are in good agreement with the experimental data and are illustrated well by the theoretical models based on experiments. It is indicated that the proposed reconstruction method has a good ability to predict mechanical properties of porous rocks. However, this method is deeply limited by the precision of the X-ray CT scanner, and can only partially reflect the characteristics of rocks, especially for the various geological condition in the engineering zone. In order to completely and objectively evaluate the various geological condition, a considerable number of X-ray micro-CT observations are still needed, which implies that the resource consumed by this

(34)

σ = a + bDs where Ds is fractal dimension, a and b are the fitting parameters. Rational fitting model is expressed as

(35)

σ = c3 + b (Ds + a) where a, b and c3 are the fitting parameters. Exponential fitting model is written as

(36)

σ = a × b Ds Extended classical Freundlich model is given by

σ = a + b × (Ds )c3

(37)

The numerical results are given in Table 8. The numerical results and the theoretical fitting curves are plotted in Fig. 28. Fig. 28 shows that UCS of porous rocks increases with increasing the fractal dimension. It is observed from Fig. 28 that the numerical results are in good agreement with the theoretical fitting curves. Moreover, it is found from Fig. 28 and Table 9 that the numerical results are closest to the Extended classical Freundlich model. In addition, when the constant a in Eq. (37) is equal to zero, Eq. (37) is rewritten as follows:

σ = b × Ds a

(38)

which is consistent with the prediction model from experiments (Jin et al., 2017). Table 8 Fractal dimension vs UCS by the reconstructed model. Porosity (100%) Fractal dimension UCS (MPa)

0.59 0.9113 199.46

1.23 0.9743 190.22

2.58 1.0644 187.32

4.89 1.1746 175.21

8.76 1.3011 150.86

383

14.93 1.4215 139.46

18.33 1.4714 110.98

24.29 1.5572 104.80

36.03 1.6736 50.76

49.52 1.7677 17.47

Engineering Geology 228 (2017) 371–384

X.-P. Zhou, N. Xiao

computer tomograms. Biomed. Tech. 24 (1), 1–21. Jia, R.Q., Zhao, H., 2010. A fast algorithm for the total variation model of image denoising. Adv. Comput. Math. 33 (2), 231–241. Jin, S., Zhang, J., Han, S., 2017. Fractal analysis of relation between strength and pore structure of hardened mortar. Constr. Build. Mater. 135, 1–7. Kittler, J., Illingworth, J., 1986. Minimum Error Thresholding. Elsevier Science Inc. Lanzón, M., Cnudde, V., Kock, T.D., Dewanckele, J., Piñero, A., 2014. X-ray tomography and chemical–physical study of a calcarenite extracted from a roman quarry in Cartagena (Spain). Eng. Geol. 171 (8), 21–30. Lewiner, T., Lopes, H., Vieira, A.W., Tavares, G., 2003. Efficient implementation of marching cubes cases with topological guarantees. J. Graph. Tools 8, 1–15. Lorensen, W.E., Cline, H.E., 1987. Marching cubes: a high resolution 3D surface construction algorithm[J]. Acm Siggraph Comput. Graph. 21 (4), 163–169. Lorensen, W.E., 1987. Marching cubes: a high resolution 3D surface construction algorithm. In: ACM SIGGRAPH Computer Graphics. 21 (4). pp. 163–169. Nielson, G., 2003. On marching cubes. IEEE Trans. Vis. Comput. Graph. 9 (3), 283–297. Pabst, W., Gregorová, E., Černý, M., 2013. Isothermal and adiabatic Young's moduli of alumina and zirconia ceramics at elevated temperatures[J]. J. Eur. Ceram. Soc. 33 (15–16), 3085–3093. Pabst, W., Gregorová, E., 2004. Mooney-type relation for the porosity dependence of the effective tensile modulus of ceramics[J]. J. Mater. Sci. 39 (9), 3213–3215. Ridler, T.W., Calvard, S., 1978. Picture thresholding using an iterative selection method. In: IEEE Transactions on Systems, Man and Cybernetics. vol. 8. pp. 630–632. Ryshkewitch, E., 1953. Compression strength of porous sintered alumina and zirconia. J. Am. Ceram. Soc. 36 (2), 65–68. Sarkar, G., Siddiqua, S., 2016. Effect of fluid chemistry on the microstructure of light backfill: an X-ray CT investigation. Eng. Geol. 202, 153–162. Schiller, K.K., 1971. Strength of porous materials. Cem. Concr. Res. 1 (4), 419–422. Schlueter, E.M., Zimmerman, R.W., 1997. The fractal dimension of pores in sedimentary rocks and its influence on permeability. Eng. Geol. 48 (3), 199–215. Tian, Z., Bu, J., 2014. Mathematical model relating uniaxial compressive behavior of manufactured sand mortar to MIP-derived pore structure parameters. Sci. World J. 2014, 736230. Torquato, S., 2002. Random Heterogeneous Materials. Springer, New York. Uday, K.V., Padmakumar, G.P., Singh, D.N., 2013. Some studies on morphology of the coarse-grained soils. Eng. Geol. 152 (1), 48–55. Vignoles, G.L., Donias, M., Mulat, C., Germain, C., Delesse, J.F., 2011. Simplified marching cubes: an efficient discretization scheme for simulations of deposition/ ablation in complex media. Comput. Mater. Sci. 50 (3), 893–902. Wagh, A.S., Singh, J.P., Poeppel, R.B., 1993. Dependence of ceramic fracture properties on porosity. J. Mater. Sci. 28 (13), 3589–3593. Wilhelms, J., Van Gelder, A., 1990. Octrees for faster isosurface generation. ACM Trans. Graph. 24 (5), 57–62. Yu, W., Jie, P., Wang, L., Wang, J., Zheng, J., Song, Y.F., et al., 2016. Characterization of typical 3D pore networks of Jiulaodong formation shale using nano-transmission Xray microscopy. Fuel 170, 84–91. Yun, T.S., Jeong, Y.J., Kim, K.Y., Min, K.B., 2013. Evaluation of rock anisotropy using 3d X-ray computed tomography. Eng. Geol. 163 (12), 11–19. Zheng, M., Zheng, X., Luo, Z.J., 1992. Fracture strength of brittle porous materials. Int. J. Fract. 58 (3), R51–R55. Zhou, C., Shu, R., Kankanhalli, M.S., 1994. Handling small features in isosurface generation using marching cubes. Comput. Graph. 18 (6), 845–848. Zhou, B., Wang, J., Zhao, B., 2015. Micromorphology characterization and reconstruction of sand particles using micro X-ray tomography and spherical harmonics. Eng. Geol. 184 (14), 126–137. Zhou, S., Yan, G., Xue, H., Guo, W., Li, X., 2016. 2D and 3D nanopore characterization of gas shale in Longmaxi formation based on FIB-SEM. Mar. Pet. Geol. 73, 174–180. Zhu, H., Wu, W., Chen, J., Ma, G., Liu, X., Zhuang, X., 2016. Integration of three dimensional discontinuous deformation analysis (DDA) with binocular photogrammetry for stability analysis of tunnels in blocky rockmass. In: Tunnelling & Underground Space Technology Incorporating Trenchless Technology Research. 51. pp. 30–40. Zhuang, X., Wang, Q., Zhu, H., 2016. Multiscale modelling of hydro-mechanical couplings in quasi-brittle materials. Int. J. Fract. 1–27.

reconstructed method is not much less than experiments. The next step is to introduce the randomness into this method so as to predict more types of porous rocks in various geological condition with less X-ray micro-CT observations. 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). Reference Almeida, M.S.C., Almeida, L.B., 2010. Blind and semi-blind deblurring of natural images. IEEE Trans. Image Process. 19 (1), 36–52. Alvarez, L., Mazorra, L., 1994. Signal and image restoration using shock filters and anisotropic diffusion. SIAM J. Numer. Anal. 31 (2), 590–605. Areias, P., Rabczuk, T., 2017. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elem. Anal. Des. 132, 27–41. Balshin, M.Y., 1949. Relation of mechanical properties of powder metals and their porosity and the ultimate properties of porous-metal ceramic materials. Dokl Askd SSSR 67 (5), 831–834. Bowyer, A., 1981. Computing Dirichlet tessellations. Comput. J. 24, 162–166. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of cracked solid. Int. J. Solids Struct. 12 (2), 81–97. Chen, X.D., Wu, S., Zhou, J., 2013. Influence of porosity on compressive and tensile strength of cement mortar[J]. Constr. Build. Mater. 40 (3), 869–874. Christe, P., Turberg, P., Labiouse, V., Meuli, R., Parriaux, A., 2011. An X-ray computed tomography-based index to characterize the quality of cataclastic carbonate rock samples. Eng. Geol. 117 (3–4), 180–188. Cline, H., Lorensen, W., 1988. Two algorithms for the three-dimensional reconstruction of tomograms. Med. Phys. 15 (3). Cueto, N., Benavente, D., Martínez-Martínez, J., García-Del-Cura, M.A., 2009. Rock fabric, pore geometry and mineralogy effects on water transport in fractured dolostones. Eng. Geol. 107 (1), 1–15. Cui, Z.D., Jia, Y.J., 2013. Analysis of electron microscope images of soil pore structure for the study of land subsidence in centrifuge model tests of high-rise building groups. Eng. Geol. 164 (12), 107–116. Delaunay, Boris, 1934. Sur la sphère vide. Bull. Acad. Sci. URSS 6, 793–800. Gent, A.N., Thomas, A.G., 1959. The deformation of foamed elastic materials. J. Appl. Polym. Sci. 1, 107–113. Gent, A.N., Thomas, A.G., 1963. Mechanics of foamed elastic materials. Rubber Chem. Technol. 36, 597–610. Gibson, L.J., Ashby, M.F., 1997. Cellular Solids: Structures and Properties. Cambridge University Press, Cambridge. Goldstein, T., Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM J. Imag. Sci. 2 (2), 323–343. Grasmair, M., 2009. Locally adaptive total variation regularization. In: The 2nd International Conference on Scale Space and Variational Methods in Computer Vision, pp. 331–342. Griffith, A.A., 1920. The theory of rupture and flow in solids. Trans. R. Soc. Lond. A 221, 63. Guan, Z., Song, C., Gu, Y., 2006. The boundary recovery and sliver elimination algorithms of three-dimensional constrained Delaunay triangulation. Int. J. Numer. Methods Eng. 68 (2), 192–209. Hashin, Z., 1988. The differential scheme and its application to cracked materials. J. Mech. Phys. Solids 36 (6), 719–734. Hasselman, D.P.H., 2010. Griffith flaws and the effect of porosity on tensile strength of brittle ceramics. J. Am. Ceram. Soc. 56 (9), 491. Herman, G.T., Liu, H.K., 2010. Three dimensional display of human organs from

384