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

Abstract

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 flow resistance leads to inaccurate results, especially when the velocity and cross angle are large. Therefore, it is necessary to find a portable method for calculation. Physical experiments of single fracture with different apertures (e = 0.77, 1.18, 1.97, 2.73 mm) were set up first to study the flow characteristics, showing obvious non-Darcian flow, which can be depicted by the Forchheimer equation when the flow velocity is sufficiently large. The computational fluid 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 coefficient > 0.99). Then, the CFD models of flexural crack with different angles from 30° to 150° were established to compute the pressure drop of flexural crack at different velocity. It was found that the local head loss of the flexural crack varied with the bending angle, and its coefficient was expressed by the deformation of the logistic equation. By using this model, as well as a frictional head loss equation fitted by Forchheimer equation, the head loss of crossed fissures with fixed fracture aperture could be easily calculated.

1. 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 applications, a quadratic equation (${J = \mathit{Av} + \mathit{Bv}}^2$) and an exponential equation (${J = \mathit{Cv}}^{-m}$) are used to under limited conditions [8,9,10,11].

The Forchheimer equation can express the characteristics of non-Darcian flow in fracture [12,13,14,15]:

$$-\nabla {\mathrm{P} = \mathrm{a}\mathrm{Q} + \mathrm{b}\mathrm{Q}}^2$$

(1)

where Q is volumetric flow velocity, a and b are model coefficients, and −∇P is the pressure gradient.

The factors affecting the flow of fissure fluid are complex and include the shape of the fissure, degree of contact, roughness, fluid viscosity, and external pressure [16,17,18,19]. Liu et al. [20] studied the law of fluid movement in a fracture network and found that the coefficient (a, b) of the Forchheimer equation decreases with an increase in the fracture aperture, and the error in model fitting can be reduced by expressing the coefficients of the Forchheimer equation as a power equation of fracture aperture [14]. Olson et al. [21] studied the effect of the ratio of the fracture opening to its height on the flow pattern, and Shu et al. [22] studied the variation in head loss with the width of the fissure and the velocity of flow when water flows through L-shaped fissures. Li et al. [23] studied the correlation between the head loss of a fluid in cross fissures and the width as well as roughness of the fissures.

This study proposes a portable calculation model where the aperture and the shape of flexural crack are considered to predict head loss directly. It is difficult to consider all factors when studying the characteristics of flow in fractures. In order to investigate the influence of each one, the synergistic effects of each condition need to be reduced, and the influence rules of each condition need to be researched separately. Fracture aperture and shape are two of the important parameters to generalize fracture development, which explains why only the two most significant parameters were considered here.

2. 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.

2.1. 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.

2.2. 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:

• 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 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.
• 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.
• 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. Range of pressure losses and water flow rate

  Experimental Results	Fracture Aperture
  	0.077 cm	0.118 cm	0.197 cm	0.273 cm
  Range of water head loss (cm)	6–92	6–80	3–50	2–20
  Range of water flow Q (mL/s)	2–26	8–50	19–93	26–97

  .

2.3. 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.

2.3.1. 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.

2.3.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:

$$\mathsf{\rho }\frac{{\mathrm{d}}_{\mathrm{k}}}{{\mathrm{d}}_{\mathrm{t}}}=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left[\left(\mathsf{\mu }+\frac{{\mathsf{\mu }}_{\mathrm{t}}}{{\mathsf{\sigma }}_{\mathrm{k}}}\right)\frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{i}}}\right]{+\mathrm{G}}_{\mathrm{k}}-\mathsf{\rho }\mathsf{\epsilon }$$

(2)

Transport equation of turbulence dissipation rate:

$$\mathsf{\rho }\frac{{\mathrm{d}}_{\mathsf{\epsilon }}}{{\mathrm{d}}_{\mathrm{t}}}=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left[\left(\mathsf{\mu }+\frac{{\mathsf{\mu }}_{\mathrm{t}}}{{\mathsf{\sigma }}_{\mathsf{\epsilon }}}\right)\frac{\partial \mathsf{\epsilon }}{\partial {\mathrm{x}}_{\mathrm{i}}}\right]{+\mathsf{\rho }\mathrm{C}}_1{\mathrm{S}\mathsf{\epsilon }-\mathsf{\rho }\mathrm{C}}_2\frac{{\mathsf{\epsilon }}^2}{\mathrm{k}+\sqrt{v\mathsf{\epsilon }}}$$

(3)

where k represents the dynamics of turbulence, ε is its dissipation rating; t is time; $\mathsf{\rho }$ is the density of the fluid; v is velocity; $\mathsf{\mu }$ is dynamic viscosity coefficient, since the experimental temperature is 12 °C, $\mathsf{\mu }$ = 1.2362; ${\mathsf{\mu }}_{\mathrm{t}}$ is turbulent dynamic viscosity; ${\mathrm{G}}_{\mathrm{k}}$ is the generation term of turbulent kinetic energy $\mathrm{k}$ due to the mean velocity gradient; $\mathrm{S}$ is the deformation velocity tensor; ${\mathrm{x}}_{\mathrm{i}}$ is coordinate component; and ω is the constant angular velocity. The values of ${\mathsf{\sigma }}_{\mathrm{k}}$ and ${\mathsf{\sigma }}_{\mathsf{\epsilon }}$ were 1 and 1.2, respectively. They represent Prandtl numbers corresponding to turbulent kinetic energy $\mathrm{k}$ and $\mathsf{\epsilon }$, respectively. The value of ${\mathrm{C}}_2$ was 1.9, whereas ${\mathrm{C}}_1$, ${\mathsf{\mu }}_{\mathrm{t}}$ can be defined as:

