The understanding of rock mechanical and transport properties on micro/mesoscopic is instructive to many engineering applications, such as coal safe mining (McKee et al. 1987; Jasinge et al. 2011), coal-bed methane recovery (Konecny and Kozusnikova 2011; Chen et al. 2013), shale gas exploitation (Spencer 1989; Kwon et al. 2001), CO₂ geological sequestration (Klinkenberg 1941; Liu and Rutqvist 2010), contamination migration control (Sadhukhan et al. 2012), etc. Natural rock, such as coal and hydrocarbon reservoir, inherently contains a large number of different scale pores and fractures, which provide the storage place and transport path for coal-bed methane, crude oil and natural gas. In recent decades, coalbed methane (CBM) has become an important energy resource around the world. The CBM is usually absorbed in the porous structure of deep coal seam. Without properly control of the methane desorption along with mining may cause severe disasters in deep coal mining, such as gas outburst. Conventionally, CBM is directly discharged into the atmosphere along with mining through ventilation which contributes a lot to global warming. In contrast, the CBM turns out to be a valuable resource if we could exploit it safely. The understanding of coalbed methane absorption–desorption mechanism in pore- scale could provide an opportunity for exploiting this valuable gas simultaneously with the coal mining (Sampath 1982). As the application of hydraulic fracturing techniques, the unconventional hydrocarbon resources such as shale gas and tight oil can be extracted with commercial rates. The study of crack propagation mechanism in pore-scale can shed light on the influence of porous structure, pre-existing cracks and mineral constitution to the crack propagation (Moore and Lockner 1995). On the contrary, the long-term steady state of fluid stored in the porous rock is the most concerned issue in CO₂ geological sequestration and contamination migration control (Liu et al. 2011). Fundamentally, the research in pore-scale of rock can elucidate the evolution mechanism of porous structure under pressure and its interaction with fluid which may make clear the macroscopic phenomenal of engineering application (Morrow 1984; Ju et al. 2008; 2013).

Obtaining digital porous structure of rock is usually through experimental methods and numerical reconstruction methods. The experimental methods mainly include X-ray computed tomography (CT) (Ketcham and Carlson 2001; Cnudde and Boone 2013), scanning electron microscope (SEM) (Jones et al. 1981), laser confocal microscopy (Fredrich et al. 1993). Among them, the X-ray CT method could provide 3D rock porous structure without destroy the rock sample. In addition, with the application of central synchrotron facility, the micro CT (μ-CT) (Arns et al. 2007) could be used to access the high resolution (several microns) images of rock porous structure. However, obtaining porous structure through the CT scanning test is laborious and expensive. First of all, the rock cores drilling could cost millions of dollars as human activity extending to geological complex regions. On the other hand, high resolution X-ray CT equipment is too expensive to be used widely. The CT scanning conducted on limited amount of core samples is incapable of uncovering the general relationship between the micro/mesoscopic porous structure and the macroscopic rock properties. Repeated drilling at different sites and more experimental analyses are required in order to get accurate rock property specifications. This dilemma drives researchers to develop a more economical way, i.e. numerical reconstruction method, to obtain the porous structure of rocks.

The reconstruction methods can be mainly divided into statistical information based method (Yeong and Torquato 1998a, b) and rock diagenetic process based method (Øren and Bakke 2002). The first one is achieved by matching the statistical information of the reconstruction model with the real rock model. This method could be applied to reconstruct any type of porous rock by matching the different type of statistical information but it is sometimes hard to reproduce the long-range connectivity of porous structure. The second one reconstructs the rock model by simulating the diagenetic process of sedimentary rock, like compaction and cementation. It is capable to reproduce the long-range connectivity of pores but the application is limited to certain type of rock with clearly distinguished grain size and specific diagenetic process.

