Geofluids
Volume 2017 (2017), Article ID 7510527, 12 pages
https://doi.org/10.1155/2017/7510527

Research Article

A Pore-Scale Simulation on Thermal-Hydromechanical Coupling Mechanism of Rock

Rui Song,¹ Mengmeng Cui,² Jianjun Liu,1,3 P. G. Ranjith,⁴ and Yun Lei1,5

¹School of Geoscience and Technology, Southwest Petroleum University, Chengdu 610500, China
²School of Petroleum and Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China
³State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China
⁴Deep Earth Energy Research Laboratory, Department of Civil Engineering, Monash University, Melbourne, VIC 3800, Australia
⁵Shenyang Research Institute, China Coal Technology & Engineering Group Corp., Shenyang 110016, China

Correspondence should be addressed to Rui Song; songrui0506@126.com and Jianjun Liu; liujj0906@sina.com

Received 24 February 2017; Accepted 28 March 2017; Published 18 April 2017

Academic Editor: Yi Wang

Copyright © 2017 Rui Song et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Thermal-hydromechanical (THM) coupling process is a key issue in geotechnical engineering emphasized by many scholars. Most existing studies are conducted at macroscale or mesoscale. This paper presents a pore-scale THM coupling study of the immiscible two-phase flow in the perfect-plastic rock. Assembled rock matrix and pore space models are reconstructed using micro-CT image. The rock deformation and fluid flow are simulated using ANSYS and CFX software, respectively, in which process the coupled physical parameters will be exchanged by ANSYS multiphysics platform at the end of each iteration. Effects of stress and temperature on the rock porosity, permeability, microstructure, and the displacing mechanism of water flooding process are analyzed and revealed.

1. Introduction

Thermal-hydromechanical (THM) coupling processes in geotechnical media play an important role in a wide range of engineering applications. Many significant issues, such as resources mining (e.g., coal, geothermal energy, natural gas, and oil) [1, 2], traffic engineering (e.g., tunnel and metro) [3], and underground repository (e.g., chemotoxic and nuclear waste and CO₂ as well as natural gas sequestration) [4, 5], attracted the scholars’ attentions greatly.

Many scientific efforts have been exerted to reveal the THM interaction mechanism in the geotechnical systems. Since most geotechnical applications are characterized by long-term operating (several tens or hundreds of years) and large scales in size (several hundreds or thousands of meters in length, width, and depth), it is impossible to conduct in situ physical experiments. Therefore, mathematical models and simulation codes are emphasized by scholars. The THM coupling model originates from the isothermal hydromechanical (HM) coupling mechanism (also named fluid-solid interaction). The first HM coupled theory is the 1D consolidation theory of soil proposed by Terzaghi, followed by Biot’s 3D consolidation theory with isothermal and elastic consolidation [6]. By introducing the nonisothermal terms to the extended Biot’s equation [7] or using the averaging approach of the mixture theory, the basic THM coupling models are established [8]. In this interacting process, the rock is regarded as continuous mass points with characteristic parameters of both fluid and solid, which are governed by momentum, mass, and energy conservation laws. In this case, the porosity and fluid saturation are adopted to represent the storage capacity and mobility of fluid. Meanwhile, the elastic modulus, passion’s ratio, cohesive strength, internal frictional angle, and other parameters of rock are used to reflect the deformability of rock. As discussed in the international DECOVALEX (DEvelopment of COupled models and their VALidation against EXperiments), coupled conservation equations can be solved by the finite element method (FEM) or discrete finite element method, and sometimes both are used to handle the problem of fluid flow in fractures [9, 10]. On this basis, extensive improvements have been achieved, mainly concerning the fluid flow equation or solid constitutive model.

Theoretically, multiphase models based on Navier-Stokes equations in porous media are applied for THM coupling analysis in geotechnical materials, which covers the immiscible or miscible multiphase flow of Newtonian fluid or non-Newtonian fluid [11]. In the same way, verities of plastic rock models, such as the ideal elastic-plastic model [12], the viscoelastoplastic model [13], Mohr-Coulomb [14], Drucker-Prager [15], and the damage model [16], have been proposed with the purpose of acquiring the rock constitutive model which better reflect the features of natural rock. Many codes have also been employed to model and calculate these coupling equations, such as COMSOL, ABAQUS, RFPA, FLAC, UDEC, ROCMAS, THAMES, FRACON, COMPASS, FROCK and CODE_BRIGHT, FRT-THM, FLAC-TOUGH, and FRACture [17–19]. However, as tremendous grids of the models are required to construct the complex and disordered pore structure of rock, the above studies have significant shortcomings that are neglecting the fluid flow in the micropore space of rock.

