A Fast Calculation Model for Local Head Loss of Non-Darcian Flow in Flexural Crack

: Local head loss caused by fracture intersection is often ignored because there has not been a simple method to calculate it until now. Relevant research shows that neglecting the local ﬂow resistance leads to inaccurate results, especially when the velocity and cross angle are large. Therefore, it is necessary to ﬁnd a portable method for calculation. Physical experiments of single fracture with di ﬀ erent apertures (e = 0.77, 1.18, 1.97, 2.73 mm) were set up ﬁrst to study the ﬂow characteristics, showing obvious non-Darcian ﬂow, which can be depicted by the Forchheimer equation when the ﬂow velocity is su ﬃ ciently large. The computational ﬂuid dynamics (CFD) software ANSYS FLUENT was used to build numeric simulation models. A good correlation between CFD simulation results and physical experiment results was found (Pearson’s correlation coe ﬃ cient > 0.99). Then, the CFD models of ﬂexural crack with di ﬀ erent angles from 30 ◦ to 150 ◦ were established to compute the pressure drop of ﬂexural crack at di ﬀ erent velocity. It was found that the local head loss of the ﬂexural crack varied with the bending angle, and its coe ﬃ cient was expressed by the deformation of the logistic equation. By using this model, as well as a frictional head loss equation ﬁtted by Forchheimer equation, the head loss of crossed ﬁssures with ﬁxed fracture aperture could be easily calculated.


Introduction
Fracture fluid in nature usually moves at a high speed. When its velocity reaches a certain value, the loss of fluid velocity is dominated by the inertial force, and this phenomenon is called the high speed non-Darcian seepage [1]. Basak et al. [2] and Madhav found that high speed non-Darcian flow is prevalent in hydraulic engineering, whereas Tartakovsky et al. [3] observed the same result in confined aquifers. Research on the law of non-Darcian seepage of high-speed fluids forms the basis of the effective prediction of complex fracture seepage in rock masses and is key to ensuring safety and ecological stability in engineering construction [4]. To describe the non-Darcian fluid in fissures, Zimmerman et al. [5] established a single sandstone fissure to determine the Reynolds number boundary for the best characterization equation of non-Darcian flow. Qian et al. [6] studied the evolution of flow from Darcian to non-Darcian in fissures under confined and unconfined condition and concluded that the critical Reynolds number decreases gradually with the increase in fissure width. Quinn et al. [7] carried out an experimental study on non-Darcian flow in fissures of a dolomite aquifer and quantitatively analyzed its hydraulic gradient and velocity.
At the microscopic scale, the Navier-Stokes equation is often used to study the relationship between the velocity and pressure of Newtonian fluids in cracks. Because of its complex calculation, however, it is difficult to apply it to engineering practice. To improve the efficiency of engineering

Experiments and Methodology
The fracture aperture was first used as the control condition to construct an indoor physical model, and head loss and velocity were observed to analyze the characteristics of the fluid in the single fissure. Following the verification of the computational fluid dynamics (CFD) model in simulation of characteristics of the fluid in the single fissure with parameters most used [24][25][26][27][28], the relationship between the shape of the fracture (corner θ) and head loss was studied using it.

Laboratory Test Device
The test device shown in Figure 1 consists of a water intake system, fissure system, pressure measurement system, and width measurement system. The water intake system was a set of water distribution devices that controlled the water intake head and ensured the stability of water inflow by setting an overflow port on the water distributor. The main fractured rock mass of the fractured system was composed of marble purchased from the stone market (50 × 30 × 1.8 cm), and its width could be adjusted to 0.77, 1.18, 1.97, and e = 2.73 mm, respectively, during the experiment. The piezometric system was composed of piezometric tubes and observation plates buried in different sections. The measurement system consisted of a device to measure the width of the fracture aperture and a flow-measuring device. The crack photos and ruler were taken by a 1080P HDMI digital microscope are magnified and read on the computer, and then the gap width of each section was accurately measured by using the image recognition technology of the supporting software, whose accuracy can be controlled at 0.01 mm. The flow-measuring device used a stopwatch, 1000 mL or 2000 mL beaker, and an electronic balance to determine the flow velocity. The accuracy of the electronic balance was 0.01 g. The climate did not change significantly, and the water temperature fluctuated slightly during the test, around 12 • C. system, which can control the water intake head and ensure the stability of water inflow; ② is fissure system, of which the fracture width can be adjusted to 0.77, 1.18, 1.97, and e = 2.73 mm; ③ is the piezometric system; ④ is fracture aperture measurement system, the major instrument is a 1080P HDMI digital microscope; ⑤ is flow measurement system.