The method we employed in this paper is a statistical information based reconstruction method which is based on the improved simulated annealing algorithm (SAA). A fractal descriptor which characterizes the complexity of pore morphology is introduced into the reconstruction procedure and act as one of the statistical information control function. For more information about the reconstruction procedure, please refer to the paper 3D Numerical Reconstruction of Well Connected Porous Structure of Rock Using Fractal Algorithms (Ju et al. 2014).

The research of rock mechanical behavior can be divided into the following three aspects: conventional experimental test on rock samples, using simplified models to explain the mechanical behaviors of rocks and directly calculate the mechanical behaviors from the porous rock structure. The conventional mechanical behavior test considers the overall rock mechanical behavior of rock samples and relates the test results to porosity, mineral content, the existence of cracks, fractures etc. Sone and Zoback (Sone and Zoback, 2013) systematically studied the mechanical behavior of shale gas reservoir rock and conclude that the understanding of elastic properties of reservoir rocks is crucial for gas resources exploitation. (Xie et al. 1997) used fractal geometry to describe the roughness of rock joints and employed photoelastic method to study the influence of fracture roughness to the rock mechanical behavior. With the experimental results available, an empirical relationship may then be proposed statistically to relate the rock mechanical properties with apparent rock physical properties, such as porosity, fracture density. Besides the conventional method, a lot of simplified models have been proposed to study the rock mechanical behavior and some models have incorporated the microstructure information. The representative simplified models include: mineral spatial average scheme model which mainly concerns the influence of mineral content and could give prediction of the upper and lower elastic modulus of rock. (O’Connell and Budiansky 1974; Walsh and Brace 1984; Han et al. 1986; Hudson et al. 1996). Liu et al. (Liu et al. 2009) considered portion of pores and fractures as soft part and other parts of rock as hard part and derive a constitutive relations between stress and rock mechanical properties; Ju et al. (Ju et al. 2013) build simplified physical and numerical models of porous sandstone which only contain sphere pores and investigated the influence of porosity on the static and dynamic rock mechanical behavior. In recent years, as the development of computer calculation capacity and accessibility of micro/mesoscopic rock porous structure, researchers started to study the rock mechanical properties directly based on the rock porous structure which called digital rock physics. Lots of researchers (Garboczi 1998; Arns et al. 2002; Makarynska et al. 2008) used finite element method to study the rock mechanical and physical properties based on the digital CT images directly. Saenger et al. (Saenger et al. 2011) developed a finite differential method to study the dynamic effective properties of porous rock. However, to the best of our knowledge, most of the digital physics numerical methods only gave estimates of elastic modulus and did not take into account the non-linear elastic rock behavior. We employed in this paper a non-linear elastic to plastic finite element method, called SciFEM, to study the influence of porous structure to rock mechanical behaviors.

The outline of the remainder of the paper is as follows: In Sect. 1, we present the reconstructed porous structure using the improved SAA. We also compare the reconstructed model with the real rock model in terms of statistical information and geometric morphology. In Sect. 2, the details how to convert the digital rock images into 3D meshed finite element models are described. The mechanical response of reconstructed model and real rock model will be compared and discussed in Sect. 3. We make concluding remarks in Sect. 4.

2.1 CT image processing

The low-permeability sandstone CT image was obtained by X-ray CT machine (ACTIS300-320/225) at State Key Laboratory of Coal Resources and Safe Mining. The rock sample was 25 mm in diameter and 10 mm in height. The original CT image was 1,024 × 1,024 × 400 voxels with resolution of 25 um. Figure 1a shows one slide of the original CT images. As we can clearly distinguish, there were four components: the black which represent pores, the dark gray maybe feldspar, the bright gray maybe Quartz and bright white maybe metal mineral.

Original image and binary image

In order to highlight the influence of porous structure, we convert the raw CT images to binary images which only contain pore phase and solid phase. In this tentative study, we consider the three different mineral phases as one rock matrix phase. One layer of CT image and its binary result are show in Fig. 1. It should be noted that some noisy parts, which defined as isolated and less than 125 voxels, were eliminated to ensure the precise of reconstruction and the following mechanical property calculation. For the convenience of the reconstruction, the center cubic box (512 × 512 × 400 voxels) of the images was extracted, binarized and set as the reference model in the reconstruction procedure.

