Effects of Contact Area and Contact Shape on Nonlinear Fluid Flow Properties of Fractures by Solving Navier-Stokes Equations

The in ﬂ uences of contact shape and contact area on nonlinear ﬂ uid ﬂ ow properties through fractures are investigated by solving Navier-Stokes equations. The evolutions of nonlinear relationships between ﬂ ow rate and hydraulic pressure drop, Forchheimer coe ﬃ cients, nonlinear factor, critical hydraulic gradient, distributions of ﬂ ow streamlines, and tracer ﬂ ow paths at di ﬀ erent times are systematically estimated. The results show that the nonlinear relationships between ﬂ ow rate and hydraulic pressure drop can be well described by Forchheimer ’ s law, in which the nonlinear term coe ﬃ cient b is approximately three orders of magnitude larger than the linear term coe ﬃ cient a . The smaller contact area corresponds to smaller variations in many aspects such as ﬂ ow rate, critical hydraulic gradient, ﬂ ow streamlines, and tracer ﬂ ow paths. The critical hydraulic gradient decreases with the increasing degree of contact shape variations while the contacts have the same mean area. The increase in hydraulic pressure drop can induce signi ﬁ cant eddies and decrease the permeability and/or conductivity of fractures. However, the distributions of streamlines and tracer ﬂ ow paths are not dramatically disturbed under a large hydraulic pressure drop.


Introduction
Accurate estimations of hydraulic properties of rock fractures are of very importance in many fields such as geothermal energy development [1,2] and enhanced oil recovery [3][4][5]. It is well known that both the contact area ratio of fractures and high hydraulic pressure drop can significantly influence the permeability and/or conductivity of rock masses (Chen et al. 2015; [6,7]). However, the effects of contact shape and contact area on nonlinear flow properties of fluids through fractures induced by high hydraulic pressure drop by solving Navier-Stokes equations have not yet been systematically investigated.
In early works, the relationships between hydraulic aperture and mechanical aperture have been studied, and mathematical expressions that incorporated the contact area ratio have been proposed. Walsh [8] concluded that e 3 /E 3 = ð1 − CÞ/ð1 + CÞ, where e is the hydraulic aperture, E is the mechanical aperture, and C is the contact area ratiothat is defined as the ratio of contact area to the area of surfaces. Zimmerman et al. [9] reported that C = 0:25 in their studies. Zimmerman and Bodvarsson [10] derived that e 3 /E 3 ≈ ð1 − 1:5σ 2 apert /E 2 +⋯Þð1 − 2CÞ. However, the quantifications of the contact area ratio are difficult because the rock mass is invisible. Wang et al. [6] carried out shearflow tests on rough fractures and calculated the contact area ratio according to the shear displacement and normal displacement using a code developed by themselves. They found that as shear displacement increases from 1 to 3 mm, the contact area ratio dramatically decreases from 0.33 to 0.09. When the shear displacement continuously increases to 9 mm, the contact area ratio fluctuates in a negligibly small range. Liu et al. [11] performed numerical simulations to model fluid flow through rough-walled fractures during shearing by solving Navier-Stokes equations. The results show that the contact area ratio increases linearly with boundary stiffness and is in the range of 0.01~0. 16 when boundary stiffness = 0~2 GPa/m, joint compression strength = 100~150 MPa, and initial normal stress = 0:52 :0 MPa. However, the contact area and locations are dependent on shear displacement and fracture surface morphologies [12], and the effect of the contact area ratio on nonlinear fluid flow properties is not estimated.
The previous studies commonly assumed that fluid flow is in the linear flow regime, in which the flow rate is linearly correlated with hydraulic pressure drop and the cubic law is 2 Lithosphere typically used for modeling fluid flow through fractures [13][14][15][16][17]. However, in the karst systems and/or in the pumping tests [18,19], the flow rate is very large and the inertial force cannot be negligible, resulting in a nonlinear relationship between flow rate and hydraulic gradient. Thus, the permeability and/or conductivity of fractures are reduced [20][21][22]. Although the effects of joint surface roughness (JRC), Reynolds number (Re), and shear displacement on nonlinear properties of fluid flow have been systematically investigated [23][24][25][26], the influences of the contact shape and contact area have not been studied in an in-depth way. The present study, first, numerically generated a series of fracture models containing contacts with different shapes and different area ratios. Then, the Navier-Stokes equations are solved based on the finite volume method (FVM), and the nonlinear relationships between the flow rate and hydraulic drop, as well as linear and nonlinear coefficients in Forchheimer's law and critical hydraulic gradient, are quantitatively analyzed. Finally, the streamlines of fluid flow under different shapes and areas of contacts and different hydraulic pressure drops are estimated, and the distributions of tracer flow paths at different times are discussed.