Experiments on rock at micro- or mesoscale are more feasible in consideration of the size of the test sample. The mesoscopic THM coupled study refers to the fluid flow test on rock core under the condition of pressure and temperature obtained from experiments or numerical simulation. Nowadays, the uniaxial/triaxial test of rock at HTHP (High Temperature, High Pressure) is widely used in the laboratory to acquire mechanical properties [20]. A multiphase flow displacement device has been added to the triaxial system to test the variation of fluid transport properties [21]. Acoustic emission device is employed to monitor the cracking development during the loading process [22]. NMR can be adopted to investigate the movable fluid in this process [23]. However, these devices cannot be used at the same time. Micro-CT is another approach for investigating the status of fluid flow or solid deformation and crack [24], but the resolution is limited by the core holder (for pressure or temperature loading) [25]. In addition, it is difficult to conduct real-time monitoring according to the time requirements for a full scan [26]. Equipment with smaller size and better measuring accuracy is essential to satisfy the requirements of microscopic experiments. Thus, it is difficult to conduct the 3D real-time investigation and research on micromesoscopic THM coupled process at HTHP, in which case most existing studies of this kind are monitored by surface imaging technology, for example, SEM [27]. The 3D printing rock-like specimen has been adopted to print the rock matrix and investigate the inner fluid flow or crack development [28], but the microstructure and mechanical properties of substitutions are different from the natural rock. Numerical simulation based on micromesoscopic modelling of rock has been regarded as a platform to study the THM coupled process of rock.

The developments of the imaging technologies, such as nano- or micro-CT technology and SEM, make it possible to investigate the rock structure and minerals distribution at the resolution of micron or nanometer [29]. On this basis, pore-scale modelling has been emphasized as an effective means to conduct fluid flow or deformation simulation [30]. In literature, the pore-scale modelling methods can be classified into two categories: pore network model and grids model [31]. The former one is characterized by topologically representative network with idealized assumptions, which is efficient for the prediction of multiphase flow but is confined to fluid simulation only [32]. The second one is reconstructing the grids or mesh models from a binarized three-dimensional image, and then simulations are conducted using LBM or Navier-Stokes (NS) equations [33–35]. The volume tracking method including volume of fraction (VOF), Level Set (LS), or Level Set Method Progressive Quasistatic (LSMPQS) method is usually adopted along with the NS equations [36, 37]. Models of this kind can reproduce the image of rock but contain large amounts of elements because of the complex microrock structure, which in consequence are computationally demanding. Another advantage of this approach is that it is able to establish rock matrix models [38], which can be assembled with the corresponding pore model and used for the multiphysics coupled simulation. However, there are still some challenges in mesh quality controlling, modelling of the mineral distribution, and the remesh algorithm of fluid-solid interface since the shape of the rock structure varies in the THM coupled process. In mesoscopic aspects, the crack of the rock is the main concern. This allows the matrix in micropore models and crack to be established and simulated using the same methods as the microscopic simulation. However, most studies relevant to crack development and THM coupled process are limited to 2D [39].

This paper presents a fully THM coupled process in the rock using the reconstructed and assembled pore-scale models of rock matrix and pore space. Then the effects of stress and temperature on the pore structure, petrophysical properties, and water flooding efficiency are analyzed.

2. Pore-Scale Modelling and Boundary Conditions

The pore-scale models of both pore space and rock matrix are generated using the algorithm proposed in our previous paper [40]. The rock samples used in this paper are drilled from the original rock sample and scanned by Zeiss Xradia MICROXCT-400 of the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation in Southwest Petroleum University. A cube of voxels is extracted from the original micro-CT images and used as an input to the reconstructing process. The detailed information of the rock images used in this paper is listed in Table 1. The distribution of pore radius of the samples used in this study is shown in Figure 1.

Table 1: Rock sample image information.

Number
Rock type
Image resolution
μm/pixel
Size/pixel
Porosity
Number of computational elements

Matrix
Pore

B1
Berea sandstone
5.345
400³
19.65%
6428360
1571640

C1
Carbonate
3.314
400³
17.12%
6630176
1369824

MS1
Synthetic sandstone
2.055
300³
34.86%
2198473
1176527