$${\mathrm{C}}_1=\mathit{max}\left[0.43,\frac{\eta }{\eta +5}\right]$$

(4)

$${\mathsf{\mu }}_{\mathrm{t}}={\mathsf{\rho }\mathrm{C}}_{\mathsf{\mu }}\frac{{\mathrm{k}}^2}{\mathsf{\epsilon }}$$

(5)

${\mathrm{C}}_{\mathsf{\mu }}$ represents model coefficients, the experimental as well as DNS data on the inertial sublayer of a channel or boundary layer flow suggest that ${\mathrm{C}}_{\mathsf{\mu }}$ = 0.09. The values of the parameters were recommended and by Launder et al. [33,34].

3. Results and Discussion

3.1. 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 $-\nabla \mathrm{P}$ corresponding to different flow Q under different fracture apertures was obtained and is analyzed in

Table 2. Regression of the linear and Forchheimer equations.

Category	Fracture Aperture (mm)	Fitting Equation	${\mathit{R}}^2$	Akaike Information Criterion (AIC)	Schwarz’s Criterion (SC)	AARE
Linear Equation	0.77	$-\nabla \mathrm{P}$ = 23.0001Q	0.9973	7.26	7.31	4.6%
	1.18	$-\nabla \mathrm{P}$ = 8.7935Q	0.8793	10.82	10.87	26.1%
	1.97	$-\nabla \mathrm{P}$ = 2.7505Q	0.8542	9.99	10.04	31.5%
	2.73	$-\nabla \mathrm{P}$ = 1.0793Q	0.8871	7.86	7.90	25.6%
Forchheimer Equation	0.77	$-\nabla \mathrm{P}$ = 21.1457Q + 0.0927Q2	0.9992	6.20	6.30	2.2%
	1.18	$-\nabla \mathrm{P}$ = 1.7394Q + 0.1761Q2	0.9957	7.61	7.71	9.1%
	1.97	$-\nabla \mathrm{P}$ = 0.4776Q + 0.0329Q2	0.9997	4.03	4.13	1.1%
	2.73	$-\nabla \mathrm{P}$ = 0.3109Q + 0.0098Q2	0.9976	4.18	4.26	3.0%

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].

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.

3.2. 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² ranging from 0.9643 to 0.9951. As is well-known, when the determination coefficient R² 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. Analysis using Pearson’s correlation coefficient.

Fracture Aperture (mm)	0.77	1.18	1.97	2.73
Pearson’s correlation coefficient	0.995	0.996	0.999	1

, 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.

3.3. 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]:

$$\frac{v^2}{2g}+z+\frac{\mathrm{P}}{\mathsf{\rho }g}=\mathit{const}$$

(6)

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]:

$$\frac{v_1^2}{{2g\mathsf{\alpha }}_1}{+z}_1+\frac{{\mathrm{P}}_1}{\mathsf{\rho }g}=\frac{v_2^2}{{2g\mathsf{\alpha }}_2}{+z}_2+\frac{{\mathrm{P}}_2}{\mathsf{\rho }g}{+h}_{\mathrm{l}}{+h}_{\mathrm{f}}$$

(7)

