Abstract
In this study, the hydromechanical behavior was experimentally investigated by conducting fracture geometry evolution-based gas flow tests on single fractured granite. The traditional fracture roughness value () of a single surface cannot adequately explain the intrinsic permeability in multiple samples. Instead, the permeability is better explained by aperture distribution, and it is positively correlated with the mean aperture size. The permeability is also affected by a combination of contact area and void space under changing confining stress. The average embedded depth increases rapidly under low-loading stress (less than 8 MPa) and slowly under high-loading stress (higher than 8 MPa). A simplified Hertz contact model was used to fit the experimental data. The model uses the reciprocal of the standard deviation of the aperture as the average curvature radius of the fracture. Good agreement was found between the experimental and theoretical results. A modified smooth parallel-plate model to describe flow friction in fractures was developed, which takes fracture contact and void space into account. A new friction model was also developed that takes both the nonlinear effect and the influence of fracture relative roughness into account by multiplying the ideal model by two power-law functions. This new model is effective in predicting the friction factor. Finally, taking into consideration the built permeability, mechanical, and friction models, a nonlinear flow model was proposed. The flow rate can be estimated based on the 3D fracture topography data, the applied loading stress, and the pressure gradient. The influence of four types of aperture distribution on the friction factor was investigated. The results showed that a high mean aperture size would result in a low friction factor. In addition, the influence of the mean aperture size on the friction factor is obvious, even under a very high-loading stress while the effect of the standard deviation is negligible when the loading stress is very high.
1. Introduction
An accurate understatement of fluid flow through fractured rock is crucial in many science and engineering situations, including extraction of geothermal energy, recovery of oil and natural gas, geologic sequestration of CO2, development of unconventional shale gas, coalbed methane reserves, and storage of radioactive waste [1–6]. One area of this field that attracts significant interest is the hydromechanical responses of rock fractures. Over the course of completing experimental and theoretical studies [7–14], researchers have proposed several models describing fluid flow and transport in fractured rocks. During the flow process, pressure head loss occurs due to flow acceleration and friction with the walls. The dimensionless friction factor is an important parameter used to describe pressure drops within rough fractures [11]. Nazridoust et al. [10] investigated the friction factor for laminar, single-phase flow through a 2D rock fracture channel by considering different apertures using computational fluid dynamics (CFD). The authors proposed a new friction factor model based on the results and best-fit regression. Chen et al. [15], Qian et al. [16], and Tzelepis et al. [17] experimentally investigated the effects of flow regime and fracture roughness on the friction factor. Furthermore, Zhang and Nemcik [12] experimentally investigated the influences of different fracture asperities on the friction factor. They suggested a model optimized for linear flow regimes with low Reynolds numbers (). More recently, Zhou et al. [18] proposed a theoretical formula form of the friction factor containing both viscous and inertial terms that is formulated by incorporating the Forchheimer equation. They also proposed a new friction factor model based on a recent phenomenological relationship for the Forchheimer coefficient. The inertial term is considered as a power function of the relative roughness.
Surface roughness is one of the most important factors affecting the hydromechanical behavior of rock fractures. The joint roughness coefficient parameter (JRC), first introduced by Barton [19], has been widely used in geotechnical and rock engineering applications. Because the JRC is estimated subjectively by visually comparing the observed roughness profiles to the 10 standard roughness profiles [20], the results always vary based on the experience level of the investigator. To overcome this shortcoming, researchers have developed methods to calculate JRC from other parameters derived from digitized profiles. These parameters include the first-derivative root-mean-square [21] as well as the roughness profile index , sometimes referred to as the profile sinuosity [22]. However, a rock fracture consists of two rough wall surfaces. Its hydromechanical behavior is controlled by interactions between the two surfaces [23]. More specifically, its mechanical properties are mainly controlled by the shape, size, and number of contacts between the surfaces, and its transport properties are determined by the separation (aperture) between the surfaces [24].
In a recent study, Xia et al. [23] introduced the composite morphology parameter to examine the nonlinear fluid flow characteristics through real rock joints under different contact conditions. is defined as the root-mean-square height of the composite morphology. Aperture distribution based on the composite morphology of rock joints, another factor affecting fracture hydromechanical behavior, has also received attention in recent years. Typically, the aperture distribution for a single fracture is generated based on the assumption of an approximately Gaussian distribution or an approximately log-normal distribution [25]. Other methods to obtain aperture distribution mainly depend on direct or indirect experiments, including casting [26, 27], NMRI [28], the X-ray CT technique [29, 30], and laser scanner-based methods [31–33]. Yeo et al. [26] carried out radial and unidirectional flow experiments with aperture replicas of natural sandstone fractures to study the effect of shear displacement at 0, 1, and 2 mm on the aperture size distribution. The experiments demonstrated that both the number of contact points and the fractional contact area decreased with increasing shear displacement. In addition, the mean aperture and standard deviation increased, and the ratio of standard deviation to mean aperture increased slightly. Koyama et al. [34] investigated the evolution of aperture and transmissivity with increasing shear displacement, taking the sample size effect into account. The fracture aperture was found to increase anisotropically during shear, with a more pronounced increase in the direction perpendicular to shear displacement. This resulted in a significant fluid flow channeling effect. Sharifzadeh et al. [32] used geographic information system (GIS) technology, which can obtain accurate joint surface morphologies, to determine and visualize the aperture distribution of a jointed granite block. Using GIS, the authors could visualize and interpret changes in aperture distributions at different normal stresses and shear displacements. Substantial studies have been conducted on the effects of mechanical behavior, including normal loading and shearing, on aperture distribution and flow behavior in rock joints; some experimental studies also have investigated the influence of aperture distribution on hydromechanical behavior. Kulatilake et al. [27] adopted the traditional Wood’s metal casting method to determine the spatial distribution of the aperture under normal stress for natural rock fractures. New functional relations were developed among aperture fractal dimensions, normal stress, mean aperture, and fluid flow rate per unit hydraulic gradient per unit width. Since the hydromechanical properties of a rock fracture depend strongly on the aperture distribution of the fracture [27], further study on this topic is merited.
Gas has the properties of low viscosity and easy diffusibility, so when gas flow through rock fractures, a higher Reynolds number occurs and nonlinear flow regimes, including transition and turbulent flow, generate more easily. Gas flow through rock fracture research can also be beneficial to understand supercritical carbon dioxide fracturing technology used in an Enhanced Geothermal System (EGS) [35, 36], since supercritical carbon dioxide has the same property of low viscosity and easy diffusibility [37]. In the present study, nitrogen flow tests were conducted through granite rock fractures under different loading stresses and pressure gradients. During the experiment, the aperture distribution was calculated based on 3D surface topography data obtained using a laser scanner. And then, based on the experimental flow data, numerical simulation was carried out and fracture geometry parameters such as contact area, void space, and embedded depth were calculated during loading. These fracture geometry parameters were then used to explain the hydromechanical behavior of the fracture. A permeability model considering both fracture topography and loading stress was developed. A simplified Hertz contact model was used to determine the relationship between average embedded depth and loading stress. A new fracture model that considers the evolution of fracture contact and void space under changing stress conditions was proposed. Based on the fracture model, a mathematical friction model was formulated to describe the full regime of fluid flow in rock fractures. The performance of the friction model was evaluated against several existing models [10, 12, 18] using experimental data and parametric analysis. A nonlinear flow model considering the fracture geometry evolution-based permeability model, mechanical model, and friction model was proposed.
2. Experimental Procedure
2.1. Sample Preparation and Surface Topography
A cube of granite rock was collected from an outcrop in Hubei Province, China. Cylindrical granite rock samples were cored using a 50 mm diameter drill and then cut into segments nearly 100 mm in length. Samples were ground and polished using a lapping machine. In order to eliminate the effect of water saturation, all the samples were baked for 2 h at a temperature of 105°C in a resistance furnace, retaining no water saturation [38]. The samples were cooled naturally to room temperature (20°C). After that, a splitting tool [13, 33], based on the theory of the Brazilian test, was used to create fractures (Figure 1(a)). Figure 1(b) shows the sample containing a single fracture. Basic geometrical information about the prepared samples is listed in Table 1.
The rough surfaces of the single fractured granite were characterized by a noncontact 3D laser profile scanner before the hydromechanical tests. The scanned data were input into the open-source software ParaView 5.6.0 (Kitware Inc., USA), a program designed for multiplatform data analysis and visualization. Figure 2(a) displays the 3D topography of both fracture sides for each sample. At first glance, it appears that the two sides of each fracture can be sealed tightly together. However, detailed assessment (Figure 2(b)) shows that the topographies of the two sides do not fit perfectly. Contact areas provide mechanical support for fracture under stress, while noncontact areas provide void spaces for fluid flow. The value of was used to evaluate the roughness of a single fracture surface [33]. It is defined as the ratio of the true length of a fracture surface trace to its projected length in the fracture plane. can be calculated by equation (1a) [33], and the JRC can be calculated by equation (1b) [39]. The results are listed in Table 1. where and represent the coordinate of the profile line in the length direction (m), and represent the coordinate of the profile line in the height direction (m), and is the length of the fracture surface parallel to the flow direction (m).
(a)
(b)
Ten evenly spaced profiles parallel to the flow direction were used to evaluate surface roughness. The resulting values are listed in Table 1. Note that the values are not identical for both sides of each sample. This is consistent with the results shown in Figure 2(b). The fracture surface of the sample TG3 is rougher than T3, which comparatively has a smooth surface.
2.2. Gas Flow Test
The gas flow test was conducted using a traditional triaxial servo test system, which is schematically demonstrated in Figure 3. The system consists of a main support frame, a loading system, a triaxial cell, gas cylinders, a pressure and flow rate monitoring unit, and a computerized data acquisition unit. The cylindrical sample is placed on the bottom seat of the triaxial cell. Then, the sample is covered with a polyurethane membrane. In addition, two clamps with O-rings are placed on the top and bottom pedestals over the membrane to avoid leakage. Subsequently, the cell is closed and filled with hydraulic oil. The hydraulic oil exerts confining stress on the sample, and the nitrogen gas cylinder provides inlet pressure. The flow outlet is kept open, meaning that the outlet pressure is approximately equal to the atmosphere pressure. A gas flow meter connected at the flow outlet is used to measure the flow rate.
Figure 4 shows the temporal development of the test conditions. During the test, different combinations of inlet pressures and confining stresses were applied to the sample. To obtain a wide range of volumetric flow rates, the inlet pressure varied from 0.005 MPa to 3.6 MPa. Note that the inlet pressure was always maintained at a lower level than the confining pressure to prevent sealing failure and sample collapse. Effective stress theory [40] was used to calculate the total confining stress (equation (2a)), assuming a Biot coefficient of 1. The pressure in the fracture was assumed to be the average of the inlet and the outlet pressures. Effective stress can be calculated from equation (2a) using the fracture pressure represented by equation (2b): where is the effective confining stress (MPa), is the total confining stress (MPa), is the Biot coefficient (-), is the pressure in the fracture (MPa), is the inlet gas pressure (MPa), and is the outlet gas pressure (MPa), which is equal to atmospheric pressure in this study.
3. Experimental Results and Discussion
3.1. Aperture Distribution
The 3D fracture topography was determined using coordinate data. To generate the aperture between two matched fracture surfaces, a geometrical model matching the size of the fracture was created. Then, the model was meshed into discrete elements at a spatial density based on the scan laser’s accuracy in the and directions. Each element has an aperture estimated from the gap between the two fracture surfaces. The initial state of the fracture (zero loading stress) is assumed to occur at the gap value at which 0.5% of the elements are in contact.
Figure 5 shows the aperture distribution for all the samples. In the present study, a Gaussian distribution was used to characterize the mean aperture size and the standard deviation of the aperture distribution. Table 2 shows the curve fitting results of the fracture surface. From Figure 5, it can be clearly seen that samples T3 and T8 are represented by tall, narrow curves, while sample TG3 is represented by a broader curve with a lower peak. Both T3 and T8 peak at approximately 6.5%, higher than the peak frequency of TG3; this indicates that T3 and T8 have more concentrated aperture distributions than sample TG3. However, the mean aperture size of sample T3, located approximately 0.4 millimeters below the value with the peak frequency, is larger than the mean aperture size of sample T8 (slightly larger than 0.25 millimeters). These results are consistent with the fitting results presented in Table 2.
3.2. Fracture Permeability
During the test, the volumetric flow rate at the outlet was directly measured. Figure 6 shows the relationship between the pressure gradient and the volumetric flow rate under different effective confining stresses. As the effective confining stress increases, the slope of the curves becomes lower and lower. This effect is probably caused by progressive fracture closure under increasing confining stress.
Several empirical equations have been proposed for nonlinear flow phenomenon in both porous and fractured mediums. One of the most widely used is the Forchheimer equation [41], a zero-intercept quadratic equation that can effectively characterize nonlinear flow processes at the macroscopic level without focusing on the transition from linear to nonlinear flow: where is the linear coefficient (kg·s-1·m-5) and is the nonlinear coefficients (kg·m-8); both are related to fluid properties and medium geometries [42]. is the pressure gradient (MPa/m), and is the volumetric flow rates (m3/s).
Since gas is strongly compressible, the flow rate within the fracture is different from the rate measured at the outlet and the Forchheimer equation is inadequate to be directly used to fit the experimental data. However, the Forchheimer equation can be written as equation (4). Based on the mass conservation theorem and the ideal gas law under isothermal condition, equations (5) and (6) can be obtained, respectively. Submitting these two equations into equation (4), one can obtain equation (7). This equation can be integrated to produce a quantitative expression (equation (8)) that can be used to fit the experimental data in the present study. The parameters (equation (9)) and (equation (10)) [43] can be obtained through the origin and join experimental data points plotted in vs. and in a three-dimensional space. where is the fluid density (kg/m3); , , and are the fluid density (kg/m3), volumetric flow rate (m3/s), and mass flow rate (kg/s) at the outlet, respectively; is the width of the fracture surface perpendicular to the flow direction (m); is the equivalent perpendicular gap of smooth parallel plates (m); is the fluid viscosity (Pa·s); is the non-Darcy coefficient (1/m); and is a nonlinear coefficient without considering fluid density (1/m5).
The fitted parameters and and the correlation coefficient are listed in Table 3. These results show that the correlation is strong—most of the fitting coefficient values are greater than 0.998. Figure 7 plots the variation of coefficients and with increasing effective confining stress. Both and increase by 4 orders of magnitude as the effective confining stress increases from 1 to 40 MPa. Coefficient , which represents the intrinsic permeability, is related to the viscous forces (viscous friction) at the fluid-solid interface [43]. Based on equation (9), the hydraulic aperture and the intrinsic permeability () can be calculated [43].
(a)
(b)
Table 4 lists the experimental permeability values calculated from coefficient , and Figure 8 shows the curves between intrinsic permeability and effective confining stress. With increasing stress, the intrinsic permeability declines rapidly at first. However, the decline quickly slows, indicating a decrease in the rate of fracture closure. The reason for this may be that the rate of increase in effective stress acting on the contact asperities becomes less as the contact area increases. Under high confining stresses (12 MPa), the permeability gradually approaches a constant value. By comparing different samples, one can easily conclude that sample T8 has the lowest intrinsic permeability, whereas sample T3 exhibits the highest intrinsic permeability. This is closely related to differences between the aperture distributions. As shown in Section 3.1, sample T8 has a concentrated aperture distribution and a small mean aperture size. This means that the void space capable of acting as flow channels is much smaller. In contrast, the high mean aperture size of T3 is associated with one of the highest intrinsic permeability values. This also implies that the permeability is positively proportional to the mean aperture size. Based on these analyses, a permeability model relating to confining stress and mean aperture size is proposed. It can be written as equation (11). Figure 9(a) gives the comparison of experimental data and predicted surface, and Table 5 (column 1) shows the fitting results. A good agreement is obtained between experimental data and the proposed model. where is the mean aperture size in the Gaussian interpolation distribution based on fracture aperture (m) and , , and are the fitting parameters.
(a)
(b)
The fitting parameter reflects the non-Darcy coefficient , which represents the nonlinear effect. When the nonlinear effect becomes negligible, that is , equation (8) reduces to Darcy’s law. As confining stress increases, both the relative roughness and tortuosity of the flow path also increase as a result of decreasing hydraulic aperture and increasing fracture contact area, respectively. This causes the fluid pressure to dissipate during the inertial process, resulting in an increase in the value of , as shown in Figure 7(b) [43].
Figure 2 and Table 1 show that TG3 is the roughest surface and T3 the smoothest. This can be seen by examining the value, which is related to the actual length of the profile line for a single fracture surface of each sample. Although TG3 is the roughest, sample T8 exhibits the lowest intrinsic permeability (Figure 8). This result is not consistent with the conclusion drawn from the analysis of the single fracture surface topography. This is due to the fact that flow channel characteristics are determined by the fracture contact state and void space structure, which depends on the characteristics of both fracture surfaces. It is clear that the hydraulic characteristics of rock fractures are largely determined by the distribution of fracture apertures [44].
4. Analysis of the Nonlinear Friction Factor
4.1. Nonlinear Friction Factor Based on the Planar Fracture Model
In the present study, because of the strong compressibility of nitrogen, the well-known Darcy-Weisbach equation [45] (equation (12)) cannot be directly used to calculate the friction factor. However, when comparing equations (12) and (4), one can obtain equation (13). Simplifying this equation and submitting equation (5) into this equation, the form for calculating the friction factor for gas flow in rock fracture under isothermal condition can be written as equation (14). The Reynolds number, , which is widely used to predict the boundary between linear and nonlinear flow behavior, is defined as the ratio of inertial forces to viscous forces within a fluid. For fluid flow in fractures, it is written as equation (15). where is the friction factor (-), is the hydraulic diameter (m) (two times the gap between the two parallel planes), is the average flow velocity through a cross-section perpendicular to the flow direction (m/s), and is the Reynolds number (-).
Figure 10(a) shows the friction factor as a function of the Reynolds number for the smooth parallel-plate model. It can be observed that the onset of nonlinearities happens at the Reynolds number around 10, and then the curves separate from the line calculated by the ideal cubic flow model and gradually turn to be horizontal with the increase of the Reynolds number. These results are similar to those obtained from flow tests in smooth and rough pipes [46]. The linear, transitional, and turbulent flow regimes are clearly seen from Figure 10(a). We can also observe that the friction factor in the transitional and turbulent flow regimes is randomly related to the increasing loading stress. This may be because the friction factor based on the planar fracture model ignores the evolution of fracture geometry. With the gradual increase of loading stress, the contact area goes up and the void space decreases, leading the flow path to become more tortuous. In other words, the flow resistance, described as the friction factor, should increase monotonously, rather than being randomly related to the loading stress as observed in Figure 10(a). So the evolution of fracture geometry during loading should be taken into account when calculating the friction factor.
(a)
(b)
4.2. Evolution of the Fracture Geometry during Loading
The deformation of the fracture was not possible to measure when using a cylindrical sample. In order to determine the evolution of the fracture geometry, a numerical simulation was conducted considering the mesofracture geometry (aperture distribution, Figure 11(a)). Equation (16) is the governing equation of the numerical model, describing a steady flow in a planer fracture. It was numerically formulated by use of the finite volume method (equation (17)). Fixed pressure as boundary conditions was applied at the left (inlet) and right (outlet) sides of the fracture (Figure 11(a)). At the outlet, the flow rate was calculated. Therefore, an average fracture permeability can be estimated through the numerically applied pressure gradient and the obtained flow rate. The local aperture changed when two fracture surfaces moved under loading. When the numerically estimated fracture permeability is equal to that measured from the experiment at a loading stage, then the relative displacement of the two fractures was considered as the deformation of the fracture at this loading stage. Table 4 compares the experimentally and numerically estimated permeability; the maximum relative error is less than 5‰. This allowed us to numerically obtain a mechanical displacement and a contact state under each specific confining stress in accordance with the experiment. Once the mechanical displacement was obtained, the contact ratio, void space volume, and embedded volume were calculated via area and volume integration of the contact and noncontact elements. where indicates the center element, indicates neighbor element of the center element, and is the connection area between the center and the neighbor element (m2).
Figure 11 shows the values for the aperture and the pressure distribution of sample T3 under effective confining stresses of 1, 8, 20, and 40 MPa. The noncontact area decreases obviously as the confining stress increases from 1 to 20 MPa. In contrast, it decreases slowly as the confining stress increases from 20 to 40 MPa. Based on the pressure distribution (Figure 11(b)), it is obvious that the flow path is tortuous rather than straight. The nonlinearity of the pressure distribution also strengthens as the stress increases. This is likely due to the associated irregular increase of the contact area, which leads to the development of more complex and tortuous flow channels.
The two factures gradually come into contact under increasing stress. Since the actual contact area supports all of the effective stress acting on the fracture, it has a significant influence on the mechanical displacement of the fracture. The void space, which determines the size of the flow channel, can have an important effect on flow characteristics within fractures. The embedded volume has an important influence on fracture stiffness in the theoretical Hertz model [47]. The Hertz model assumes that the asperities would become imbedded in the fracture surface, allowing displacement to increase. Based on this assumption, the embedded volume is considered to be the total volume of the asperities embedded into each side of the fracture.
Figure 12 is a schematic diagram of aperture distribution. When the displacement equals (), all elements with an aperture size smaller than come into contact, while all elements with an aperture size larger than remain open. In the diagram, the contact area is represented by the shape (equation (18)). The void space can be calculated by taking the static moment of shape to the axial (equation (19)). The embedded volume can be calculated by taking the static moment of shape to the axial (equation (20)). Finally, the contact ratio and average embedded depth can be calculated using the contact area, embedded volume, and normal projection area of the fracture (equations (21) and (22)). where is the contact area (m2), is the displacement under certain loading stress (m), is the normal projection area of an element (m2), is the aperture size of an element (m), is the void space (m3), is the maximum aperture size (m), is the embedded volume (m3), is the contact ratio (-), is the normal projection area of the fracture (m2), and is the average embedded depth (m).
Figure 13(a) plots mechanical displacement as a function of loading stress. Rapid increases in normal displacement can be easily noticed at low stresses. However, once the stress reaches approximately 12 MPa, the rate of increase becomes much slower. This behavior is consistent with many other experimental studies of nonlinear deformation mechanisms of rock joints under normal stress [48–50]. It is worth noting that sample T8 has a much smaller displacement than samples T3 and TG3. The aperture distribution in Figure 5 shows that the proportion of small-sized apertures for sample T8 is much greater than that for the two other samples. Since the mean aperture size for T8 is the smallest, its displacement is the smallest under each loading stress. A similar relationship between displacement and the proportion of small-sized apertures exists for samples T3 and TG3. This leads us to conclude that the proportion of small-sized apertures influences the normal displacement significantly and that the normal displacement under loading stress is positively correlated with the mean aperture in this study. Figure 13(b) shows the average embedded depth versus loading stress. With an increase in loading stress, the average embedded depth increases quickly under small loading stresses and slowly under large loading stresses. This increasing pattern is similar to the pattern observed for fracture normal displacement.
(a)
(b)
Figure 14(a) shows the contact ratios of different samples as a function of the effective confining stress. As the confining stress increases from 1 MPa, the contact area increases rapidly. However, as the stress continues to increase, the rate of increase slows due to the increase in the contact area. Figure 14(b) demonstrates the void space as a function of the effective confining stress, which can be used to evaluate the flow ability of the fracture to some extent. This decreasing trend shows a high rate of decrease at low confining stresses followed by a low rate at high confining stresses. The increase in contact area causes the flow path to become more torturous, and the reduction of void space leads to narrower flow channels. These two factors jointly lead to the observed decrease in intrinsic permeability.
(a)
(b)
4.3. Nonlinear Friction Factor Based on the Planar Fracture Model considering Evolution of the Fracture Geometry
Since the contact area affects the tortuosity of the flow path and the void space determines the size of the space the fluids flow through, it is important to take both factors into account when calculating the friction factor. In this study, a modified smooth parallel-plate model for flow friction in fractures was proposed (Figure 15). In this model, discrete contact areas in a fracture are merged to form mechanical supports on either side of an open channel [51]. The channel allows fluid flow and has the same volume as the total void space of the fracture.
(a)
(b)
In the smooth parallel-plate model, shown in Figure 15(a), the hydraulic aperture under a certain confining stress is calculated using equation (9). In the new model, effective fracture width , geometrical aperture , and hydraulic aperture are related to the contact ratio and the total void space. They are calculated using equations (23)–(25). The friction factor and Reynolds number in the new model can be calculated using equations (26) and (27). where is the width of the fracture surface perpendicular to the flow direction in the new model (m), is the geometrical aperture (m), and is the hydraulic aperture of the new model (m), shown in Figure 15(b).
The Re-f curves derived from the new physical model are markedly influenced by the effective confining stresses (Figure 10(b)). Specifically, the friction factor calculated based on the new physical model considering fracture geometry evolution increases monotonously with the loading stress. This is more theoretically logical because increase in confining stress causes increasing fracture closure, leading to a greater contact area and a reduction in the void space, resulting in a more complex flow channel and inducing extra pressure loss.
5. Nonlinear Flow Model considering the Fracture Geometry and Loading
5.1. Nonlinear Friction Model
In previous studies [12, 18], a dimensionless parameter called the relative roughness coefficient was used to describe the complexity of the flow channel. It is defined as the ratio of the initial maximum asperity to the hydraulic aperture. However, the initial maximum asperity is a constant and cannot reflect the evolution of the flow channel. Thus, a coefficient that reflects the dynamic evolution of relative roughness should be developed. In the present study, this problem was solved by including the geometrical aperture in the calculation of the relative roughness coefficient. The relative roughness coefficient is expressed as the ratio of the absolute roughness, assumed to be (Figure 15(b)), and the geometrical aperture of the flow channel : where is the relative roughness coefficient (-).
Figure 16 shows all of the friction data as a function of the Reynolds number and the relative roughness. All the experimental data are distributed over a surface, allowing us to see the influence of relative roughness on the friction factor. Overall, the friction factor increases exponentially with increases in the relative roughness. T8 has the overall roughest flow channel, whereas the flow channel of T3 is relatively smooth. These results differ from, and are more realistic than, those obtained using the value alone (Table 1) because an value calculated based on a single fracture surface topography cannot describe the tortuosity of a contact rock fracture. In addition, is not adequate to describe the dynamic evolution of fracture roughness under changing stress.
From the above analysis, it can be concluded that flow friction in a rock fracture is related to two factors: the Reynolds number, which determines the flow regime, and the fracture void geometry, characterized by fracture roughness. Despite the progress achieved by previous researchers such as Nazridoust et al. [10] (equation (29)), Zhang and Nemcik [12] (equation (30)), and Zhou et al. [18] (equation (31)), the relationships between friction factor and flow regime and fracture roughness for fluid flow in rock fractures remain unclear. where is the peak asperity height (m).
In the present study, both nonlinear effects and the influence of relative roughness are considered by multiplying the ideal parallel-plate model by two power-law functions (equation (32)). When the Reynolds number is small and the fracture is absolutely smooth, equation (32) becomes the ideal friction model (). When substituting equations (27) and (28) into equation (32), the friction factor model intuitively reflects the effects of fracture contact and void space evolution can be obtained (equation (33)): where , , , and are the fitting parameters.
The proposed model was used to fit the experimental data. Table 5 (column 2) presents the fitting results. The correlation coefficient is approximately 0.973, signifying a high degree of correlation. Figure 17 compares the experimental and the predicted data. Both show good agreement, implying that the proposed model is reasonable for the prediction of flow friction in a rock fracture.
The performance of the proposed model was investigated by using a parametric sensitivity study and comparing the results with Nazridoust et al.’s model [10], Zhang and Nemcik’s model [12], and Zhou et al.’s model [18] (Figure 18). Nazridoust et al.’s model employs a noncontact fracture and ignores fracture roughness. Because of this, the model produces a single line and fails to describe the influence of fracture roughness. Zhang and Nemcik’s model validates the idea that for real and contact rock fractures, the friction factor is larger than that predicted by the ideal parallel-plate model, even for very low Reynolds numbers (). The results are represented by a set of straight lines parallel to the ideal friction model. However, the model is not able to characterize the nonlinear friction effect (Figure 18(a)) because the Reynolds number is quite small (). Zhou et al.’s model considers both fracture roughness and nonlinear friction effects, but it tends to underestimate the friction factor for low Reynolds numbers (), where all the Re-f curves lie together, close to the ideal friction model. All of the above results differ from those obtained using the new model (Figure 18(b)). At low Reynolds numbers, the lines are parallel to those representing the ideal friction model. However, as the Reynolds number increases, the lines turn gradually to the horizontal direction, representing the gradual transition to a turbulent flow regime observed experimentally [15, 33].
(a) Nazridoust et al.’s, Zhang and Nemcik’s, and Zhou et al.’s models
(b) Proposed model
5.2. Mechanical Model
In this study, the Hertz contact theory [47] was used to explain the deformation mechanism of the fracture. As normal displacement increases, more asperities become imbedded into fracture surfaces, resulting in a gradual increase in the contact area. This leads to higher fracture stiffness and a slower increase in the displacement rate under large loading stresses. The Hertz contact theory predicts that when the interaction between asperities is ignored, the relationship between normal loading stress and deformation for a single asperity can be calculated using equation (34). Here, the Hertz contact model was used to describe the average fracture embedded depth and normal loading stress. To implement the model, a variable that can indicate the average curvature radius of the fracture should be determined. The standard deviation of aperture distribution can indicate the closeness of the spatial fit between two rough surfaces of one fracture. A low standard deviation indicates that apertures tend to be close to the mean aperture size, which implies a large curvature radius and small average embedded depth. In contrast, a high standard deviation indicates that aperture sizes are spread out over a wider range of values, reflecting a smaller curvature radius. Based on this relationship, the reciprocal of the standard deviation was used to characterize the curvature radius of the fracture. By integrating the average embedded depth and reciprocal of the standard deviation into the Hertz contact model, a simplified model (equation (35)) was obtained. Table 5 (column 3) shows the fitting results, and Figure 9(b) demonstrates the fitting surface of the normal loading stress and compares it with experimental data. Although the average embedded depth data were obtained numerically and thus have associated measurement errors, a good agreement between experimental data and the model was achieved (correlation coefficient ). The fitting results prove that the numerical simulation based on the flow result is valid and that the experimental data are reasonable. where is the effective elastic modulus (MPa), is the average curvature radius (m), is the asperity deformation (m), is the standard deviation of aperture (m), and , , and are the fitting parameters.
5.3. Effect of Fracture Geometry Parameters on Nonlinear Flow
Based on analysis, a nonlinear flow model derived from the cubic flow model was proposed in the present study. When taking into consideration the mechanical model (equation (36a)), the permeability model (equation (36b)), the friction model (equation (36c)), and the fracture geometry parameters, the cubic flow model can be modified to have equation (36d). Once the 3D fracture topography data was obtained, the mean aperture size and standard deviation were calculated through the Gaussian distribution model fitting. Assuming that a certain stress was loaded on the fracture, the average embedded depth under the given loading stress was obtained using the mechanical model (equation (36a)). In addition, the permeability was calculated by use of the mean aperture size and the given stress (equation (36b)). Once the average embedded depth was obtained, the displacement was computed using equation (20). Afterwards, the contact ratio and the void space were determined through equations (18) and (19), respectively. Then, both the void space and the contact ratio served to determine the geometrical aperture while the permeability and the contact ratio were used to assess the hydraulic aperture in the new physical model. Once these two parameters were known, the relative roughness was estimated through equation (28). Therefore, the nonlinear friction factor under a given Reynolds number could be obtained using the newly built friction model (equation (36c)). In the end, all the parameters (, , and ) used in the nonlinear flow model (equation (36d)) for predicting the flow rate are obtained. Consequently, the 3D fracture topography data, the given loading stress, and the applied pressure gradient can be used to infer the nonlinear flow rate.
Figure 19 shows the intrinsic permeability under different mean aperture sizes as a function of the loading stress. It can be clearly seen that the permeability first decreases sharply during the low-loading stress stage (less than 10 MPa) and is then followed by a slow and steady decreasing trend. High mean aperture size results in high permeability, which is consistent with the experimental results in Section 3.2. Figure 20 investigates the effect of fracture geometry parameters and on the relative roughness . Figure 20(a) shows that the relative roughness would drop with the increase of when the remains constant. This is because higher means both a higher (equation (36a)) and a higher displacement (equation (20)), thus leading to a lower and a higher . According to equation (28), it is possible to obtain a lower . Regarding the effect of on the as shown in Figure 20(b), since a higher means a higher in the new physical model, the obtained through equation (28) would be lower when both the loading stress and remain unchanged. Figure 21 presents the friction factor as a function of the loading stress ranging from 2 MPa to 100 MPa. It is similar to the experimental results based on the proposed model in the present study in Figure 10(b). As explained in Section 4.3, a higher loading stress causes a higher fracture closure, which leads to a greater contact area and a reduction in the void space, resulting in a more complex flow channel and an extra pressure loss.
(a)
(b)
From the analysis in Section 3.1, one can find out that sample T8 has a low mean aperture size and a low standard deviation and T3 has a high and a low , while TG3 has a high and a high . Figure 22 shows the relationship between the friction factor and the fracture geometry parameters and under loading stresses of 20 MPa and 100 MPa. In this case, four types of aperture distribution were described based on experimental information regarding the granite fracture aperture distribution. Type I has a relatively low (0.2 mm) and a low (0.05 mm), which refers to the distribution characteristics of sample T8. Type II has a relatively low (0.2 mm) and a high (0.2 mm). Type III has a relatively high (0.4 mm) and a low (0.05 mm), which refers to the distribution characteristics of sample T3. Type IV has a relatively high (0.4 mm) and a high (0.2 mm), which corresponds to the distribution characteristics of sample TG3. From Figure 22(a), it can be noticed that type I has the highest friction factor, which means that more pressure loss happens when the fluid flows through the fracture. This conclusion is consistent with the one obtained from the experimental results in Section 5.1 supporting that T8 has the roughest flow channel and the highest friction factor under the same loading stress. The friction factor of type III and IV aperture distribution is obviously lower than that of types I and II, which implies that the friction factor is sensitive to the mean aperture size . Figure 22(b) shows that with the increase of the loading stress from 20 to 100 MPa, the effect of the standard deviation on the friction factor weakens considerably, which means that under a high-loading stress, the effect of the standard deviation on the friction factor diminishes. However, the influence of the mean aperture size on the friction factor is still obvious even when the loading stress is high.
(a)
(b)
6. Conclusions
In this study, the effect of the fracture geometry, which determines aperture distribution, fracture contact, void space, and embedded depth, on the nonlinear flow and the mechanical behavior of a single-fractured granite rock was examined by performing laboratory flow tests with effective confining stresses ranging from 1 to 40 MPa and inlet pressures ranging from 0.005 to 3.6 MPa. The conclusions are summarized here. (1)The aperture distribution can be used to explain the permeability. The intrinsic permeability is positively correlated with the mean aperture size. T8 has the smallest mean aperture size and the lowest intrinsic permeability, while T3 has the highest values. This is not consistent with conclusions drawn from the traditional surface roughness based on the topography of a single fracture surface (), which indicates that TG3 is the roughest and T3 the smoothest fracture(2)The proportion of small-sized apertures influences the normal displacement significantly, and the higher the mean aperture size, the larger the normal displacement for a given stress. The average embedded depth increases rapidly under small loading stresses and slowly under large loading stresses. A simplified mechanical model based on the Hertz contact theory was used to fit the embedded depth data. The model uses the reciprocal of the standard deviation of aperture to characterize the average curvature radius of the fracture. Good agreement (with a correlation coefficient of 0.935) was obtained. The fitting result demonstrates that the numerical simulation based on the flow result is valid and the experimental data are reasonable(3)The linear, transitional, and turbulent flow regimes are clearly seen from the figure. The ideal cubic flow model underestimates the friction factor. In the new physical model, which considers both fracture contact and void space, the friction factor data points are located above the ideal friction line. This implies that our method for calculating the friction factor and Reynolds number is reasonable. Based on the experimental results, a new friction model was proposed, in which both nonlinear effects and the influence of fracture roughness are considered by multiplying two power-law functions additionally. A comparison of the proposed and existing models using parametric analysis shows that the new friction factor model can more effectively characterize total flow behavior in a rock fracture(4)By integrating the built permeability model, the friction model, and the mechanical model into the cubic flow model, a nonlinear flow model was proposed. The flow rate was estimated based on the 3D fracture topography data, the given loading stress, and the applied pressure gradient. When the mean aperture size remains constant, the relative roughness decreases with the increase of the standard deviation, and when the standard deviation is unchanged, the relative roughness decreases with the increase of the mean aperture size. The investigation of influence of four types of aperture distribution on the friction factor showed that a high mean aperture size would result in a low friction factor. In addition, the influence of the mean aperture size on the friction factor is obvious, even under a very high-loading stress while the effect of the standard deviation is negligible when the loading stress is very high
Data Availability
All the original data (data points in Figure 6) and analysed data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The work presented in this paper is funded by the National Natural Science Foundation of China (No. 51774056), the Chongqing Research Program of Basic Research and Frontier Technology (No. cstc2017jcyjB0252), the National Science Fund for Distinguished Young Scholars of China (No. 51625401), and the Chongqing University Postgraduates’ Innovation Project (CYB17047). We would like to thank LetPub (http://www.LetPub.com) for providing linguistic assistance during the preparation of this manuscript.