S5
Sandstone from Shengli Oilfield, China
2.51
400³
12.11%
7031198
968802

S6
Sandstone from Shengli Oilfield, China
5.01
200³
40.34%
4772656
3227344

Figure 1: Pore radius distribution of images used in this paper.

As is shown in Figure 2(a), the reconstructed rock matrix and pore space models are assembled. In the simulation, the deformation of rock matrix is analyzed using ANSYS, while the fluid flow in the micropore of rock is simulated by CFX. The initial boundary conditions are applied in the Workbench platform. Taking the sample MS1 as an example, its boundary conditions of pore space and rock matrix are presented in Figures 2(b) and 2(c), respectively. The front and back surface of matrix along direction are fixed. The other surfaces are imposed on confining pressure. The pressure gradient is applied between the front and back surface of pore space along direction. The interface between rock matrix and pore is defined as FSI boundary. Other surfaces of pore space are defined as impermeable boundary. Constant temperature boundaries are defined for both solid and fluid parts.

(a) Assembled rock matrix and pore space model

(b) Boundary condition of pore space

(c) Boundary condition of rock matrix

Figure 2: Reconstructed model and boundary conditions of sample MS1.

Rock matrix is assumed as isotropic, homogenous, and ideal elastic-plastic, thus only a limited range of stress and temperature values are simulated in this study. The rock properties used in the simulation are presented in Table 2, in which the mechanical parameters of rock are tested by microindentation test of rock sample. And the fluid properties at 273 K are listed in Table 3.

Table 2: Rock properties of different sample.

Sample number
Property

Density
/kg·m⁻³
Elastic modulus
/GPa
Poisson’s Ratio
Yielding strength
/MPa
Thermal expansion/×10-5°C⁻¹

B1
2100
18.43
0.225
92.9
3.5

C1
2700
76.26
0.24
250
6

MS1
2300
14.19
0.31
81
5

S5
2675
20.13
0.28
73
6

S6
2500
9.35
0.29
67
6

Table 3: Fluid properties used in the simulation at the temperature of 273 K.

Fluid
Density
(kg/m³)
Viscosity
(cP)
Interfacial tension
(mN/m)
Contact angle ()

Drainage
Water flooding

Water
890
48
1

Oil
1200
1

3. Mathematical Model of THM Coupling in Rock

The mathematical model of multiphase flow in deformable rock contains two parts: governing equations of fluid flow and solid deformation.

3.1. Governing Equations of Fluid Flow

VOF (volume of fraction) model in CFX software is used to simulate the immiscible water and oil in the reservoir. The continuity equation for the th phase is [41]

where is volume fraction of fluid in the cell and is fluid’s density. When , the equation is just the continuity equation for the single flow.

The properties in the transport equations are determined by the volume fraction of the component phases in each cell:

All other properties (e.g., viscosity) are computed in this manner.

Navier-Stokes equation is used as the conservation of momentum for the fluid flow [41]:

The energy equation is shared by all the phases in the control volume and can be described as [41]

where is the effective thermal conductivity of phase and the energy and temperature are calculated by the mass-weighted average value of all phases:

where and represent the energy and temperature of phase, respectively.

The interfacial tension between two immiscible phases is unneglectable in the micropores of rock, which would lead to high capillary force. Here, the continuum surface force (CSF) model proposed by Brackbill et al. in 1992 [42] is used as follows:

In the CSF model, the phase interface curvature can be calculated by the local gradients of phase interface normal, which is determined by the volume fraction gradient of phase:

By the divergence theorem, the force on the interface can be transferred into the volume force. It has the following form:

where is the surface tension coefficient and is defined in terms of the divergence of the unit normal () of phase interface.

Considering the wall adhesion effect, the contact angle between the solid surface and the fluid is adopted to modify the unit normal () of phase interface nearby the surface

where is the unit vectors normal to the wall and is tangential to the wall, respectively, and is the contact angle. In the simulation process, the structured mesh model is divided into different parts. The contact angle, which follows a uniform distribution with a range of given interval, is assigned randomly to each part to obtain uniformly wet system. The wettability of the model is determined by the mathematical expectation of the given interval of the contact angle.

Using CFX software, the outlet flow rate can be acquired. Then the absolute permeability is calculated in the following term [43]:

Then the relative permeability is given by [43]