2.2 Rock model reconstruction method

The framework of the reconstruction method was based on the SAA which was first introduced by (Yeong and Torquato 1998a, b). At the beginning, a randomly distributed black-and-white voxels structure with the same black and white voxels number of that in the reference model. The statistical information was extracted from both the reference model and the initial reconstruction model and quantified in terms of control function values. Then, continually tries of position interchange of a pair of black and white voxels in the reconstructed model were performed to make sure its morphology evolved similar with the reference model. When a trial position interchanging was performed, the control function values of the new model were calculated. Let f₀(r _i ) the control function values of the reference model, and let f(r _i ) and f′(r _i ) the control function values of the reconstruction model before and after the trial position interchange, respectively. The difference between the reconstruction models with the reference model, which called “Energy”, can be calculated by the following formulations respectively:

where E represented the “Energy” difference between the reconstruction model and reference model before the interchanging; similarly, \( E^{\prime} \) represented the “Energy” difference between the reconstruction model and reference model after the interchanging. The trial position interchanging was accepted with the probability p(ΔE) via the Metropolis criteria:

where T was decreasing slowly with the evolution of the reconstructed model. The initial value of T should be set at a relative high value that ensured the beginning acceptance rate of interchanging is about 75 %.

Besides the traditionally employed two point probability function and linear-path control function, we developed a fractal dimension description control function to better characterize the complex pore morphology. Besides, we added a pre-conditioning method in the later stage of the reconstruction procedure to increase the reconstruction efficiency.

Based on the same statistical information extracted in the reference rock model but different random sequence of the interchanging of voxels, two models were reconstructed. As shown in Fig. 2, the reconstructed models have different porous structure with the real rock model. But they have almost the same statistical information of the porous structure, as shown in Fig. 3.

Visualization of reference and reconstruction models

Comparison of function values. Notesa1, a2, a3 show the comparison of two-point probability control function value of the two reconstructed model with the real rock model in X, Y, Z direction; b1, b2, b3 show the comparison of linear-path control function value in x, y, z direction; c1, c2, c3 show the comparison of fractal descriptor control function value in xy, yz, zx direction

To further testify the reconstruction accuracy, we compared the pore number, average pore volume, average pore area and average pore specific surface area of the reference model and the reconstructed model. These parameters are listed in Table 1. As we can see in Table 1, the pore number of the reconstructed models is less than the reference model. This may be caused by the fact that there are some small isolated pores in the reference model which less appear in the reconstructed model. The pore specific surface areas are similar between the reference model and the reconstructed model, reflecting that the pore geometry of the models are similar.

We also employed the fractal dimension concept to analyze and compare the pore geometry between the reconstructed model and the reference model. As we known, the area and volume of regular shape pores have the relationship as follows:

Where A _i and V _i are the area and volume of each pore. For irregular shape pores, i.e. pores in rock, this relationship can be revised as:

where D is the fractal dimension which is bigger or equal to 2 and less than 3. The results of applying fractal dimension analysis to the pore space are shown in Fig. 4.

Fractal dimension analysis

Until to now, the morphology of porous structures has been established. In the next section, we will use our finite element method to calculate the non-linear elastic to plastic mechanical behaviors of the reference and reconstructed porous structures.