Laboratory Test Method
Taking the fracture aperture and water pressure as control conditions, the system was reorganized when the fracture aperture was changed. The construction and operation of the device were as follows: 1. Assembling: Once the stone had been cut into two pieces in a professional stone crushing workshop, we measured and arranged it in the laboratory and installed the inlet and outlet pipe as well as pressure tubes. The fracture aperture was measured after some necessary adjustment and fixing of the stones, and then all of the interface was sealed with acrylic glue. When the test  1 is the water intake system, which can control the water intake head and ensure the stability of water inflow; 2 is fissure system, of which the fracture width can be adjusted to 0.77, 1.18, 1.97, and e = 2.73 mm; 3 is the piezometric system; 4 is fracture aperture measurement system, the major instrument is a 1080P HDMI digital microscope; 5 is flow measurement system.

Laboratory Test Method
Taking the fracture aperture and water pressure as control conditions, the system was reorganized when the fracture aperture was changed. The construction and operation of the device were as follows:

1.
Assembling: Once the stone had been cut into two pieces in a professional stone crushing workshop, we measured and arranged it in the laboratory and installed the inlet and outlet pipe as well as pressure tubes. The fracture aperture was measured after some necessary adjustment and fixing of the stones, and then all of the interface was sealed with acrylic glue. When the test Connection: After checking the sealing of the cracked plate, the water inlet was connected to the water intake system, and the pressure gauge pipe was connected to the pressure measurement system. It was ensured that no bubbles were in the pressure gauge pipe during the experiments. 3.
Test: Water pressure was adjusted through moving the height of the water tank, and more than 5 min was needed to wait for piezometric head stabilization. The fluid flow was measured by the volumetric method mentioned previously. The range of test conditions is shown in Table 1.

Numeric Simulation Model
CFD numerical calculation module can be used to simulate the laminar and turbulent flows of Newtonian and non-Newtonian fluids in hydraulics. Thakur et al. [29] used CFD simulations to study the influence of the size of Newtonian and non-Newtonian fluids on their flow fields and power consumption, Aubin et al. [30] modeled and analyzed turbulence in a stirred tank based on a CFD model, and Apsley et al. [31] used a CFD model to study the effects of roughness on turbulent flow. ANSYS FLUENT 19.1 (ANSYS, Canonsburg, PA, USA) was employed to construct a CFD model to simulate the characteristics of flow in fractures in this study.

Conceptual Model
A single fracture model with length of 100 mm, height of 20 mm, and aperture of 1 mm was built in FLUENT, while the bending angle of flexural crack was set to 0, 30, 45, 60, 90, 120, and 135 degrees. Because the viscous force had little influence on the characteristics of motion of high speed flow, the model was set up based on turbulent state of an ideal fluid, and the velocity at the inlet as well as a pressure of zero at the outlet were set as boundary conditions. The structure of the numerical model is shown in Figure 2.
Water 2020, 12, x FOR PEER REVIEW 4 of 15 was complete, the sealant was cleaned and the width of the gap readjusted to another one. The stones were then fixed and sealed again. We planned to set the gap widths within the range of 0.05-0.10, 0.10-0,15, 0.15-0.2, and 0.2-0.3 cm, while the actual measured gap widths were 0.77, 1.18, 1.97, 2.73 mm, respectively. 2. Connection: After checking the sealing of the cracked plate, the water inlet was connected to the water intake system, and the pressure gauge pipe was connected to the pressure measurement system. It was ensured that no bubbles were in the pressure gauge pipe during the experiments. 3. Test: Water pressure was adjusted through moving the height of the water tank, and more than 5 min was needed to wait for piezometric head stabilization. The fluid flow was measured by the volumetric method mentioned previously. The range of test conditions is shown in Table 1.

