Abstract

The knowledge of flow phenomena in fractured rocks is very important for groundwater-resources management in hydrogeological engineering. The most commonly used tool to approximate the non-Darcy behavior of the flow velocity is the well-known Forchheimer equation, deploying the “inertial” coefficient that can be estimated experimentally. Unfortunately, the factor of roughness is imperfectly considered in the literature. In order to do this, we designed and manufactured a seepage apparatus that can provide different roughness and aperture in the test; the rough fracture surface is established combining JRC and 3D printing technology. A series of hydraulic tests covering various flows were performed. Experimental data suggest that Forchheimer coefficients are to some extent affected by roughness and aperture. At last, favorable semiempirical Forchheimer equation which can consider fracture aperture and roughness was firstly derived. It is believed that such studies will be quite useful in identifying the limits of applicability of the well-known “cubic law,” in further improving theoretical/numerical models associated with fluid flow through a rough fracture.

1. Introduction

Understanding fluid flow behaviors in fractured rock aquifers is of great concern in numerous industrial and scientific fields, such as water resources management, sustainable urban drainage, contaminant pollution control, hazardous wastes isolation, and geothermal fluid or hydrocarbon exploitation.

It has been observed that fluid flow behavior through a rock mass depends on the geological-structural and geological-technical (such as degree of fracturing, orientation, persistence, weathering, moisture conditions and seepage aperture and filling, and roughness) [1]. As we all known, quantity and proximity of fractures pervading a rock mass define the degree of fracturing of the rock and, as a consequence, its permeability. Evidently, the higher the persistence of discontinuities, the stronger the permeability. The weathering process determines an increase of porosity. When it comes to a single fracture, aperture and roughness will be the dominant factors. An adequate knowledge of fluid flow through a rough-walled fracture is a starting point for a better interpretation of fluid flow and solute transport in fractured rock aquifers.

A wide number of groundwater flow models and numerical tools were developed based on the assumption of Darcy’s flow through individual fractures. It has been well recognized, however, that the linear Darcy’s law is not always adequate to describe the flow behaviors in natural fractures. Typical characteristics of natural rock fractures include rough walls and asperity contact [2, 3], and non-Darcy flow may occur as a result of nonnegligible inertial losses. Previous experimental work demonstrated that Darcy’s law fails to predict pressure drops in fractures when inertial effects are relevant before the fully developed turbulence [47]. In the post-Darcy regimes where inertial effects are significant, two equations are invariably used to describe pressure drop as a function of average velocity: Forchheimer and Ergun equation [7, 8]. Other seemingly different correlations may be traced back or manipulated to fit the basic forms of these two equations [9].

To account for such inertial losses, a widely used equation referred to as Forchheimer’s law was developed. Qian et al. [10] used well-controlled laboratory experiments to investigate the flow and transport in a fracture under non-Darcy flow conditions, and it is found that the Forchheimer equation fits the experimental - relationship nearly perfectly, whereas the Darcy equation is inadequate in this respect. Cherubini et al. [11] investigated the nonlinear flow by analyzing hydraulic tests on an artificially created fractured rock sample, and the experimental result shows that the relationship between the flow rate and head gradient matches the Forchheimer equation and describes a strong inertial regime well. Javadi et al. [12] performed both laminar and turbulent flow simulations for a wide range of flow rates in an artificial three-dimensional fracture. They developed a new geometrical model for nonlinear fluid flow through rough fractures, which suggested a polynomial expression, like the Forchheimer law, to describe the dependence of pressure drop on discharge. Considering that the mechanisms of non-Darcy two-phase flows in fractures are not well understood and no general model has been presented to describe them, relying on the full cubic law, Radilla et al. [13] carried out a series of flow experiments and presented a model to solve this preceding issue.

To incorporate the Forchheimer’s law in the analytical or numerical solutions for the simulation of non-Darcy flow, the determination of the phenomenological coefficients in Forchheimer equation becomes indispensable. For many practical problems, however, it may be difficult to directly determine the coefficients through experimental tests. In this circumstance, parametric expressions have to be used instead. For porous media, numerous theoretical and empirical expressions of coefficients have been developed, which are typically characterized as a function of the particle diameter and the medium porosity [14]. For fluid flow in rock fractures, however, the study on the parametric expressions for these two coefficients and their physical background is yet insufficient and rarely reported in the literature.