where is the total single-phase flow rate through the model and is the total flow rate of phase in multiphase conditions with the same imposed pressure drop. And both and can be acquired by the Fluent software.

3.2. Governing Equations of Rock Matrix Deformation

The three-dimensional equilibrium differential equation is

Here is the stress tensor and is the body force.

The three-dimensional geometric equations of the rock matrix are

Here is strain and is the displacement component.

The elastic physical equations are

and here is elastic modulus and is Poisson’s ratio.

4. THM Coupling Simulation

Based on the rock mesh model, the THM coupling mechanism in rock and its influence on water flooding process in the petroleum industry are analyzed using both ANSYS and CFX software. The fluid used in the single-phase flow simulation is water, and both oil and water are used for two-phase flow. In the CFX solver, a laminar flow is assumed. A transient model is used with the second-order backward Euler scheme. Automatic timestep and a convergence criterion of 10⁻⁶ are used. In the ANSYS solver, a transient structural solver is used to apply the boundary conditions of solid part and the default solver control is used. The mechanical input file will be generated and used as input to the ANSYS multifield solver.

5. Single-Phase Flow

5.1. Influences of Effective Pressure on Porosity and Permeability

Based on the structured mesh models of samples B1, C1, MS1, and S6, the evolution mechanism of effective pressure is analyzed. Taking sample MS1 as an example, under the condition of  Pa,  MPa, and the constant temperature of 273 K, the deformation displacement of both matrix model and interface of pore model is presented in Figures 3(a) and 3(b). Due to the complexity and inhomogeneity of microstructure of rock, the strain distribution is also characterized by inhomogeneity. Meanwhile, the fluid pressure and velocity distribution are shown in Figures 3(c) and 3(d). It indicates that water flows mainly along the channels with larger size and better connectivity to reduce the flowing resistance.

(a) Deformation displacement of matrix

(b) Deformation displacement of pore space

(c) Distribution of fluid pressure (Pa)

(d) Distribution of fluid velocity (m/s)

Figure 3: THM simulation of model MS1 under the condition that the confining pressure is 20 MPa and the pore pressure is 1000 Pa.

As is shown in Figure 4, the porosity and permeability decrease along with the rising of the confining pressure under the condition of constant pore pressure (i.e., the rising of the effective pressure). The rate of porosity drop becomes slower when the confining pressure is beyond 30 MPa, which means the plastic deformation occurs inside the rock matrix. The reason lies in that rock is assumed as isotropic, homogenous, and ideal elastic-plastic, and thus the strain would not change when the load is beyond its yield strength. Meanwhile, it shows a negative correlation between the drop rate of porosity and the elastic modulus of rock.

(a) versus confining pressure

(b) versus confining pressure

Figure 4: Variation of porosity and permeability with the rising of confining pressure and constant pore pressure. refers to the ratio of the porosity under the simulation condition to the origin sample, so it is with referring the permeability.

In addition, comparative analysis on the variation of permeability with the porosity under the same load is shown in Figure 5. It is found that the decline of the permeability is larger than the porosity. Taking model S6 as an example, when the porosity declines by 10%, the permeability drop ratio reaches almost 20%. This indicates that the permeability of rock is more sensitive to the confining pressure than the porosity.

Figure 5: Permeability variation versus porosity variation.

5.2. Effects of Temperature on Porosity and Permeability

In this section, the effects of temperature on the porosity and permeability under the condition of constant confining pressure (20 MPa) and pore pressure (5 MPa) are analyzed. The rock is assumed to be elastic under the temperature, and only the thermal expansion of rock in the temperature range of 20°C, 100°C is analyzed. The plastic deformation and the crack development under the temperature in the simulation are not considered. The thermal strain of model MS1 is presented in Figure 6, characterized by nonhomogeneous distribution.

Figure 6: Thermal strain distribution of model MS1 under the temperature of 100°C.

As is shown in Figure 7, the porosity and permeability of rock decrease with the rising of the temperature under the condition of constant confining and pore pressure, which is in contrast with the traditional experimental benchmark data, especially in a high temperature. This is because the rock matrix expands with the rising of the temperature, but the expansion would also lead to microcrack, which will enlarge the porosity and permeability of rock, too. When the temperature is varied from 20°C to 100°C, the porosity drop is less than 2%, but the maximum permeability drop is about 5%, which means the permeability is also more sensitive to the temperature than the porosity.

(a) versus temperature

(b) versus temperature

Figure 7: Variation of porosity and permeability with temperature.