Numeric Simulation Model
CFD numerical calculation module can be used to simulate the laminar and turbulent flows of Newtonian and non-Newtonian fluids in hydraulics. Thakur et al. [29] used CFD simulations to study the influence of the size of Newtonian and non-Newtonian fluids on their flow fields and power consumption, Aubin et al. [30] modeled and analyzed turbulence in a stirred tank based on a CFD model, and Apsley et al. [31] used a CFD model to study the effects of roughness on turbulent flow. ANSYS FLUENT 19.1 (ANSYS, Canonsburg, PA, USA) was employed to construct a CFD model to simulate the characteristics of flow in fractures in this study.

Conceptual Model
A single fracture model with length of 100 mm, height of 20 mm, and aperture of 1 mm was built in FLUENT, while the bending angle of flexural crack was set to 0, 30, 45, 60, 90, 120, and 135 degrees. Because the viscous force had little influence on the characteristics of motion of high speed flow, the model was set up based on turbulent state of an ideal fluid, and the velocity at the inlet as well as a pressure of zero at the outlet were set as boundary conditions. The structure of the numerical model is shown in Figure 2.

Mathematical Model
The CFD model was established based on the size and parameters of the physical model. High speed non-Darcian flow in the field of physical experiments was simulated by using the realizable k-ε model, which was proposed by Launder and Spalding in 1972 and subsequently constructed by Shih et al. [32] in 1995. Based on the time-averaged continuity equation, the equations of turbulent dynamic transport and the turbulence dissipation rating were established.
Turbulence dynamic transport equation: Transport equation of turbulence dissipation rate: where k represents the dynamics of turbulence, ε is its dissipation rating; t is time; is the density of the fluid; v is velocity; µ is dynamic viscosity coefficient, since the experimental temperature is 12 • C, µ = 1.2362; µ t is turbulent dynamic viscosity; G k is the generation term of turbulent kinetic energy k due to the mean velocity gradient; S is the deformation velocity tensor; x i is coordinate component; and ω is the constant angular velocity. The values of σ k and σ ε were 1 and 1.2, respectively. They represent Prandtl numbers corresponding to turbulent kinetic energy k and ε, respectively. The value of C 2 was 1.9, whereas C 1 , µ t can be defined as: C µ represents model coefficients, the experimental as well as DNS data on the inertial sublayer of a channel or boundary layer flow suggest that C µ = 0.09. The values of the parameters were recommended and by Launder et al. [33,34].

Data Analysis of Laboratory Physical Model Test
In the physical experiment of single fracture, aperture e (0.77, 1.18, 1.97, 2.73 mm) was used as the test variable. The pressure gradient −∇P corresponding to different flow Q under different fracture apertures was obtained and is analyzed in Table 2 and plotted in Figure 3.    All of the scatterplots of the volume of flow and the pressure gradient are similar to a parabola except at e = 0.77 mm, which is in accordance with the results that reported by Shu et al. [22], Wang et al. [35], and Akinbodewa et al. [36]. They show obvious non-Darcian flow in the test, since the shape of the line is not linear but parabolic. The transition from linear to nonlinear behavior is like that described by Andrade et al. [37].