Roughness has a large influence on fluid flow and transport through tight, rough-walled fractures where non-Darcy flow is particularly easy to occur [15]. Qian et al. [10, 16] created the surface roughness by gluing small Plexiglas plates to one of the two vertical walls, keeping the second vertical wall smooth. Zhang and Nemcik [17] investigated the flow behavior in fractures, which were created in the laboratory by splitting initially intact samples of a fine grain sandstone block into halves.

More recently, the developed three-dimensional flow models were used to simulate fluid flow through various random synthetic rough-walled fractures. The rough-walled fractures were created by combining random fields of aperture and the mean wall topography or midsurface, which quantifies undulation about the fracture plane. Javadi et al. [12] performed both laminar and turbulent flow simulations for a wide range of flow rates in an artificial three-dimensional fracture.

The primary motivation of the present study was to experimentally evaluate the Forchheimer equation coefficients for non-Darcy flow in Forchheimer equation in rough-walled fractures, where effects of fracture roughness, fracture aperture, and flow regime were considered and effects of orientation and fracture networks are not, as the focus is not on bulk flow, but rather on the properties of individual discrete fractures. A series of hydraulic tests were performed under different roughness and fractures apertures, with the fracture formed by two planes and the flow taking place in the horizontal direction. The aperture varied from 1 to 8 mm, whereas the “fracture roughness” was created by combining JRC and 3D printing technology, where the Joint Roughness Coefficient (JRC) means a dimensionless measure of fracture surface roughness ranging from 0 to 20 (Figure 2). Based on the experimental observations, a semiempirical equation of the Forchheimer’s coefficients and dependent on hydraulic aperture and roughness was proposed.

2. Theories and Experimental Methodology

2.1. Theories Relevant to This Study

Models of fluid flow in porous media are widely used in multiple engineering and scientific researches. The traditional linear equation for flow in porous media is based on Darcy’s law:where is the fluid viscosity, is the flow velocity, is the permeability, and is the hydrodynamic pressure.

It is well known that Darcy’s law is not sufficient for describing accurately high Reynolds number flows. An inertial term can be added to Darcy’s equation, known as the Forchheimer term to account for the nonlinear behavior of the pressure difference versus velocity data. Forchheimer’s equation provides a general relation including this nonlinear effect:where and are the Forchheimer coefficients describing pressure losses due to viscous and inertial dissipation mechanisms, respectively [2, 18]. Where and , ( is the aperture of the idealized parallel smooth fracture; is the fluid density) is the intrinsic permeability, and is the non-Darcy coefficient or inertial resistance coefficient dependent on the geometrical properties of the medium [19]. Using dimensional analysis, Schrauf and Evans [20] rewrote (2) as where and are dimensionless coefficients. By comparing (3) to (2), one can obtain and .

2.2. Sample Preparation

In this study, fracture specimens with different joint roughness coefficient were fabricated by 3D printing technology according to the JRC standard profile curve proposed by Barton and Choubey [21], and a high-speed non-Darcy seepage experiment system is established based on this. The Darcy and non-Darcy seepage experiments of rough fracture with various fracture aperture are carried out in this system. The method of fracture generation is presented in the following three steps.

2.2.1. Rock Fracture Printing

The JRC profile curve proposed by Barton is scanned and plotted according to the original size. After obtaining the digital JRC standard curve, the 3d fracture model is built by 3D modeling of CAD software. Considering the seal ability and convenience of adjusting the crack width, 15 mm cushion layer is set on both sides of the rough fracture, and 25 mm flow transition zone is added in the inlet and outlet of the smooth plate. The 3D CAD model is transformed into STL format file and three-dimensional fracture plate model is made based on 3D printing technology (precision 0.1 mm) as shown in Figure 1.

2.2.2. Production of Cement Specimens

The outer diameter of 110 mm and wall thickness of 3.2 mm of a circular PVC pipe were used as a sample mold, with rough fracture plate in the center. Copper pipe was embedded in the mold upstream, central, and downstream for water pressure monitoring. Specimen pouring groove is shown in Figure 3. The test piece has a diameter of 10.36 cm and a total length of 35 cm. The rough fracture portion has a width of 8 cm and a length of 30 cm.

2.2.3. Test Device Assembly

