# A Coupled Darcy-Forchheimer Flow Model in Fractured Porous Media
## 5 Nonlinear Seepage Analysis of Tunnel  
_Model Setting and Calculating Parameters_  
Due to the complex geological conditions and large water pressure of buried fracture tunnels with rich water depth, the traditional calculation of cubic law would result in a large deviation from the field test. Therefore, it is necessary to analyze the seepage field of the tunnel surrounding the rock by using the Darcy-Forchheimer coupled flow method in the fractured rock mass.  
As shown in Figure 6, a numerical model of the buried tunnel with rich water depth is divided into 10,000 quadrilateral elements, with a total number of 10,201 nodes and a size of 0.8 m. The length and height of the calculation model are 80 m and 80 m, respectively, according to the geological conditions of the site and the influence of the boundary. The tunnel is a circular tunnel with a radius of 3 m and the coordinates of the center of the circle are (40,40). The tunnel does not consider the support measures and implements internal drainage, so the tunnel water pressure \(p\) is 0 Pa. The calculation parameters were selected from the engineering geology report. The water density was set as 1000 kg/m\({}^{3}\), the permeability coefficient of rock mass was \(1\times 10^{-8}\) m\({}^{2}\), the porosity was 0.3, and the fracture permeability coefficient was \(1\times 10^{-5}\) m\({}^{2}\). The viscosity of the fluid is \(10\times 10^{-3}\) Pa \(\times\) s.  
Figure 5: Relationship of error and mesh size for intersecting fracture case.  
In the simulation, it is assumed that the water supply of the surrounding rock far site is sufficient, and the underground water surface does not decrease with the tunnel drainage. Therefore, the groundwater level remains constant, and the water pressure is 0 Pa. The seepage field of the tunnel is linearly distributed along the numerical direction under the action of gravity and the gradient is fluid gravity. Therefore, the water pressure at the bottom is 8 \(\times\) 10\({}^{5}\) Pa (80 m head), and the water pressure on both sides is linear, as shown in Figure 6.  
To simplify the calculation, the tunnel model does not consider the geological structure, such as strata and faults, and only simulates the influence of random fractures in surrounding rock. The fracture model was generated by the random Mento-Carlo method [33; 34; 35; 36]. The position of the fracture obeys uniform distribution, and the length obeys exponential distribution,  
\[f(l)=\lambda e^{(-\lambda l)} \tag{20}\]  
where \(f(l)\) is the probability density function of fracture length, and \(\lambda\) is the exponential distribution parameter. The length of fracture is in the range of \([l_{min},\,l_{max}]\). Another important parameter is the direction of the fracture. In the three-dimensional fracture, it refers to the inclination and dip angle of the fracture, while in the two-dimensional fracture system, it can be expressed by the angle of the fracture \(\theta\). The angle \(\theta\) of fracture can be simulated by the von-Mises distribution,  
\[f(\theta)=\frac{e^{\kappa\cos{(\theta-\mu)}}}{2\pi I_{0}(\kappa)} \tag{21}\]  
\[I_{n}(\kappa)=\frac{1}{\pi}{\int_{0}^{\pi}\exp(\kappa\cos(x))\cos(nx)dx} \tag{22}\]  
where \(\mu\) is the mean value of fracture inclination and \(\kappa\) is the von-Mises distribution coefficient, which indicates the non-uniform distribution of fracture.