Comparison between Experiment and Simulation
The linear equation does not perform well when used to fit the relationship between the volume of flow and the pressure gradient except at e = 0.77 mm (Absolute Average Relative Error, AARE L = 4.6-31.5%), while the Forchheimer equation fits well with larger determination coefficients (AARE F = 1.1-9.1%). The sum of squares of residual errors of the Forchheimer equation was one to three orders of magnitude smaller than that of the linear Equation and has smaller values of Akaike information criterion (AIC) and Schwarz's criterion (SC). Thus, the Forchheimer equation is more suitable for describing the law of flow movement in fissures.

Comparison between Experiment and Simulation
The realizable k-ε model was used to calculate the pressure drop in a single fracture with different apertures. The results of the CFD model were compared with those of the laboratory test at the same fracture aperture. Figure 4 shows that all results of the CFD model were close to those of the physical test, with the determination coefficient R 2 ranging from 0.9643 to 0.9951. As is well-known, when the Water 2020, 12, 232 7 of 15 determination coefficient R 2 of a non-intercepting fitting curve is close to one, the two sets of data are similar.
The realizable k-ε model was used to calculate the pressure drop in a single fracture with different apertures. The results of the CFD model were compared with those of the laboratory test at the same fracture aperture. Figure 4 shows that all results of the CFD model were close to those of the physical test, with the determination coefficient R 2 ranging from 0.9643 to 0.9951. As is wellknown, when the determination coefficient R 2 of a non-intercepting fitting curve is close to one, the two sets of data are similar. To validate the above conclusion, Pearson's correlation coefficient was employed for further analysis. As shown in Table 3, all Pearson's correlation coefficients were close to one, indicating a strong correlation between the CDF modeling results and the numerical results, which is consistent with Figure 4. This means that the CFD model accurately simulated the flow of fracture fluid under specified conditions, consistent with the results of Wang et al. [38]. From the above analysis, it can be inferred that the CFD model was useful for studying the characteristics of flow in fractures, and it has potential in flexural crack modelling.

Effect of Fracture Shape on Head Loss
It is known that the local head loss changes when the fracture aperture changes. In order to study the relationship between the fracture shape (corner θ) and the local head loss, the fracture width is should be fixed first.
Bernoulli's equation represents the conservation of mechanical energy of an ideal fluid and can represent the process of constant energy change of the fluid per unit mass, the relation of element flow on any section, as shown in Equation (5) [39]: Based on this, the energy equation of the real fluid with energy loss at the cross-section of upstream and downstream can be written as follows [28]: To validate the above conclusion, Pearson's correlation coefficient was employed for further analysis. As shown in Table 3, all Pearson's correlation coefficients were close to one, indicating a strong correlation between the CDF modeling results and the numerical results, which is consistent with Figure 4. This means that the CFD model accurately simulated the flow of fracture fluid under specified conditions, consistent with the results of Wang et al. [38]. From the above analysis, it can be inferred that the CFD model was useful for studying the characteristics of flow in fractures, and it has potential in flexural crack modelling.