The concrete sample (in Figure 3(b)) was fixed to the test stand, as shown in Figure 4. During the test, by setting different thickness gasket in the smooth edges (cushion area) reserved between the upper and lower plates of the sample to simulate different fracture aperture, the fracture aperture was set to 1, 2, 3, 5, and 8 mm, respectively. The upstream pressure was controlled by a pump and a regulator. The water tank is prepared between the inlet and sample to maintain the stability of the water flow within the specimen and, ultimately, the circulating water supply system which is made up of the storage tank, pumps, test stand, and other equipment (Figure 4). Due to the large flow, delta weir was used as a flow measurement device.

2.3. Experiment System Test

In order to verify the accuracy of the test equipment and the measurement system, we conduct the fracture flow test under the condition of laminar flow. The results for JRC = 0~2 and JRC = 16~18 low velocity seepage test in a rough fracture with aperture set as 3 mm is shown in Figure 5. It can be seen that when the hydraulic gradient is low (), there is a slight deviation from Darcy’s law between flow and hydraulic gradient at JRC = 16–18. Overall, the flow and hydraulic gradient satisfy a linear relationship, indicating that the fracture flow in this scenario is in the laminar flow stage. The water flow conforms to the linear Darcy's law; the second term of the right side of (2) can be ignored. The results of the relevant test show the accuracy of the results obtained from another angle.

3. Experiment Results

3.1. Relation between Flow within Smooth Fracture and Crack Opening

In order to explore the correlation between velocity, hydraulic gradient, and fracture aperture under high hydraulic gradient, high-speed non-Darcy seepage tests in a smooth fracture are carried out before carrying out rough fracture non-Darcy seepage test.

Figure 6 shows the velocity and hydraulic gradient curve for 4 fracture apertures. The relationship between the velocity and the hydraulic gradient for the large fracture aperture and high velocity can be well described by a power function, and the flow velocity and the hydraulic gradient have deviated from the linear relationship [2224]. Through fitting the test data point, the transform Izbash equation (4) can be used to analyze the results. Izbash’s equation describes the relationship between the fluid flow velocity and the hydraulic gradient is given a theoretical background arriving from a drag model. The exponent in Izbash’s equation [25, 26] has been obtained over a wide range of Reynolds’ numbers.where , are the fitting coefficients.

The result of analysis showed that the correlation coefficient of each curve is more than 0.989; = 0.98~3.30, value increases as the fracture aperture becomes larger. For nonlinear flow yet to approach fully turbulent state, changes between 1 and 0.5, = 0.5 for all apertures in this test representing a fully turbulent flow.

The flow rate is usually expressed in terms of the dimensionless Reynolds number, which quantifies the relative strength of inertia forces as compared to viscous forces. For fluid flow through fractures, Reynolds number [27, 28] can be defined aswhere is the characteristic length (hydraulic radius) of the idealized parallel smooth fracture.

The law of non-Darcy seepage in smooth fractures can be deduced according to the semiempirical theory of Prandtl turbulence [22]. The relationship between the average velocity of fractures and the hydraulic gradient and fracture aperture is obtained by mathematical derivation.where is the average velocity of the section, is the Carmen constant, and is the undetermined coefficient.

As shown in (6), by assigning , the product of the average velocity and the hydraulic gradient can be described as a function of the fracture aperture , as shown in the following:where , , are the fitting coefficient, which can be obtained from the test data.

According to the relationship shown in (7), a fitted curve (as shown in Figure 7) could be obtained by sorting out the test data in Figure 6. The results show that the correlation coefficient is 0.991 in our experiments. The values of , , and are estimated as 0.223, 0.44, and 51.47, respectively.

As shown in Figure 7, the velocity and hydraulic gradient in smooth fracture can be expressed as the expression of the fracture aperture in the case of complete turbulence.

3.2. Relation between JRC and Flow within Roughness Fracture
3.2.1. Basic Experimental Data

Based on the results of Figure 7, it can be concluded that is the expression of fracture aperture . In reality, the fracture in rock crack propagation is impossible to be smooth. The fracture surface of different rock cracks is different, so whether the correlation of smooth fracture flow can be directly applied to rough fracture is still an open question. Therefore, the fracture roughness is introduced, and, based on this, it is of great theoretical value to explore whether there is some relevance in velocity, hydraulic gradient, fracture aperture, and roughness.