Based on the simulation, the permeability variation of model MS1 with the effective pressure and temperature along , , and direction is presented in Figure 8. In the simulation, the confining pressure is variable and the pore pressure is constant to be 5 MPa. It can be found that the permeability and its drop rate vary for different direction with the rising of effective pressure and temperature. Considering the isotropic assumption of rock, the complex and disorder structure of rock is the main reason to this phenomenon.

(a) Permeability versus effective pressure and temperature along direction

(b) Permeability versus effective pressure and temperature along direction

(c) Permeability versus effective pressure and temperature along direction

Figure 8: Permeability variation of model MS1 with effective pressure and temperature along , , and direction.

6. Two-Phase Flow

Considering that the oil solubility in water is small enough to neglect, the VOF (volume of fraction) model is used to simulate the immiscible displacement process between water and oil in the reservoir. The core sample is initially saturated with water to represent the original stratum without oil. Then the first-cycle oil flooding is proceeded to represent the formation of oil, in which process the core sample becomes more oil-wet. After that, water injection is simulated to represent the water flooding development of the reservoir. The oil distribution of MS1 after the first cycle of oil flooding process and the second cycle of water flooding process is shown in Figure 9.

(a) After first cycle of oil flooding process

(b) After second cycle of water flooding process

Figure 9: Oil volume fraction of MS1 after oil flooding and water flooding process; the red parts represent the oil phase.

6.1. Influences of Effective Pressure on Water Flooding Efficiency

The effects of confining pressure on water flooding efficiency under the condition of constant pore pressure (1000 Pa) and temperature (20°C) are analyzed. The variations of relative permeability curves of model S5 and MS1 are presented in Figure 10. It is found that, with the rising of the confining pressure (i.e., the rising of the effective pressure), the relative permeability of both water and oil decreases, which means the decline of the fluid mobility. In this case, the residual oil saturation increases. The reason lies in that the size of pore space (especially the throat) decreases with the rising of the effective pressure, which leads to a higher capillary pressure and a lower oil recovery. Thus, higher pressure of induced water contributes to EOR in the water flooding process.

(a) Relative permeability curves of S5 for different confining pressure

(b) Relative permeability curves of MS1 for different confining pressure

Figure 10: Effects of stress on the relative permeability.

6.2. Effects of Temperature on Water Flooding Efficiency

As is shown in Figure 11, with the rising of temperature, the relative permeability of both water and oil increases. Though the pore size reduces caused by the expansion of the rock matrix with the rising of temperature, the decrease of water and oil viscosity improves the fluid mobility. Meanwhile, the decline of oil-water mobility ratio reduces the fingering effect in the displacing process, which promotes the sweep efficiency and oil recovery as well. Thus, high temperature of induced water is beneficial to enhance oil recovery, especially for heavy oil with high viscosity.

(a) Relative permeability curves of S5 for different temperature

(b) Relative permeability curves of MS1 for different temperature

Figure 11: Effects of temperature on the relative permeability.

7. Conclusion

In this paper, a pore-scale study on the thermal-hydromechanical coupling simulation of porous rock is conducted. Based on the structured mesh models of rock matrix and pore space, the effects of stress and temperature on the microstructure, porosity, permeability, and relative permeability in the linear elastic and linear thermal-expanding process are analyzed. The results indicate that the rising of effective pressure or temperature would lead to the decline of the porosity and permeability, and the drop ratio of permeability is larger than that of porosity. The relative permeability of oil and water decreases with the increasing of the effective pressure, so it is with the oil recovery. However, the relative permeability of the two phases and oil recovery increase as a result of the fluid mobility improvement by the rising of temperature. Thus, high temperature and high pressure of induced water are beneficial to enhance oil recovery, especially for heavy oil with high viscosity. Though the petroleum industry is the main concern in this paper, the outcomes can be applied to other kinds of THM coupling process of porous media.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This paper is financially supported by National Science and Technology Major Project of China under Grant no. 2017ZX05013001-002.

References

J. Rutqvist, “Fractured rock stress-permeability relationships from in situ data and effects of temperature and chemical-mechanical couplings,” Geofluids, vol. 15, no. 1-2, pp. 48–66, 2015.

S. A. Bea, U. K. Mayer, and K. T. B. Macquarrie, “Reactive transport and thermo-hydro-mechanical coupling in deep sedimentary basins affected by glaciation cycles: model development, verification, and illustrative example,” Geofluids, vol. 16, no. 2, pp. 279–300, 2016.