Governing Equations of Fluid Flow through Fractures
The fluid flow through fractures is governed by the Navier-Stokes equations, by assuming that the fluid is incompressible and is a kind of Newtonian fluids. The tensor form of Navier-Stokes equations is written as shown in [27,28] ρ ∂u ∂t where u is the flow velocity tensor, ρ is the fluid density, p is the hydraulic pressure, T is the shear stress tensor, t is the time, and f is the body force tensor. The flow rate Q can be calculated by multiplying the cross-sectional area A by u. The units of Q, A, and u are m 3 /s, m 2 , and m/s, respectively. It is difficult and time-consuming to solve the Navier-Stokes equations, because a series of partial differential equations should be solved by iterations. Therefore, some useful empirical functions are proposed such as Forchheimer's law [22,27,29,30], which can be written as where ΔL is the length of the model and a and b are the linear and nonlinear coefficients, respectively. When Q is sufficiently small, the values of bQ 2 are much smaller than aQ. Thus, equation (2) can be reduced to Equation (3) is similar to Darcy's law, in which the flow rate is linearly correlated with hydraulic pressure drop. For fluid flow through fractures, equation (3) can be expressed as the following cubic law [31,32]: where w is the width of fractures, e is the hydraulic aperture, and μ is the dynamic viscosity. When bQ 2 is not sufficiently small with respect to aQ, the terms bQ 2 cannot be negligible, resulting in a nonlinear relationship between Q and ΔP/ΔL. A nonlinear factor E is utilized to quantify the nonlinearity of fluid flow, defined as E represents the percentage of hydraulic pressure drop induced by inertial forces to the total hydraulic pressure drop. Although E = 0:01 and E = 0:5 have been used for quantifying the onset of nonlinear flow of fluids in some previous works [33,34], it is commonly adopted that E = 0:1 is the threshold [20,21]. In the present study, E = 0:1 is chosen as the transition of fluids from a linear regime to a nonlinear regime. The hydraulic gradient (J) that corresponds to E = 0:1 is the critical hydraulic gradient (J c ). J is a nondimensional parameter, which is defined as the ratio of hydraulic head difference to the length of fractures. Note that the Re can also characterize the fluid flow properties of single fractures yet cannot interpret the flow state of complex fracture networks. However, J has a known value as long as the hydraulic water heads applied on the opposite boundaries are set. Therefore, J is suitable to characterize the flow state for both single fractures and fracture networks and is adopted and utilized.

Numerical Simulation Results and Analysis
3.1. Numerical Fracture Models Containing Contacts. Ten two-dimensional square numerical fracture models are    4 Lithosphere and 25%, respectively. The side length of each contact (L c ) is calculated as follows: where m and n are coefficients and R is the random number that ranges from 0 to 1. The values of L mean , m, and n are listed in Table 1. A parameter K is used to represent the deviation of side length of contacts from the mean values, which is calculated as follows: A larger K represents that the contacts are more scatteredly distributed, while a larger L mean indicates a larger contact area.
The fracture models as shown in Figure 1 are then imported into ANSYS ICEM for meshing with a maximum side length of 0.2 mm. The meshed files are imported into ANSYS FLUENT to model fluid flow. The hydraulic pressure difference (ΔP) is assigned to be 0.01, 0.05, 0.1, 0.3, and 0.5 Pa, respectively. When the fluid flow is in a steady state, the flow rate is calculated, and the relationships between the flow rate and hydraulic pressure drop can be quantitatively analyzed. The water at a room temperature is selected as the fluid, flowing from the left side of the model to the right side. The density is 998.2 kg/m 3 and the dynamic viscosity is 0.001003 Pa·s.