In the above-mentioned smooth parallel plate high-speed fracture flow test, the minimum Reynolds number is 2769. Obviously it is in turbulent state. In rough fractures, the hydraulic radius is clearly less than the smooth fracture due to the roughness. In order to describe the relationship between flow velocity and hydraulic gradient in complete flow, the Forchheimer equation for Darcy and non-Darcy effect is still used as benchmark.

Based on the high-speed non-Darcy seepage test, the high-speed non-Darcy seepage tests with fracture aperture of 1, 2, 3, 5, and 8 mm are carried out to obtain the relationship between the hydraulic gradient and velocity in the fracture. Reynolds numbers are obtained in each case based on (5). According to the knowledge of the classic pipeline flow, when Reynolds number is greater than 2300, the flow can be considered as turbulence. So the non-Darcy turbulent can be distinguished.

For different JRC and fracture apertures, a series of repeated experiments are carried out, respectively, to obtain the flow rate under the corresponding hydraulic gradient. The experimental results show that when fracture apertures is small, the repeated experiments show greater deviation. When the aperture is 1 mm, the maximum deviation is 6.27%; when the aperture is 8 mm, the maximum deviation is only 2.34%. Correspondingly, the maximum deviations are 5.71%, 3.62%, and 2.89% respectively, when the aperture is 2 mm, 3 mm, and 5 mm. For the reason why the experiment deviation is larger when the aperture is smaller, it should be related to the reason that the aperture is adjusted under high water pressure.

From the relationship between the hydraulic gradient and the velocity in Figures 8(a)8(e), the relationship between the hydraulic gradient and the flow velocity obtained from each test deviates from the linear relationship, indicating that the flow of cracks does not follow Darcy’s law. It is non-Darcy flow. For the same fracture aperture, the greater the JRC is, the greater hydraulic gradient it will require and the stronger the degree of deviation from the linear relationship is.

3.2.2. Analysis Based on Power Function

Based on the least-squares method, the fitting curve is obtained according to the power function of (4). The correlation coefficient is 0.98967~0.99708. The fitting coefficient and correlation coefficient are summarized in Table 1. The exponent coefficient is between 0.60388 and 0.7815.

Data presented in Table 1 are the results for a 3 mm aperture fracture. The flow state includes the flow from the transition zone to the turbulence zone. Therefore, under different roughness conditions, value cannot be unified 0.5. After several fittings, when the hydraulic gradient index = 0.60~0.68, the correlation coefficient can be maintained above 0.9.

Modeled on (7), the relationship between (the ratio of flow rate to hydraulic gradient ) and JRC was plotted in Figure 9. It can be observed that there is a negative power function relationship between and JRC.

The fracture width of the test specimen was 0.08 m. According to the flow rate calculation equation, the empirical relationship between the unit discharge , the hydraulic gradient , and the rock joint roughness coefficient JRC can be obtained from the flow calculation equation:where is unit discharge (m3/s), is the hydraulic gradient, JRC is the rock joint roughness coefficient, respectively, and , , are the empirical coefficients obtained from test. Empirical coefficient = 3.209875~3.132625, = 0.30518~0.38066, and = 0.60~0.68.

According to (7), (8) is used to establish the empirical relationship between flow rate, hydraulic gradient, and roughness. However, due to the change of hydraulic gradient during the test process, the flow state in the fracture is transformed from transition state to turbulent state. The hydraulic gradient index cannot be used as a certain value. Although the empirical coefficient does not change much, it is impossible to fit the equation by using a certain value.

3.2.3. Forchheimer Equation Analysis

In order to further study seepage in rough fractures under high velocity conditions, the Forchheimer equation is introduced to analyze the experimental data.

The results obtained by fitting the Forchheimer equation are summarized in Table 2. The correlation coefficient of the experimental data which is fitted based on the Forchheimer equation is more than 0.98, and the fitting effect is better. It can be seen from Figure 10 that, as the flow velocity increases, the influence of inertial force becomes more and more obvious, and the relationship between hydraulic gradient and seepage velocity increasingly deviates from the linear relationship.

In order to analyze the effect of the fracture aperture on the seepage movement, the relationship between the monomial coefficient and the quadratic coefficient of the Forchheimer equation is summarized in Table 2. The fracture aperture is plotted, respectively, as shown in Figure 11. It can be concluded that the curve is consistent with the general relationship between the coefficient and the fracture aperture in (3). The monomial coefficient and the quadratic coefficient decrease with the increase of the fracture aperture. Overall, the conclusion is consistent with the results of Qian et al. [16] who uses ups and downs of the triangle and rectangle to characterize roughness. It is shown that the high-speed non-Darcy seepage in different roughness fractures still follows the traditional empirical equation.