Z. X. Sun, X. Zhang, Y. Xu et al., “Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model,” Energy, vol. 120, pp. 20–33, 2017.

D. Wu, Y. Zhang, R. Zhao, T. Deng, and Z. Zheng, “A coupled thermal-hydraulic-mechanical application for subway tunnel,” Computers & Geotechnics, vol. 84, pp. 174–182, 2017.

O. Kolditz, S. Bauer, C. Beyer et al., “A systematic benchmarking approach for geologic CO₂ injection and storage,” Environmental Earth Sciences, vol. 67, no. 2, pp. 613–632, 2012.

F. Dupray, C. Li, and L. Laloui, “THM coupling sensitivity analysis in geological nuclear waste storage,” Engineering Geology, vol. 163, pp. 113–121, 2013.

M. F. Kanfar, Z. Chen, and S. S. Rahman, “Fully coupled 3D anisotropic conductive-convective porothermoelasticity modeling for inclined boreholes,” Geothermics, vol. 61, pp. 135–148, 2016.

L. Jing, “A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 3, pp. 283–353, 2003.

M. Hassanizadeh and W. G. Gray, “General conservation equations for multi-phase systems: 1. Averaging procedure,” Advances in Water Resources, vol. 2, pp. 131–144, 1979.

L. Jing, O. Stephansson, C. F. Tsang, L. J. Knight, and F. Kautsky, DECOVALEX II Project Executive Summary, No. SKI-R–99-24, Swedish Nuclear Power Inspectorate, 1999.

M. Sahimi, Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches, John Wiley & Sons, New York, NY, USA, 2011.

K. R. Rajagopal, A. Z. Szeri, and W. Troy, “An existence theorem for the flow of a non-newtonian fluid past an infinite porous plate,” International Journal of Non-Linear Mechanics, vol. 21, no. 4, pp. 279–289, 1986.

J. Rutqvist, D. Barr, R. Datta et al., “Coupled thermal-hydrological-mechanical analyses of the Yucca Mountain Drift Scale Test—comparison of field measurements to predictions of four different numerical models,” International Journal of Rock Mechanics and Mining Sciences, vol. 42, no. 5-6, pp. 680–697, 2005.

T. V. Gerya and D. A. Yuen, “Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems,” Physics of the Earth and Planetary Interiors, vol. 163, no. 1–4, pp. 83–105, 2007.

P.-Z. Pan, X.-T. Feng, X.-H. Huang, Q. Cui, and H. Zhou, “Coupled THM processes in EDZ of crystalline rocks using an elasto-plastic cellular automaton,” Environmental Geology, vol. 57, no. 6, pp. 1299–1311, 2009.

L. Laloui, M. Nuth, and L. Vulliet, “Experimental and numerical investigations of the behaviour of a heat exchanger pile,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 30, no. 8, pp. 763–781, 2006.

L. C. Li, C. A. Tang, S. Y. Wang, and J. Yu, “A coupled thermo-hydrologic-mechanical damage model and associated application in a stability analysis on a rock pillar,” Tunnelling and Underground Space Technology, vol. 34, pp. 38–53, 2013.

O. Stephanson, L. Jing, and C. F. Tsang, Coupled Thermo-hydro-mechanical Processes of Fractured Media: Mathematical and Experimental Studies, vol. 79, Elsevier, Amsterdam, The Netherlands, 1997.

J. Taron, D. Elsworth, and K.-B. Min, “Numerical simulation of thermal-hydrologic-mechanical-chemical processes in deformable, fractured porous media,” International Journal of Rock Mechanics & Mining Sciences, vol. 46, no. 5, pp. 842–854, 2009.

B. Vásárhelyi and A. Bobet, “Modeling of crack initiation, propagation and coalescence in uniaxial compression,” Rock Mechanics and Rock Engineering, vol. 33, no. 2, pp. 119–139, 2000.

H. T. Hall, “Some high-pressure, high-temperature apparatus design considerations: equipment for use at 100 000 atmospheres and 3000°C,” Review of Scientific Instruments, vol. 29, no. 4, pp. 267–275, 1958.

G. J. Fischer and M. S. Paterson, “Measurement of permeability and storage capacity in rocks during deformation at high temperature and pressure,” International Geophysics, vol. 51, pp. 213–252, 1992.

