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 CO2 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.
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
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.
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 isHere is the stress tensor and is the body force.
The three-dimensional geometric equations of the rock matrix areHere is strain and is the displacement component.
The elastic physical equations areand 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−6 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)
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
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.
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.
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
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
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
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
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
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.