Effect of Fracture Shape on Head Loss
It is known that the local head loss changes when the fracture aperture changes. In order to study the relationship between the fracture shape (corner θ) and the local head loss, the fracture width is should be fixed first.
Bernoulli's equation represents the conservation of mechanical energy of an ideal fluid and can represent the process of constant energy change of the fluid per unit mass, the relation of element flow on any section, as shown in Equation (5) [39]: Based on this, the energy equation of the real fluid with energy loss at the cross-section of upstream and downstream can be written as follows [28]: Water 2020, 12, 232 where v is flow velocity in the fracture, (m·s −1 ); g is gravitational acceleration, 9.8 m/s 2 ; h f is the local head loss caused by the shape of the fissure (θ); h l is frictional head loss; "1" and "2" are cross-sectional symbols of upstream and downstream, respectively, and α is the total flow correction factor. For fully developed turbulent flow, α ≈ 1. The total head loss from upstream to downstream is set to ∆H, which can be written as: To study the local head loss caused by the shape of the fracture, a comparative numerical model was constructed with a fracture aperture of 1 mm, length of 100 mm, width of 20 mm, and the value of θ at zero. Under the same conditions, we established seven numerical models with θ = 30, 45, 60, 90, 120, 135, and 150 degrees. The contour maps of pressure at flexural crack with different θ can be seen in Figure 5.
where v is flow velocity in the fracture, (m·s −1 ); g is gravitational acceleration, 9.8 m/s 2 ; h f is the local head loss caused by the shape of the fissure (θ); h l is frictional head loss; "1" and "2" are crosssectional symbols of upstream and downstream, respectively, and α is the total flow correction factor. For fully developed turbulent flow, α ≈ 1. The total head loss from upstream to downstream is set to ∆H, which can be written as: To study the local head loss caused by the shape of the fracture, a comparative numerical model was constructed with a fracture aperture of 1 mm, length of 100 mm, width of 20 mm, and the value of θ at zero. Under the same conditions, we established seven numerical models with θ = 30, 45, 60, 90, 120, 135, and 150 degrees. The contour maps of pressure at flexural crack with different θ can be seen in Figure 5. In the CFD numerical model test, except for the shape of the fracture (θ), the boundary conditions and parameters of each model were the same, because of which the frictional head loss of each model was the same (h l 0 = h l i ). When θ = 0, the local head loss caused by fracture shape was 0 (h f 0 = 0). By comparing the results of these models, the local head loss caused at different values of θ at different flow velocities was calculated.
∆H i − ∆H 0 = h fi (9) where ∆H 0 is the head loss when θ = 0, and ∆H i is that when θ = i. The local head loss of each flexural crack at different angles (θ) is shown in Figure 6. It is clear that local head loss caused by fracture shape increased with the flow velocity and the crossed angle. In the CFD numerical model test, except for the shape of the fracture (θ), the boundary conditions and parameters of each model were the same, because of which the frictional head loss of each model was the same (h l 0 = h l i ). When θ = 0, the local head loss caused by fracture shape was 0 (h f 0 = 0). By comparing the results of these models, the local head loss caused at different values of θ at different flow velocities was calculated.
where ΔH0 is the head loss when θ = 0, and ΔHi is that when θ = i.
The local head loss of each flexural crack at different angles (θ) is shown in Figure 6. It is clear that local head loss caused by fracture shape increased with the flow velocity and the crossed angle. The velocity distribution in the water-crossing section is more uneven in non-Darcian flow, and there is shear stress between the layers of flow [40]. When the fluid flowed through the turning point of the flexural crack in a state of turbulence where the inertial force played a dominant role, the flow could not change direction abruptly along the fracture shape, with the consequence that the mainstream of fluids was separated from the side-wall of the fracture, and vortices were generated [41]. The eddy current and mainstream were superimposed, and spiral motion occurred downstream of the corner. In this motion, the generation, collision, merging, splitting, and disappearance of the eddy current consumed a large amount of the fluid's mechanical energy, which was the local head loss caused by the shape of the fracture.
It has been found that the larger θ is, the more active the motion of the eddy current is, and the greater the head loss is [42]. This is consistent with the results of this study. Some researchers have observed that initial pressure on the fluid controls its instantaneous flow rate in pipes. When the shape of the pipe remains constant, the local resistance loss increases with the Reynolds number [43], which agrees with the results of this study. Therefore, flow velocity and fracture shape are the important factors affecting local head loss. The velocity distribution in the water-crossing section is more uneven in non-Darcian flow, and there is shear stress between the layers of flow [40]. When the fluid flowed through the turning point of the flexural crack in a state of turbulence where the inertial force played a dominant role, the flow could not change direction abruptly along the fracture shape, with the consequence that the mainstream of fluids was separated from the side-wall of the fracture, and vortices were generated [41]. The eddy current and mainstream were superimposed, and spiral motion occurred downstream of the corner. In this motion, the generation, collision, merging, splitting, and disappearance of the eddy current consumed a large amount of the fluid's mechanical energy, which was the local head loss caused by the shape of the fracture.
It has been found that the larger θ is, the more active the motion of the eddy current is, and the greater the head loss is [42]. This is consistent with the results of this study. Some researchers have observed that initial pressure on the fluid controls its instantaneous flow rate in pipes. When the shape of the pipe remains constant, the local resistance loss increases with the Reynolds number [43], which agrees with the results of this study. Therefore, flow velocity and fracture shape are the important factors affecting local head loss.