T.-F. Wong, C. David, and W. Zhu, “The transition from brittle faulting to cataclastic flow in porous sandstones: mechanical deformation,” Journal of Geophysical Research B: Solid Earth, vol. 102, no. 2, pp. 3009–3025, 1997.

H. Li, H. Guo, S. Zhou, Z. Meng, and X. Wang, “NMR analysis of movable remaining oil of low pemeability reservoir,” Journal of Southwest Petroleum University (Science & Technology Edition), vol. 38, no. 1, pp. 119–127, 2016 (Chinese).

S. Iglauer, A. Paluszny, C. H. Pentland, and M. J. Blunt, “Residual CO₂ imaged with X-ray micro-tomography,” Geophysical Research Letters, vol. 38, no. 21, pp. 1440–1441, 2011.

T. Bultreys, M. A. Boone, M. N. Boone et al., “Fast laboratory-based micro-computed tomography for pore-scale research: illustrative experiments and perspectives on the future,” Advances in Water Resources, vol. 95, pp. 341–351, 2015.

P. Liu, Y. Ju, P. G. Ranjith, Z. Zheng, and J. Chen, “Experimental investigation of the effects of heterogeneity and geostress difference on the 3D growth and distribution of hydrofracturing cracks in unconventional reservoir rocks,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 541–554, 2016.

K. Hadley, “Comparison of calculated and observed crack densities and seismic velocities in westerly granite,” Journal of Geophysical Research, vol. 81, no. 20, pp. 3484–3494, 1976.

F. M. Ezzein and R. J. Bathurst, “A transparent sand for geotechnical laboratory modeling,” Geotechnical Testing Journal, vol. 34, no. 6, pp. 590–601, 2011.

T. Bultreys, W. De Boever, and V. Cnudde, “Imaging and image-based fluid transport modeling at the pore scale in geological materials: a practical introduction to the current state-of-the-art,” Earth-Science Reviews, vol. 155, pp. 93–128, 2016.

R. Hilfer and T. Zauner, “High-precision synthetic computed tomography of reconstructed porous media,” Physical Review E, vol. 84, no. 6, Article ID 062301, 2011.

M. J. Blunt, B. Bijeljic, H. Dong et al., “Pore-scale imaging and modelling,” Advances in Water Resources, vol. 51, pp. 197–216, 2013.

P. V. Marcke, B. Verleye, J. Carmeliet, D. Roose, and R. Swennen, “An improved pore network model for the computation of the saturated permeability of porous rock,” Transport in Porous Media, vol. 85, no. 2, pp. 451–476, 2010.

P. Mohammadmoradi and A. Kantzas, “Pore-scale permeability calculation using CFD and DSMC techniques,” Journal of Petroleum Science and Engineering, vol. 146, pp. 515–525, 2016.

Z. Liu and H. Wu, “Pore-scale study on flow and heat transfer in 3D reconstructed porous media using micro-tomography images,” Applied Thermal Engineering, vol. 100, pp. 602–610, 2016.

A. Ferrari and I. Lunati, “Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy,” Advances in Water Resources, vol. 57, no. 9, pp. 19–31, 2013.

P. Mohammadmoradi and A. Kantzas, “Petrophysical characterization of porous media starting from micro-tomographic images,” Advances in Water Resources, vol. 94, pp. 200–216, 2016.

J. Liu, R. Song, and M. Cui, “Numerical simulation on hydromechanical coupling in porous media adopting three-dimensional pore-scale model,” The Scientific World Journal, vol. 2014, Article ID 140206, 8 pages, 2014.

S. Q. Yang, Strength Failure and Crack Evolution Behavior of Rock Materials Containing Pre-Existing Fissures, Science Press, Beijing, China, 2015.

R. Song, J. Liu, and M. Cui, “A new method to reconstruct structured mesh model from micro-computed tomography images of porous media and its application,” International Journal of Heat and Mass Transfer, vol. 109, pp. 705–715, 2017.

ANSYS User's Guide, ANSYS, 2012.

J. U. Brackbill, D. B. Kothe, and C. Zemach, “A continuum method for modeling surface tension,” Journal of Computational Physics, vol. 100, no. 2, pp. 335–354, 1992.

H. A. Sheldon and A. Ord, “Evolution of porosity, permeability and fluid pressure in dilatant faults post-failure: implications for fluid flow and mineralization,” Geofluids, vol. 5, no. 4, pp. 272–288, 2005.