where v is flow velocity in the fracture, (m·s⁻¹); g is gravitational acceleration, 9.8 m/s²; $h_{\mathrm{f}}$ is the local head loss caused by the shape of the fissure (θ); $h_{\mathrm{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 $\mathsf{\Delta }H$, which can be written as:

$$\mathsf{\Delta }{H=h}_{\mathrm{l}}{+h}_{\mathrm{f}}$$

(8)

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_{{\mathrm{l}}_0}{ = h}_{{\mathrm{l}}_{\mathrm{i}}}$). When θ = 0, the local head loss caused by fracture shape was 0 ($h_{{\mathrm{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.

$$\mathsf{\Delta }{H}_{\mathrm{i}} - \mathsf{\Delta }{H}_0{ = h}_{\mathrm{f}}{}_{\mathrm{i}}$$

(9)

where ΔH₀ 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.

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.

3.4. 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. Pairwise correlation analysis between h_f 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 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_{\mathrm{f}}\left(v,\mathsf{\theta }\right)$ should consider all of them, and can be expressed as [43]:

$$h_{\mathrm{f}}=\xi \frac{v^2}{2g}$$

(10)

$\xi$ 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:

$$\xi =\frac{{2\mathit{gh}}_{\mathrm{f}}}{v^2}$$

(11)

${\xi }^{\ast }$ is the average value of $\xi$ at different velocities. As shown in Figure 7, the value of ${\xi }^{\ast }$ 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:

$$y = \frac{\mathrm{b}}{{1 + \mathrm{Ce}}^{-\mathrm{at}}}$$

(12)

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:

$${\xi \left(\mathsf{\theta }\right)=\mathrm{A}}_2+\frac{{\mathrm{A}}_1{-\mathrm{A}}_2}{1+{\left(\frac{\mathsf{\theta }}{\mathrm{b}}\right)}^p}$$

(13)

The crossed angle (θ) was the independent variable with a range of 0–180 degrees. ${\mathrm{A}}_1$, ${\mathrm{A}}_2$, $\mathrm{b}$, and p are the characteristic parameters of the equation. For our study, they were ${\mathrm{A}}_1$ = 0.1673, ${\mathrm{A}}_2$ = 3.7346, $\mathrm{b}$ = 88.0791, and p = 3.7179.

Based on the Equation (12), we can quickly calculate the local head loss:

$$h_{\mathrm{f}}(\mathsf{\theta },v)=\frac{v^2}{2g}\left(3.7346+\frac{0.1673-3.7346}{1+(\mathsf{\theta }/{88.0791)}^{3.7179}}\right)$$

(14)

where g is the gravity coefficient, with the value of 9.8 m/s².

3.5. 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² > 0.99, as follows:

$$J = 0.6803v + {1.1035v}^2$$

(15)

Considering the distance of frictional head loss, the calculation equation of frictional head loss can be obtained as:

$$h_{\mathrm{l}}=\left(0.6803v+{1.1035v}^2\right)\times L$$

(16)

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:

$${\mathrm{Q}}_{1-2}{=\mathrm{Q}}_{2-\mathrm{a}-3}{+\mathrm{Q}}_{2-\mathrm{b}-3}{=\mathrm{Q}}_{3-4}$$

(17)

Assuming that the flow velocity at the entrance (1#) and exit (4#) is v, the flow velocity in each flow path is:

$$v_{1-2}=2v_{2-\mathrm{a}-3}=2v_{2-\mathrm{b}-3}{=v}_{3-4}=v$$

(18)

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:

∆H_1-4 = h_l(1-4) + h_f(1-4)

(19)

where:

h_l(1-4) = h_l(1-2) + h_l(2-a) + h_l(a-3) + h_l(2-b) + h_l(b-3) + h_l(3-4)

(20)

h_f(1-4) = h_f(1-2-a) + h_f(2-a-3) + h_f(a-3-4) + h_f(1-2-b) + h_f(2-b-3) + h_f(b-3-4)

(21)

Among the Equations (19) and (20), h_l(1-2), h_l(2-a), h_l(a-3), h_l(2-b), h_l(b-3), and h_l(3-4) refer to the frictional head loss of the fracture segment of (1-2), (2-a), (a-3), (2-b), (b-3), and (3-4); h_f(1-2-a), h_f(2-a-3), h_f(a-3-4), h_f(1-2-b), h_f(2-b-3), and h_f(b-3-4)refer to the local head loss at the positions where the direction of flow changes. Parameters of each fracture or intersection point are listed in

Table 5. Parameters of each fracture for frictional head loss calculation.

Fracture	Flow Velocity (m/s)	Length L (v) (m)	Frictional Head Loss (m)	Sum of Frictional Head Loss (m)
hl(1-2)	v	0.05	0.05 (av + bv2)	0.2av + 0.15bv2 (a = 0.6803, b = 1.1035)
hl(2-a)	0.5v	0.05	0.05 (0.5av + 0.25bv2)
hl(a-3)	0.5v	0.05	0.05 (0.5av + 0.25bv2)
hl(2-b)	0.5v	0.05	0.05 (0.5av + 0.25bv2)
hl(b-3)	0.5v	0.05	0.05 (0.5av + 0.25bv2)
hl(3-4)	v	0.05	0.05 (av + bv2)

and

Table 6. Parameters of each intersection point for local head loss calculation.

Intersection Point	Flow Velocity (m/s)	Bending Angle θ (°)	$\mathit{\xi }$	Sum of Local Head Loss (m)
hf(1-2-a)	0.5v	45	0.4387	0.7250v2/g
hf(2-a-3)	0.5v	90	2.0224
hf(a-3-4)	0.5v	45	0.4387
hf(1-2-b)	0.5v	45	0.4387
hf(2-b-3)	0.5v	90	2.0224
hf(b-3-4)	0.5v	45	0.4387

.

The head losses calculated are listed in

Table 7. Total head loss calculation of network fracture.

Flow Velocity (m/s)	0.2	0.4	0.8	1.2	1.5	2	3	5
Frictional head loss (m)	0.0338	0.0809	0.2148	0.4016	0.5765	0.9342	1.8979	4.8184
Local head loss (m)	0.0030	0.0118	0.0473	0.1065	0.1664	0.2959	0.6658	1.8494
Total head loss (m)	0.0368	0.0927	0.2621	0.5082	0.7430	1.2301	2.5637	6.6678

.

From 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.

4. Conclusions

• 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 $\mathsf{\Delta }H\left(\mathsf{\theta },v\right)$.
• 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.