Nonlinear Relationships between Flow Rate and
Hydraulic Pressure Drop. Figure 2 shows that with the increment of Q, ΔP increases nonlinearly at an increasing rate. The nonlinear relationships between Q and ΔP follow the Forchheimer's law as shown in equation (2), in which the intercept of the quadratic functions is zero. As L mean increases, the contact area increases while the void spaces for fluid flow decrease, resulting in decreases in Q when ΔP is fixed. Taking L mean = 3 mm and K = 0 (Figure 2(e)) and L mean = 5 mm and K = 0 (Figure 2(j)) as the examples, when ΔP = 0:5 Pa, Q decreases from 0.00104 m 3 /s to 0.0004895 m 3 /s, by a rate of 52.93%. When L mean = 3 mm, Q increases in a limited range with decreasing K from 2 to 0, whereas when L mean = 5 mm, Q increases significantly (from 0.0002557 m 3 /s to 0.0004895 m 3 /s by a rate of 91.44%) with decreasing K from 4 to 0. This indicates that when the contact area represented by L mean is small, the change in contact shape distributions represented by K will not change robustly, because the small area of each contact would not influence the flow paths to a large extent. However, when L mean is large, the void spaces between contacts may be narrowed and the flow paths change significantly. As a result, energy losses and Q decreases obviously.
The variations in Forchheimer coefficients a and b of functions as shown in Figure 2 are presented in Figure 3 and Table 1. When L mean = 3 mm, a changes slightly in the range of 252~264 Pa·s·m -4 , whereas a changes significantly from 800 Pa·s·m -4 to 1320 Pa·s·m -4 when L mean = 5 mm. The variations in b are similar to those of a for L mean = 3 and 5 mm. Therefore, the larger L mean induces a more significant change in both a and b. The values of b are generally 2~3 orders of magnitude larger than those of a, which are Pa·s·m -4 to 841 Pa·s·m -4 by a rate of 5.13% and b increases from 4:52 × 10 5 Pa · s 2 · m −7 to 1:09 × 10 6 Pa · s 2 · m −7 by a rate of 141.15%. However, when K increases from 3 to 4, a increases from 841 Pa·s·m -4 to 1320 Pa·s·m -4 by a rate of  6 Lithosphere 56.60% and b increases from 1:09 × 10 6 Pa · s 2 · m −7 to 2:48 × 10 6 Pa · s 2 · m −7 by a rate of 127.52%. It indicates that both a and b for a relatively large L mean are more sensitive to a larger K, since the larger K can significantly influence the void space distributions and flow paths.

Determination of Critical Hydraulic Gradient.
The variations in E versus J are shown in Figure 4. With the increment of J, the inertial force increases following quadratic functions with a nonzero intercept, resulting in a stronger nonlinearity of fluid flow represented by a larger E. The E~J curves change to a smaller extent when L mean = 3 mm compared with those when L mean = 5 mm, which are similar to those of ΔP~Q curves, a~K curves, and b~K curves as presented in Figures 2 and 3. For L mean = 5 mm, E decreases from 0.324 to 0.217 when J = 0:0005, by a rate of 33.02%. The variations in J c versus K for L mean = 3 mm and 5 mm are plotted in Figure 5. With the increment of K, J c fluctuates slightly and then decreases significantly and finally fluctuates slightly again. J c varies from 3:36 × 10 −5 to 4:19 × 10 −5 for K = 0~2 and L mean = 3 mm, while J c varies from 8:59 × 10 −5 to 1:77 × 10 −4 for K = 0~4 and L mean = 5 mm. This indicates that J c is more sensitive to a larger L mean . As L mean increases, the mean aperture providing flow paths for fluids decreases, resulting in a larger J c . The same conclusions have also been reported by Liu et al. [33]. Figure 6 shows the streamline distributions for different L mean , K, and ΔP. Comparing Figure 6(a) with Figure 6(e), it clearly exhibits that the eddies are not formed at ΔP = 0:01 Pa with fluid smoothly flowing by passing the contacts yet are formed at ΔP = 0:5 Pa behind the contacts. These eddies induce energy losses and decrease the permeability and/or conductivity. Figures 6(a) and 6(b) show that although the flow paths are changed between contacts, the distributions of flow paths are not influenced robustly as a whole. However, the flow paths are significantly changed for a relatively large L mean such as L mean = 5 mm as shown in Figures 6(c)

Conclusions
The present work studied the effects of the contact shape and contact area on nonlinear flow properties of fluids through fractures by numerically solving Navier-Stokes equations. A total of ten numerical models with different contact shapes and shapes are generated and utilized for modeling fluid flow. The nonlinear relationships between the flow rate and hydraulic drop, as well as the variations in Forchheimer coefficients and critical hydraulic gradient, are characterized. Finally, the streamline distributions and tracer flow paths are presented and analyzed.

Lithosphere
The results show that with the increment of hydraulic pressure drop from 0.01 Pa to 0.5 Pa, the flow rate increases nonlinearly, following Forchheimer's law. Under a fixed hydraulic pressure drop, the flow rate does not change significantly versus K when the mean side length of contacts is small (i.e., 3 mm), yet the flow rate decreases dramatically versus K when the mean side length of contacts is large (i.e., 5 mm). Here, K represents the maximum deviation of side lengths of contacts from the mean value. As a result, the variations in Forchheimer coefficients a and b, nonlinear   The present study is aimed at investigating the effects of the contact area and contact shape on nonlinear flow properties of fluids through contacted fractures. Some assumptions are made that the contact is rectangular-shaped and located uniformly along both xand y-directions. However, these assumptions greatly deviate from natural cases where the contacts are formed by shearing or coupled thermo-hydromechanical-chemical (THMC) processes. Therefore, in the future works, we will focus on the nonlinear fluid flow properties of natural fractures that have rough surface and contacts in a natural state.   10 Lithosphere