(1) Relationship with Monomial Coefficient. According to the physical meaning of the Forchheimer equation, the monomial coefficient and the quadratic coefficient represent the influence of the viscous and inertial forces on the flow in the fracture flow, respectively. In order to further analyze the relationship between the monomial coefficient, the quadratic coefficient of the Forchheimer equation, and the fracture roughness, the relationship between the monomial coefficient and the JRC curve is plotted in Figure 12.

According to the basic characteristics of (3), the monomial coefficient is related to the reciprocal of the square of the fracture aperture. Since the influence of the JRC is taken into account, the actual fracture aperture is changed. Therefore, the ordinate is set to , in Figures 12(a) and 12(b). The abscissa is set to JRC, so as to analyze the corresponding relationship between JRC and monomial coefficient and fracture aperture.

As shown in Figures 12(a) and 12(b), due to the influence of roughness, which leads to a changing effective fracture aperture, the relationship between and reciprocal squares is difficult to be maintained, and the analysis of test data reveals if is used as -axis, when JRC is relatively small (⩽10). The relevant numerical analysis can be distributed in the narrow band; as the JRC becomes larger (12), is used as -axis, the relevant value can be distributed in a narrower band. At this point, for narrower bands, linear relationship can be used for segment description. At the meanwhile, the figure shows that, as JRC increases, the fracture aperture becomes smaller and the relationship between and is more sensitive. Obviously, it is consistent with the actual situation.

(2) Relation with Second-Order Coefficient. The quadratic term of the Forchheimer equation describes the extent to which the fluid flow deviates from the linear law. Therefore, the quadratic coefficient obtained by fitting Forchheimer equation plays an important role in the analysis of the extent of the non-Darcy flow deviating from the linear flow. It can be obtained from (3) thatThe non-Darcy influence coefficient of high-speed non-Darcy fluid is , where is the quadratic coefficient of Forchheimer obtained by fitting the equation and is the acceleration of gravity.

From the point of view of the trend, JRC is linearly related to the quadratic term coefficient for different fracture apertures. Under the same fracture aperture condition, with the increase of joint roughness coefficient, the non-Darcy influence coefficient of rough fracture increases; in the same joint roughness coefficient fissure, the non-Darcy influence coefficient decreases with the increase of the fracture aperture. It is shown that the extent of the non-Darcy in the high-speed seepage condition is inversely proportional to the fracture aperture, which is proportional to the joint roughness coefficient in the rough rock fracture.

As discussed above, the non-Darcy influence coefficient decreases with the increase of the fracture aperture for the same roughness fracture specimen. The unit of the non-Darcy influence coefficient is m−1 and unit of the fracture aperture is m; it is inferred that the non-Darcy coefficient of the fissure is inversely proportional to the fracture aperture. The curves of the relationship in Figure 11(b) are fitted for analysis under conditions of different joint roughness coefficient. Fit the equation withwhere represents the non-Darcy influence coefficient of the fracture, and the fitting results are summarized in Table 3.

As can be seen from the data in Table 3, with the increase of joint roughness coefficient, coefficient increases. The relationship between and the joint roughness coefficient JRC is plotted in Figure 13. The relationship between these two can be described as a linear relationship. Based on the least square method, the fitting relationship iswhere is the linear fitting coefficient.

The correlation coefficient is 0.92124 and the linear coefficient is = 0.042; the empirical relationship between the non-Darcy influence coefficient and the fracture and the joint roughness coefficient JRC can be obtained as follows:In the equation, is the non-Darcy influence coefficient, JRC is the joint roughness coefficient of fracture surface, and is the fracture aperture.

According to (12), as the joint roughness increases, the non-Darcy influence coefficient increases proportionally. With the increase of the fracture aperture, the non-Darcy influence coefficient is inversely proportional, so the ratio of joint roughness coefficient to the fracture aperture can be considered an important parameter to measure the high velocity non-Darcy characteristic of rough fracture flow.

According to the curve in Sections 3.2.3 and 3.2.3 and (3), the empirical equation of hydraulic gradient and velocity takes the effect of JRC into consideration under the Darcy-Non-Darcy seepage condition.

