Abstract
Studies of rock stress sensitivity are mainly focused on experimental and data processing methods, and the mechanism cannot be adequately explained using specific pore shape models. This study, based on a random pore network simulation, explains the rock stress sensitivity mechanism for the first time. Based on the network model theory, the hydraulic conductivity equation, the dimensionless radius equation, and the effective stress equation for partially saturated rock are used to generate a three-dimensional random pore network model based on the QT platform. The simulation results show that the influence of the effective stress on the dimensionless radius becomes more significant as the aspect ratio decreases, and the relationship between dimensionless radius and effective stress can be effectively interpreted through different combinations of pore shapes. Moreover, the mechanism behind permeability stress sensitivity can be explained by establishing the relationship between permeability and effective stress.
1. Introduction
Rock stress sensitivity refers to the changes in rock petrophysical parameters caused by effective stress, including porosity stress sensitivity [1–3], permeability stress sensitivity [4, 5], and stress sensitivity of electrical resistivity [1, 6], among which permeability stress sensitivity is most frequently discussed [4, 5]. The study of permeability stress sensitivity is mainly focused on experimental and data processing methods [5], which lack in-depth analysis of the micro-pore structure eigenvalue.
Although some researchers have previously built theoretical models based on specific pore shapes, such as circular [7], oval [8], conical, hyperbolic triangular [9], and hyperbolic quadrilateral [10], these models, in which only specific pores are considered, are not sufficient to explain rock stress sensitivity. Therefore, pore network models that consider various pore shapes are preferable to study the mechanisms applicable to the micro-pore structure [11].
Pore network models include pore space-based imaging network models [12] and regular pore network models [11]. Use of the pore space-based imaging network model is relatively computationally expensive [12], whereas regular pore network models set the pore size and distribution and adopt certain distribution functions by integrating percolation theory, which results in the pore network model featuring the same complexity as that of a real rock [13]. As a result, this study, based on a QT platform, focuses on the dimensionless radius-randomized stress change by adopting a C++ program to generate randomized pore network models and then integrating the different pore types, proportions, and micro-pore networks. The permeability stress sensitivity mechanism is then explained by establishing the relationship between permeability and the dimensionless radius function.
2. Establishing the Randomized Pore Network Model
2.1. Theoretical Basis of Network Models
2.1.1. Percolation Theory
Percolation theory is central to the study of randomized pore network models, which is mainly used to describe random structures and connectivity. It is used to study the fluid distribution and flow in a random and disorderly medium, as in randomized pore network models using statistical methods. The randomized pore network model employed in this study is designed as follows: the lines in the model symbolize throats of a certain volume and flow resistance with different distribution modes such as a uniform distribution and the intersections of the lines symbolize the throats without volume and flow resistance that only function as connections; therefore, the calculation of percolation parameters should be mainly focused on the linearity calculation without considering the node. Such a technique, based on the linear distribution of the network, studies the impact of the microdistribution on the macroproperties of the porous media. Since the percolation theory complies with the probability theory and the statistical accuracy depends on the sample size, the samples should meet the requirements for a reliable statistical result, and therefore the network simulation based on percolation theory requires that the three-dimensional network model has nodes of at least 20 × 20 × 20 [13]. Using these conditions and the allowable computer performance, the models can be applicable to a broader range of situations. Within the calculation limits of the computer, the more nodes the models have, the more macroattributes they will reflect.
2.1.2. Similarity Principle between Water and Electricity
The network simulation is based on the similarity principle between water and electricity. By analyzing the flow in the porous medium, the network structure of the circuit can be used to conduct a simulation analysis. The current in the circuit follows Ohm’s law:where is the current, is the resistance, and is the voltage.
By analyzing the network circuit, we assume that the circular cross section pipes filled with conduction fluid are made of resistors: where is the pipe radius, is the pipe electrical resistivity, and is the resistor or pipe length.
Assuming that flow is laminar, the simplest cylindrical pipe is taken as an example and the Poiseuille equation is used:where is the flow, is the fluid viscosity, and is the differential pressure at both ends of the pipe.
The circuit network and fluid pore network are all made up of pipes. According to (2) and (3), Ohm’s law and the Poiseuille equation both explain the relationship between the flow of a current or fluid and the differential pressure at both ends of the pipes (i.e., voltage and flow differential pressure values), the relationship between pipes, and the differential pressure of both ends of pipes (they are voltages and flowing pressure differential value), and the relationship among pipes. The flow relationship is similar to the similarity principle between water and electricity.
Based on this, the Poiseuille equation can be simplified to Ohm’s law: where .
It can be seen from the similarity principle between water and electricity above that the fluid in the porous network is similar to the current in the circuit network; therefore, we can conduct the fluid flow simulation in the pore network model based on the current in the circuit network and use the analytical method for the current to analyze the simulated fluid model.
2.1.3. Kirchhoff Law
Kirchhoff’s circuit laws are divided into (1) Kirchhoff’s current law, whereby the current inflows of any node in the circuit are equal to the current outflows, and (2) Kirchhoff voltage law, whereby the voltage in any circuit loop shall be zero after completing the loop in the direction of current flow. For the pore network model, Kirchhoff’s current law is adopted to build the equation set of the model.
2.2. Key Points of the Program
2.2.1. The Equation Deduction of the Effective Stress for Partially Saturated Rock
This study is focused on the quasi-static, two-phase flow, random pore network model. The fluid flow in the model is fully controlled by the capillary pressure, assuming that the fluid is incompressible and the influence of a viscous force is ignored. There are several assumptions reflected in the physical model. The oil-water or gas-water interface is relatively static; namely, the fluid distribution of the corresponding pore throat unit will change once displacement occurs. This assumption is the same as the low-speed percolation of most porous media. Furthermore, in the network simulation we generally assume that the displacement process is instantly completed and then balanced; namely, the nonwetting fluid pressure () is greater than or equal to the inlet capillary pressure of a throat (), and therefore the displacement action will start and the nonwetting fluid will displace the wetting fluid in the throat. The throats without displacement are filled with wetting fluid and the pressure is kept constant (zero). Under such circumstances, the capillary pressure of the displaced throat is as follows:
Throats with displacement contain a type of flowing fluid whereas the wetting fluid at the corner is taken as a static and nonflowing fluid; its pressure value remains at 0, the flow inside the pipes can be deemed as single-phase flow, and the fluid pressure () equals the nonwetting displacement pressure (). For throats that have not been displaced, the fluid inside is the wetting phase (water) and the pressure is 0. There is no pressure difference between the nonwetting and wetting phases and capillary pressure does not exist. Therefore, the effective stress expression is as follows:where is the effective stress and is the overburden pressure.
This theory was proposed by Bishop [14], who experimentally defined the general effective stress for a partially saturated porous medium, as follows:
When wetting-phase saturation values are at the extreme endpoints (0 and 1), the product of wetting-phase saturation and capillary pressure () is 0, and when a relatively low value is obtained, for example, 0.1 or less, the peak value of will occur. However, this value is still very small, so in most cases it is negligible when compared with the effective stress [14]. Therefore, the expression above can be simplified as
This expression is the effective stress equation for two-phase fluids and is identical to (6); however, the value of the coefficient has been a topic for debate [15, 16]. If the value of the coefficient is 0, the effective stress is the confining pressure. If the effective stress coefficient is 1, the effective stress is the net stress. To simplify this study, the effective stress coefficient is set as 0, under which circumstances the pore fluid pressure is assumed to impart no pressure on the effective stress or has much less impact than the confining pressure. Therefore, this study on rock micro-pore structure is focused on the change in rock pore structure under different confining pressures.
2.2.2. Key Equations
The key equation of this method involves the hydraulic conductivity at the throats and the dimensionless radius equation. The flow equation for a hyperbolic triangle throat [9] iswhere FT is the dimensionless conductivity per unit length of unit area, which is solved by finite difference. The hydraulic conductivity equation is
Because there is no dimensionless radius equation for hyperbolic triangle throats in the literature, here we present a simple deduction based on the formula from Gangi [17]:
Moreover, it is known from another Gangi formula that [17]and considering thatit can be deduced that the dimensionless radius equation for hyperbolic triangle throats is
In a similar manner, the hydraulic conductivity equation and dimensionless radius equations for circular throats, oval throats, conical throats, hyperbolic triangle throats, and star throats can be deduced, shown in Table 1.
In these equations, is the dimensionless radius, which is the ratio of the pore throat radius under effective stress to the pore throat radius under no effective stress, FT is the dimensionless conductivity per unit length of unit area, which is solved by finite difference, is the Poisson ratio, is the length of the semiminor axis, is the length of the semimajor axis, is the aspect ratio, , and is Young’s modulus.
The dimensionless radius equations for hyperbolic triangle throats and star-shaped throats are similar; only the coefficient and index are different. Therefore, their conclusions are similar, which is illustrated by Figure 1 and in the following discussion.
2.2.3. Coefficient Matrix Solution
A key part of the network model simulation is building operation expressions, the basic idea of which is based on Kirchhoff’s law, nodes, and line conductivity. Noble and Daniels [18] presented the equation for a network of nodes: whereas .
For a simple network, the coefficient matrix can be obtained through a structural analysis; however, the simulation method of the model is based on statistical theory, which dictates that the number of nodes and lines in the network model needs to be large enough to ensure the reliability of the simulation. As a result, when the number of nodes and connections is large enough, the model cannot be solved by simple algebra, and the Cholesky decomposition is required in order to solve the coefficient matrix . However, the Cholesky decomposition affects the precision of the coefficient matrix, mainly due to the round-off error, and furthermore the program code is complicated [11]. To this end, we use the iterative method to solve the matrix.
The principal of the iterative method, which gradually approaches the real solution, is to take the assumed value as the solution and perform continuous iterations until it meets the convergence condition and obtains the solution of the equation. For the convergent system set, the deviation obtained from each displacement will decrease and its solution becomes closer to the real solution. The iteration method can also automatically adjust the occasional calculation error that occurs in the iteration. The method includes simple iteration and super-relaxed iteration and, in the following section, the simplest two-dimensional square network model is used as an example to introduce the solution process.
The basis of the simple iteration, also known as successive iteration, is to construct the fixed-point equation in order to obtain the approximate solution. The simple iteration method is solved via the following steps.
(1) Assign Initial Values for Nodes. As the ultimate solution is not associated with the initial given values, the initial value can be randomly set, yet the convergence rate depends on the accuracy of the initial values. We assume that the flow direction is from left to right in the square network, the left end of the nodes is assigned with voltage , and the right end of the nodes is assigned with voltage , while other nodes are assigned with a voltage of 0.
(2) Establish the Equation. Taking the 0 node in Figure 2 as an example, the following equation can be obtained according to Kirchhoff’s current law:
This equation can be used for every node in the network and, in this way, we obtain the system of linear equations whose quantities are identical to the node quantities:
To conduct the iteration solution process, (17) is changed as follows:
(3) Iteration Solution. The input/output current value in the network is calculated after each iteration. The iteration is only deemed completed if the input/output currents are equal or the difference is within the error. After the iteration is finished, the input/output current is solved, the resistance is calculated using Ohm’s law, and then the other parameters of the network model, such as electrical resistivity, are solved.
As the convergence rate of the simple iteration is low, we use the super-relaxed iteration, and the program obtains the adjacent node voltage values from the last step calculation during the node calculation. For example, when we calculate the points in the square network, the voltage of the left point and the point below will be replaced by the voltage values obtained from the last step:
The method is called Gauss-Seidel iteration, and it replaces the new values obtained from the previous last step in order to speed up the convergence rates. The increment can be written as
To speed up convergence and introduce the relaxing factor (obtaining the values 1-2), the above expression can be written as follows:
(4) Initial Radius Assignment. Considering the randomness of the rock pore radius and the strong heterogeneity of the reservoir rocks, the network models used in our research determine the pipe radius by random functions in order to generate randomized and nonhomogeneous network models. This study selects radius as per the logarithmic mean distribution and uniform distribution, of which the normalized standard deviations are 0.05, 0.30, 0.55, 0.80, and 1.05 and 0.05, 0.30, and 0.55, respectively. The radius is generated through the following random numbers:where is the random number, is the maximum throat radius, and is the minimum throat radius.
It can be seen from the above equation that the radius of any throat in the network is between the maximum throat radius and minimum throat radius. It should be particularly noted that the hydraulic radii (i.e., twice the volume to surface ratio of the pores, , where is the total pore volume and is the total wetted surface area) all remain constant (40 μm). Under the conditions of log uniform distribution and uniform distribution, the hydraulic radius and normalized standard deviation should meet the requirements in Table 2. When the hydraulic radius is set at 40 μm, the normalized standard deviation and the largest/smallest throat radius in Table 3 should be used.
2.3. Program Implementation
2.3.1. Input/Output Parameters
Based on the theory of the network model and the key programming points, the C++ language is used to compile random network model procedures on the QT platform. Regarding the grid size (, , and direction), the maximum/minimum throat radii, maximum/minimum aspect ratios, formation water viscosity, grid unit length (pipe length), connectivity probability, and the proportion of pores of each shape, aspect ratios are inputted through program interfaces.
When running the program, the first output is a 3D cube random network model, whose size is set at 100 × 100 × 100 (Figure 3), followed by the relevant parameters of this program (Table 4). Finally, it produces the aspect ratios of different shaped throats and graphs of the dimensionless radius and the effective stress.
2.3.2. Program Result Analysis
In order to analyze the influence of the aspect ratio on the calculation result, the program is set according to a specific throat, and the change in dimensionless radius due to effective stress is studied by setting different aspect ratios. See Figures 4–7 for the comparison results of the program valuation.
It is observed from Figures 4–7 that the smaller the aspect ratio, the more significant the influence of the effective stress on the dimensionless radius. Taking star-shaped throats as an example, we can see in the dimensionless radius equation that the smaller the aspect ratio is, the more the dimensionless radius will be influenced by the stress change. The star-shaped throats experience deformation with the change in stress (Figure 8), and when the aspect ratio is smaller, the throat approaches a crack shape, and the influence of the crack on stress is more evident.
(a)
(b)
The relationship described above is between the dimensionless radius and the effective stress in a specific throat and does not consider the range in aspect ratio values or the combination of the aspect ratio and the proportion. Therefore, taking the star-shaped pore as an example, the plan and operation results shown in Figure 9 consider the specific throat, aspect ratio, and proportion combinations; however, there will certainly be different plans for oval, cone, and trilateral throats.
In addition to the different aspect ratios and proportion combinations for one specific pore shape, other relationships exist between dimensionless radius and effective stress with combinations of two-, three-, or four-pore throat types. Due to limited space, the combination of three different throat types is taken as an example, and the relationship between dimensionless radius and effective stress is effectively explained by a specific combination plan (Figure 10).
3. Explanation for Rock Stress Sensitivity
The analysis above indicates that the pore structure of a rock changes with the degree of effective stress, which induces changes in the rock physical properties. Therefore, a crucial factor in studies of rock stress sensitivity, particularly permeability stress sensitivity, is understanding the response of the pore radius to effective stress; namely,
The relationship between the pore radius derived from the model and the dimensionless radius is as follows:which establishes the function for the pore radius (radius) and the effective stress (). It can be seen from Table 1 that the hydraulic power conductivity () is a function of the radius (radius) and the input/output flow is a function of the hydraulic power (). Furthermore, the permeability calculation function (which is based on Darcy’s Law) is related to the inflow or outflow, from which we establish the relationship between permeability and effective stress (). Therefore, our model, by fitting data from the top of the rock core through a combination of different pore shape models (see Figure 11), can effectively explain the mechanisms behind rock permeability stress sensitivity.
4. Conclusion
(1)The effective stress equation for partially saturated rock is derived and verified, and the core operation expressions of the network model are obtained using iteration methods.(2)The network model simulation shows that a circular throat is the special case for an ellipse, and the smaller the aspect ratio, the greater the effect of stress on the dimensionless radius.(3)The relationship between the experimental dimensionless radius and the effective stress can be explained through different pore shape combinations based on the network model.(4)The permeability stress sensitivity is effectively explained using the network simulation.
Competing Interests
The authors declare that they have no competing interests.