Establishment of Fast Calculation of Local Head Loss
A pairwise correlation analysis between the independent variables (θ and v) and the dependent variable (local head loss h f ) was quantitatively represented using the Pearson coefficient, as shown in Table 4. Table 4. Pairwise correlation analysis between h f and θ as well as v.

Parameter Facture Shape (θ) Flow Velocity (v)
Local head loss (h f ) Pearson's correlation coefficient 0.296 0.758 Significant (bilateral) P 0.003 0 From the above tables, it is clear that Pearson's correlation coefficient had the value P v > P θ , indicating that the influence of v on local head loss was greater than that of θ. As the flow velocity and fracture shape are positively correlated with head loss, the calculation model, h f (v, θ) should consider all of them, and can be expressed as [43]: ξ is the coefficient of local head loss that is usually related to the geometry of the member and the nature of the fluid flowing through the member, and it can be rewritten as: ξ * is the average value of ξ at different velocities. As shown in Figure 7, the value of ξ * is nonlinearly positively correlated with θ, and error caused by changes in the flow velocity is minimal. Therefore, the calculation model for the local head loss coefficient is a function related only to the shape of the fracture but not to velocity.
Water 2020, 12, x FOR PEER REVIEW 10 of 15

Establishment of Fast Calculation of Local Head Loss
A pairwise correlation analysis between the independent variables (θ and v) and the dependent variable (local head loss hf) was quantitatively represented using the Pearson coefficient, as shown in Table 4. Table 4. Pairwise correlation analysis between hf and θ as well as v.

Parameter
Facture Shape (θ) Flow Velocity (v) Local head loss (hf) Pearson's correlation coefficient 0.296 0.758 Significant (bilateral) P 0.003 0 From the above tables, it is clear that Pearson's correlation coefficient had the value Pv > Pθ, indicating that the influence of v on local head loss was greater than that of θ. As the flow velocity and fracture shape are positively correlated with head loss, the calculation model, h f v, should consider all of them, and can be expressed as [43]: ξ is the coefficient of local head loss that is usually related to the geometry of the member and the nature of the fluid flowing through the member, and it can be rewritten as: ξ  is the average value of ξ at different velocities. As shown in Figure 7, the value of ξ  is nonlinearly positively correlated with θ, and error caused by changes in the flow velocity is minimal. Therefore, the calculation model for the local head loss coefficient is a function related only to the shape of the fracture but not to velocity. The curve characteristic observed is S-shaped, which is somewhat similar to a logistic function. The logistic function was first used to reflect a nutritional relationship to population mathematical model in 1838, but it has been found useful in many fields. The original logistic function is: (12) where t is the time variable, and a, b, and C are the model constants. The curve characteristic observed is S-shaped, which is somewhat similar to a logistic function. The logistic function was first used to reflect a nutritional relationship to population mathematical model in 1838, but it has been found useful in many fields. The original logistic function is: where t is the time variable, and a, b, and C are the model constants. The deformation of logistic equation [44] was used to fit the head loss coefficient as follows: The crossed angle (θ) was the independent variable with a range of 0-180 degrees. A 1 , A 2 , b, and p are the characteristic parameters of the equation. For our study, they were A 1 = 0.1673, A 2 = 3.7346, b = 88.0791, and p = 3.7179.
Based on the Equation (12), we can quickly calculate the local head loss: where g is the gravity coefficient, with the value of 9.8 m/s 2 .