Applying the MIMICS^® software (http://biomedical.materialise.com/mimics), we converted the geometrical models into finite element models. First the images were input into the MIMICS^® software and a 3D entity was constructed using the built-in algorithm “Create Masks” and “Calculate 3D Models”. During this process, the pore and rock matrix position information was recorded. Secondly, the entity was uniformly meshed into tetrahedron elements using “Remesh” and “Create volume mesh”, see Fig. 5a. Each triangle in the meshes must obey the following rules to ensure the quality of the meshes: (1) the maximum side length must be shorter than the represented length of 10 voxels, and (2) the minimum height/bottom side length ratio must be larger than 0.4. Finally, each tetrahedron element is attributed with either rock matrix or porous material properties, according to its position information recorded in the 3D entity. Through this method, the pore phase and matrix phase was distinguished by different material properties in the FEM calculation, see Fig. 5b. This method has the advantage in efficiently meshing the entity since it does not need to consider the complex interface of the pore phase and matrix phase in the meshing stage. The morphology information of porous structure was preserved by limiting the maximum size of the element. However, the meshing would inevitably lose some detailed information of the interface. The meshed model information of the reference model and the two reconstructed models were subsequently converted to the format that can be used in the SciFEA, i.e. the self-developed nonlinear software (Ju et al. 2014).

Meshes inside of the entity and the finite element model with different material properties applied to different phases

Compared to most previous studies which only concerned the elastic stage of rock mechanical behaviors, the elastic to plastic transformation of porous structure was considered in the SciFEA software. The Drucker–Prager criterion was set as the elastic to plastic determination function. After the stress state exceeded the criterion, the stress of an element would be set to decrease rather than increase along with strain increasing, and finally down to the residual strength. In this way, the stress change corresponding to the strain change was related to the stress history of that element. In the calculation, a small increase of stain was applied to the model and the elastic or plastic index was calculated to determine whether elastic strain or plastic strain calculation function should be used. Then the calculation was iterated until the new balance state of the model.

The uniaxial compression experiment was simulated using the SciFEM software. The displacement of nodes at the bottom of the model was set at zero in all directions. The displacement loading in Z direction was applied to nodes at the top of the model. The material properties of the two types of elements were given in Table 2. It should be noted that the pore phase was considered as bonding component rather than voids.

The calculation results are presented in Fig. 6. The overall mechanical responses of the reconstructed models are similar with that of the real rock model. From Figs. 6 b1–b3 we could find that the stress value in matrix phase is clearly higher than that of pore phase. It is not hard to understand that the matrix phase which has higher elasticity modulus undertakes most of the loading. As show in Fig. 6 c1–c3, the deformations of the models are quite uniform in elastic stage. However, with the increasing of the loading, the stress in some elements exceeded its elastic limit and plastic deformation happened. As show in Figs. 6 d1–d3 the stress distributions are non-uniform with higher values in the middle. As display in Figs. 6 e1–e3 the higher deformation occurred in the corner of the models. Anyhow, the reconstructed models showed consistency with the real rock model in deformation and stress distribution, before and after elastic stress limit. This illustrated that the reconstructed models have the similar mechanical properties with the real rock model and could be used in the prediction of the low-permeability sandstone mechanical properties. The stress–strain curves of the reconstructed models and the real rock model were plotted in Fig. 7. The overall mechanical properties, i.e. elastic modulus, Poisson’s ratio and strength, for the reconstruction model and the reference model are listed in Table 3.

Comparison of the mechanical responses between the reconstructed models and real rock model. Notesa Finite element model images, b contour of stresses in Z direction before failure, c displacement in Z direction before failure, d contour of stresses in Z direction after failure, e displacement in Z direction after failure. The subscript “1”, “2”, “3” represent the reference rock model, the first reconstructed model and the second reconstructed model, respectively

Stress-strain curves of the reference rock model and the reconstructed models under uniaxial compression test

In this paper, we develop a numerical reconstruction method to simulate a large-size porous structure of low-permeability sandstones. The comparison indicates that the reconstructed models have the similar statistical information, geometrical morphology with those of the reference model.

The geometric porous structures models are converted into finite element models for calculating their mechanical properties. The numerical results show a great similarity between the reconstructed models and the reference model.

The similarity in geometry and mechanical properties between the reconstructed models and the real rock model illustrate the feasibility of using a numerical reconstructed model to predict the mechanical properties of the low-permeability sandstone.