Clearly, due to the influence of roughness, the relationship between the monomial coefficient and the fracture aperture in (3) changes. When JRC is different, the head drop value of the point caused by the viscous force term and the inertia term change. Obviously when the crack propagation behavior occurs, according to the fracture aperture and roughness, the corresponding equation can be used.

3.3. Discussion

A high-speed non-Darcy infiltration experiment is developed. Based on the expression of the Forchheimer coefficient, a series of experiments are carried out to evaluate the influence of fault morphology and opening on non-Darcy flow. The experimental model involves two key geometric parameters, the fracture aperture () used to describe the intrinsic flow capacity and the roughness (JRC) used to describe the fracture morphology. The current X-ray computed tomography (CT) has also been successfully used for roughness parameter values based on contour measurements of the roughness height of the surface of the fracture wall in the laboratory [15].

It is notable that the focus of this study is on the laboratory scale evaluation of the Forchheimer equation coefficients, which has not yet been attempted to elevate the laboratory results to the field scale of natural fracturing systems. The flow and hydrodynamic gradient distribution of fractured geological media on a wide scale is a good research topic and has been extensively explored [10, 16] and remains to be explored. However, the proposed semiempirical Forchheimer equation in this study also reflects the effects of aperture and roughness parameters on hydraulic gradients and flow rates, and it is of important reference value for deterministic faults or flow headage and crack tip during crack propagation due to engineering disturbances and effective stress determination, which may provide important guidance for further development or expansion to field scale flow.

4. Conclusion

In this study, JRC is used to characterize the joint roughness, and the high-speed non-Darcy seepage test is carried out in 25 fractures of geotechnical soils combining 3D printing technique. It is mainly used to evaluate the Forchheimer equation coefficients in rough fractures. The empirical relationship between the coefficients A, B of Forchheimer equation and the aperture and roughness of fracture is first studied, which can be used to describe the high velocity non-Darcy flow in rock fractures.

Based on the results of the high-speed non-Darcy seepage test, the obtained pressure gradient and velocity curve show that the quadratic coefficient of the Forchheimer equation can be described as a function of the fracture aperture in the complete turbulence stage. Forchheimer’s law can be used to describe Non-Darcy flow behavior caused by inertia effect in deformable thick wall fractures.

On the basis of the smooth fracture test, the high-speed non-Darcy seepage was carried out. The influence of roughness on high-speed non-Darcy flow is analyzed by 25 experiments. The results show that the hydraulic gradient index = 0.60~0.68; the correlation coefficient can be maintained at 0.9 or more when considering the fracture roughness, according to the power exponent. If the Forchheimer equation is used to fit the curve, the correlation coefficient is more than 0.98, and it is clear that the Forchheimer equation has a strong adaptability to the high velocity non-Darcy flow.

The influence of aperture and roughness on Forchheimer coefficient is analyzed. It is observed that the influence of the fracture aperture on the nonlinear coefficients and is in an accordance with the existing conclusions. The roughness has a significant effect on the nonlinear coefficients and , and the JRC is in linear relationship with the coefficient . However, JRC also has an effect on an actual fracture aperture. Therefore, the equation requires appropriate correction. JRC and coefficient is linearly related. The analysis shows that when the fracture aperture is small, the impact of JRC is more significant.

The methods and results may be of interest to industrial geologists, hydrologists and engineers who focus on rock aquifer systems, and modelers who seek to validate more comprehensive numerical treatments.

Symbols

, :Forchheimer coefficients (–)
:Fitting coefficients (–)
:Undetermined coefficient (–)
:Dimensionless coefficients (–)
:Aperture of fracture (L)
:Pore pressure (ML−1T−2)
:Acceleration of gravity (LT−2)
, :Fitting coefficients (–)
, , :Fitting coefficients (–)
:Flow velocity (LT−1)
:Average velocity of the section (LT−1)
:Carmen constant (–)
:Hydraulic gradient (–)
:Permeability (LT−1)
:Flow rate (L3T−1)
:Unit discharge (LT−1)
:Fluid density (ML−3)
:Empirical coefficients (–)
:Kinetic viscosity (ML−1 T−1)
:Non-Darcy influence coefficient (L−1).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The research described in this paper was supported by the National Natural Science Foundation (51779083) and Fundamental Research Funds for the Central Universities (2016B08114).