The Use of Fast Calculation Equation of Head Loss
During fluid flow, the head loss includes the sum of the frictional head loss and the local head loss listed in Equation (7), while in Section 3.1, we found that "∇P-Q" can be well characterized by the Forchheimer equation. Therefore, the frictional head loss can be calculated according to the Forchheimer equation obtained previously. From numeric simulation, the fitting Forchheimer equation between hydraulic gradient (J) and flow velocity with fracture aperture of 1 mm and transverse width of 2 cm was obtained, with R 2 > 0.99, as follows: Considering the distance of frictional head loss, the calculation equation of frictional head loss can be obtained as: where L is the length of the fracture. The above quantitative equations (Equations (13) and (15)) can be applied to the calculation of head loss when the crack with aperture is 1 mm and transverse width is 2 cm. To apply this, a network fissure with four nodes 1#-4# was set as shown in Figure 8. The deformation of logistic equation [44] was used to fit the head loss coefficient as follows: The crossed angle (θ) was the independent variable with a range of 0-180 degrees. A 1 , A 2 , b, and p are the characteristic parameters of the equation. For our study, they were A 1 = 0.1673, A 2 = 3.7346, b = 88.0791, and p = 3.7179.

The Use of Fast Calculation Equation of Head Loss
During fluid flow, the head loss includes the sum of the frictional head loss and the local head loss listed in Equation (7), while in Section 3.1, we found that "∇P-Q" can be well characterized by the Forchheimer equation. Therefore, the frictional head loss can be calculated according to the Forchheimer equation obtained previously. From numeric simulation, the fitting Forchheimer equation between hydraulic gradient (J) and flow velocity with fracture aperture of 1 mm and transverse width of 2 cm was obtained, with R 2 > 0.99, as follows: Considering the distance of frictional head loss, the calculation equation of frictional head loss can be obtained as: where L is the length of the fracture. The above quantitative equations (Equations (13) and (15)) can be applied to the calculation of head loss when the crack with aperture is 1 mm and transverse width is 2 cm. To apply this, a network fissure with four nodes 1#-4# was set as shown in Figure 8. According to Kirchhoff's law, the flow equation at the node was established: According to Kirchhoff's law, the flow equation at the node was established: Assuming that the flow velocity at the entrance (1#) and exit (4#) is v, the flow velocity in each flow path is: In Figure 7, We can divide the network fracture into multiple single fractures and local loss occurrence points, so the head loss of the fracture network is the sum of them: where: Among the Equations (19) and (20)    The head losses calculated are listed in Table 7.  Table 7, it can be seen that local head loss accounts for a small proportion when the flow velocity is small, but it grows fast as the flow velocity increases because local head loss is proportional to the square of speed. Water inrush from the underground projects, such as tunnels, is always very fast, so the local head loss is necessary to be taken into consideration when doing water inflow predication. In addition, parameters of Equations (13) and (15) need to be recalibrated before use because they are obtained from the specific conditions when the fracture aperture is 1 mm and its width is 20 mm.

•
Based on the physical experiment, non-Darcian characteristics were found when the speed was high enough in fracture. The Forchheimer equation fits the experimental data well between pressure gradient and volume of flow. The fracture width and fluid velocity, tested in the physical experiment, were important factors affecting the fluid flow pattern.

•
The non-Darcian flow could be well represented by CFD model in a single fracture. It was employed to simulate the behavior of flow in flexural crack and obtain the result that local head loss was closely related to flow velocity and fracture shape, which can be represented as ∆H(θ, v).

•
The local head loss coefficient is almost independent of speed, while it has a very strong relationship with the bending angle θ of flexural crack. It can be treated as a logistic function of θ, such as ζ (θ).

•
Friction head loss and local head loss can be calculated separately when the fracture system is complicated, but it is strongly recommended to take the local head loss into consideration when the flow speed is high, as it accounts for a large proportion of the total head loss.