# Water Inrush Mechanism of Karst Collapse Column in Coal Seam Floor Based on a Variable Mass Seepage Mechanical Model

Authors: Yang Feng; Lin Hai-xiao; Shen Ling-yan; Yan An-zhi; Wu Yong-jiang
Journal: Geofluids
 (2022)

---

## Abstract

Water inrush of collapse column in the coal seam floor in the karst area is a frequent hydrogeological disaster, and it is becoming one of the significant hidden dangers threatening the safety of coal mine production in China. Based on the self-developed variable mass seepage test system of broken rock mass, different Talbot gradation seepage tests are carried out in the laboratory. A three-flow field coupled nonlinear seepage model is established based on the variable mass seepage theory and fluid mechanics theory of fractured rock mass. The three flow fields are Darcy laminar flow in the aquifer, variable mass seepage in the collapse column, and Navier-Stokes turbulent flow in the roadway. Combined with a project, a numerical calculation model is established to study the water-conducting channel’s formation process and the fluid flow pattern change. The research achievements indicate that the pore pressure gradient and seepage velocity are both 95% fit with the Forchheimer equation, and the seepage of the broken rock mass conforms to the Forchheimer equation. The seepage and migration of particles caused changes in pressure and flow rate, the pressure at the boundary between the aquifer and the collapse column decreased by 0.26 MPa, and the flow rate increased from 0.01 mm/s to 0.064 mm/s. The formation of the water-conducting channel can be divided into two stages, the seepage channel stage and the water-conducting channel stage. The initial fractures continue to develop with the seepage and migration of particles until the water-conducting channel is formed. At the bottom of the collapse column, the Forchheimer number increases from 0.72 to more than 1, and the flow changes from laminar flow to nonlinear seepage. The Forchheimer number at the center of the collapse column fluctuates between 7.5 and 55 at the top of the collapse column. The research results can be treated as an essential basis and reference for the water inrush control of the collapse column and the reasonable prediction of the water inrush amount.


---

## Body

## 1. Introduction

The collapse columns are widely distributed in the northern mining areas of China [1, 2]. The water inrush from the collapse column has the characteristics of adequate concealment, long lag time, and strong suddenness, which causes great harm to the safety of underground engineering [3, 4]. Under the influence of mining and the seepage-erosion action of water, the microcracks of the collapse column gradually expand and penetrate to form a water-conducting channel between the aquifer and the roadway [5, 6]. The seepage of fluid in the collapse column has a nonlinear characteristic [7, 8]. Therefore, the research on the nonlinear seepage model of the collapse column with variable mass has essential theoretical and practical significance for revealing the mechanism of water inrush and the reasonable prediction of water inrush.Currently, the research on the water inrush mechanism and seepage characteristics of the collapse column mainly focuses on theoretical analysis, numerical simulation, and experiments. In terms of theoretical analysis, Yin et al. [9, 10] established the mechanical model of the “thick-walled cylinder,” dividing the collapse column outburst into two modes: top, bottom, and side walls, and established the water inrush criterion of the corresponding mode. Xu et al. [11] used the elastic-plastic theory to analyze the stress distribution state of the collapse column and its surrounding area and judge the possibility of water inrush from the collapse column. In terms of numerical simulation. Zhu et al. [12, 13] took the rock damage evolution as the main line, considered the heterogeneity of rock mass, the coupling effect of seepage and deformation, and proposed the constitutive relationship of the rock failure process under the condition of fluid-solid coupling. Regarding numerical simulation, Yao et al. [14–16] established the variable mass fluid-solid coupling dynamics theory of broken rock mass and got the under different factors. In the aspect of experimental test, Chen et al. [17] conducted a variable mass seepage test of broken rock mass and found that with the action of the confined water pressure at the bottom, a large number of fine particles are lost inside the broken rock mass, resulting in increased porosity and enhanced permeability. Wang et al. [18] considered the influence of water pressure and compaction degree on the water inrush mechanism of karst collapse columns with different buried depths. Zhang et al. [19] conducted a penetration test and found that the collapsed column’s delayed water inrush was caused by the compaction of the broken rock mass caused by mining. Yao et al. [20] comprehensively considered the proportion of broken rock mass, stress state, and other conditions and divided the change of permeability characteristics into three stages.It is evident that the research on the formation process of the water-conducting channel and the change of the fluid flow pattern during the water inrush process of the collapse column is not systematic enough. Therefore, through the broken rock mass seepage test, the hydrodynamic equation of fluid flow in the collapse column is determined, and then, the nonlinear seepage model coupled with three flow fields is established. Combined with the actual engineering, the finite element software COMSOL Multiphysics is used to establish a numerical calculation model to study the water-conducting channel’s formation process and the fluid flow state change. The research results provide guidance for water inrush prevention and water volume prediction of collapse columns.


## 2. Variable Mass Seepage Test of Broken Rock Mass

### 2.1. Test Equipment and Scheme

The variable mass permeability test system of broken rock mass is shown in Figure1. It mainly includes a data acquisition system, water supply, and pressure control system, and variable mass infiltration device considering particle loss.Figure 1
Physical diagram of the variable mass seepage test system for broken rock mass.We select the rock in the collapse column with the possibility of water inrush, making rock particles of different particle sizes. The particle size combination and particle size effect impact the test results [21, 22]. The maximum diameter of the sample cannot exceed 1/5 of the container diameter, and the maximum particle size should be 5 mm. The Talbot continuous gradation formula is used for configuration [23] and is determined by
(1)Pi=diD∗n×100%,wherePi is the proportion of rock particles whose diameter is smaller than di that, %. di is the rock particle diameter, mm. D∗ is the maximum diameter of rock particles, mm. n is the Talbol power exponent.Four samples withn=0.2,0.4,0.6,and0.8 are configured, respectively, 180 g/species, 3 groups of each sample, and the average test results. The mass distribution under different ratios is shown in Table 1.Table 1
Rock particle mass (g) under different TPEs.
Rock particle size (mm)TPE value and corresponding particle mass0.20.40.60.80 ~ 0.2586.0741.1619.689.410.25~0.512.8013.1510.156.980.5 ~ 114.7017.3615.3912.141 ~ 216.8922.8923.3221.142 ~ 526.2441.8750.2253.705 ~ 1 023.3043.5861.2576.62The axial displacement is 4 mm, the confining pressure is 3 MPa, and the seepage pressure is usually 0.2-0.5 of the confining pressure, which is 1 MPa in this study. After determining the test parameters, start the next step.


### 2.2. Test Principle

The flow chart of the variable mass seepage test of broken rock mass is shown in Figure2. Before the test, we record the sample height under the sample stack. During the test, the sample’s height, osmotic pressure, and flow rate at different grades were recorded. After applying confining pressure and axial displacement, according to the drainage method, the volume of water discharged from the sample is affected by the confining pressure the reduced sample volume caused by the confining pressure, which ∆V is collected and recorded by the outlet graduated cylinder.Figure 2
Flow chart of the variable mass seepage test of broken rock mass.The volume of the sample under the combined action of axial pressure and confining pressure before the penetration test,V0′, is calculated by
(2)V0′=πr02h0−Δh−ΔV,wherer0 is the initial radius of the inner sleeve (cm), h0 is the initial height of the broken rock mass sample (cm), Δhis the axial displacement (cm), and ΔV is the discharge volume of permeate under the action of confining pressure (cm3).The porosity of the sample can be calculated by(3)φ0=v0′−v0v0.The equation calculates the cross-sectional seepage area of the sample.(4)A=v0′h.The volume flow through the broken rock mass per unit time isQ; the seepage velocity v under the osmotic pressure is the volume flow to the cross-sectional area.
(5)v=Qπr02−v0/h.The equation can calculate the pore pressure gradient during the seepage stabilization period [24].
(6)Gp=P2−P1h0−Δh=−P1h0−Δh,whereP1 and P2 are the water pressure at the upper and lower ends, respectively. The water pressure at the lower end is the same as the atmospheric pressure, P2=0.


### 2.3. Analysis of Test Results

According to the test results, the scatter diagrams and fitting curves of the four kinds of graded broken rock masses are drawn, as shown in Figure3.Figure 3
Pole pressure gradient and infiltration speed fit curve.
(a)
n=0.2(b)
n=0.4(c)
n=0.6(d)
n=0.8The scatter plot, fitting curve of the pore pressure gradient, and seepage velocity of the broken rock are shown in Figure3. Under different gradations, the R2 value obtained by fitting the pore pressure gradient scatter plot and seepage velocity according to the Forchheimer relation is above 0.96, while the R2 value obtained by fitting Darcy’s law is in the range of 0.65~0.75. The relationship between pore pressure gradient and seepage velocity in the metamorphic seepage process of broken rock under different gradations is more in line with Forchheimer’s law than Darcy’s law [25, 26].


## 2.1. Test Equipment and Scheme

The variable mass permeability test system of broken rock mass is shown in Figure1. It mainly includes a data acquisition system, water supply, and pressure control system, and variable mass infiltration device considering particle loss.Figure 1
Physical diagram of the variable mass seepage test system for broken rock mass.We select the rock in the collapse column with the possibility of water inrush, making rock particles of different particle sizes. The particle size combination and particle size effect impact the test results [21, 22]. The maximum diameter of the sample cannot exceed 1/5 of the container diameter, and the maximum particle size should be 5 mm. The Talbot continuous gradation formula is used for configuration [23] and is determined by
(1)Pi=diD∗n×100%,wherePi is the proportion of rock particles whose diameter is smaller than di that, %. di is the rock particle diameter, mm. D∗ is the maximum diameter of rock particles, mm. n is the Talbol power exponent.Four samples withn=0.2,0.4,0.6,and0.8 are configured, respectively, 180 g/species, 3 groups of each sample, and the average test results. The mass distribution under different ratios is shown in Table 1.Table 1
Rock particle mass (g) under different TPEs.
Rock particle size (mm)TPE value and corresponding particle mass0.20.40.60.80 ~ 0.2586.0741.1619.689.410.25~0.512.8013.1510.156.980.5 ~ 114.7017.3615.3912.141 ~ 216.8922.8923.3221.142 ~ 526.2441.8750.2253.705 ~ 1 023.3043.5861.2576.62The axial displacement is 4 mm, the confining pressure is 3 MPa, and the seepage pressure is usually 0.2-0.5 of the confining pressure, which is 1 MPa in this study. After determining the test parameters, start the next step.


## 2.2. Test Principle

The flow chart of the variable mass seepage test of broken rock mass is shown in Figure2. Before the test, we record the sample height under the sample stack. During the test, the sample’s height, osmotic pressure, and flow rate at different grades were recorded. After applying confining pressure and axial displacement, according to the drainage method, the volume of water discharged from the sample is affected by the confining pressure the reduced sample volume caused by the confining pressure, which ∆V is collected and recorded by the outlet graduated cylinder.Figure 2
Flow chart of the variable mass seepage test of broken rock mass.The volume of the sample under the combined action of axial pressure and confining pressure before the penetration test,V0′, is calculated by
(2)V0′=πr02h0−Δh−ΔV,wherer0 is the initial radius of the inner sleeve (cm), h0 is the initial height of the broken rock mass sample (cm), Δhis the axial displacement (cm), and ΔV is the discharge volume of permeate under the action of confining pressure (cm3).The porosity of the sample can be calculated by(3)φ0=v0′−v0v0.The equation calculates the cross-sectional seepage area of the sample.(4)A=v0′h.The volume flow through the broken rock mass per unit time isQ; the seepage velocity v under the osmotic pressure is the volume flow to the cross-sectional area.
(5)v=Qπr02−v0/h.The equation can calculate the pore pressure gradient during the seepage stabilization period [24].
(6)Gp=P2−P1h0−Δh=−P1h0−Δh,whereP1 and P2 are the water pressure at the upper and lower ends, respectively. The water pressure at the lower end is the same as the atmospheric pressure, P2=0.


## 2.3. Analysis of Test Results

According to the test results, the scatter diagrams and fitting curves of the four kinds of graded broken rock masses are drawn, as shown in Figure3.Figure 3
Pole pressure gradient and infiltration speed fit curve.
(a)
n=0.2(b)
n=0.4(c)
n=0.6(d)
n=0.8The scatter plot, fitting curve of the pore pressure gradient, and seepage velocity of the broken rock are shown in Figure3. Under different gradations, the R2 value obtained by fitting the pore pressure gradient scatter plot and seepage velocity according to the Forchheimer relation is above 0.96, while the R2 value obtained by fitting Darcy’s law is in the range of 0.65~0.75. The relationship between pore pressure gradient and seepage velocity in the metamorphic seepage process of broken rock under different gradations is more in line with Forchheimer’s law than Darcy’s law [25, 26].


## 3. Three-Flow Field Coupled Nonlinear Seepage Model

### 3.1. Darcy Laminar Flow in Aquifer

The groundwater moves in the confined aquifer, the general movement speed is small, and the fluid velocity and pressure gradient satisfy the linear Darcy laminar flow. The Darcy equations are represented by(7)∇⋅ρu=Qm,(8)q=−kμ∇p,wherek is the permeability, p is the fluid pressure, ρ is the fluid density, Qm is the source-sink term, q is the Darcy velocity, and μ is the dynamic viscosity.


### 3.2. Variable Mass Seepage of Collapse Column

The collapsed column has three parts: rock mass skeleton, fluid medium, and fine particles. The variable mass seepage of the collapse column is based on the following assumptions:(1)
Collapse column is regarded as the dual medium of pore and fissure(2)
Suspended particles are incompressible and move at the same speed as the fluid(3)
The fluid flow in the collapse column is not affected by filler particles(4)
The influence of stress on the deformation of the collapse column is ignoredThe collapsed column is affected by mining, resulting in random fractures. The micro-element is composed of rock matrix and fractures. The model is shown in Figure4.Figure 4
Schematic diagram of microelement.The fracture opening isa, the matrix block length is s, s>>a, the microelement fracture rate is 3·a/s, and the matrix porosity is set as φ. The microporosity is obtained by
(9)nf=φ+3as.


#### 3.2.1. Forchheimer Equation

The fluid movement in the fractured rock with different gradations conforms to the Forchheimer equation [27, 28]. Therefore, the Brinkman equation is used as the fluid motion equation in the collapse column. The finite element software COMSOL Multiphysics is described by Forchheimer’s improved Brinkman equation, and the expression is shown in
(10)qμk+βfq+Qbrn=∇−pI+μn∇q+∇qT+F,(11)ρ∇⋅q=Qbr,whereβf is the non-Darcy factor, n is the porosity, F is thevolume force affecting the fluid, and I is the unit matrix.According to the literature [29], the non-Darcy coefficient is obtained by empirical
(12)βf=0.005n5.5k.The relationship between the Forchheimer and non-Darcy coefficient satisfies(13)cf=βfk.


#### 3.2.2. Mass Conservation Equation

For dual-medium variable mass seepage, the diffusion of particles is neglected, and the mass conservation equation of particles in the collapse column satisfies equation (14).Conservation of fluid mass satisfies:(14)3s∂ac∂t+∂cφ∂t+∇⋅cq¯=3s∂a∂t+∂φ∂t,(15)3s∂a1−c∂t+∂1−cφ∂t+∇⋅1−cq¯=0,whereC is the particle volume concentration, q¯ is the seepage velocity, and t is the seepage time.


#### 3.2.3. Evolution Equation of Porosity and Fracture Rate

According to the literature [30], the evolution of unit matrix porosity and fracture rate under erosion satisfies
(16)∂φ∂t=λ1ρsφmax−φCq,(17)∂a∂t=λ2ρsamax−aCq,whereλ1 is the pore dissolution coefficient and λ2 is the fracture dissolution coefficient.


#### 3.2.4. Auxiliary Equations

In the process of variable mass seepage, particle seepage migration causes the porosity and permeability of the collapse column to change. After the formation of the aqueduct in the collapse column, water inrush occurs, and the porosity and permeability are stable. The change in permeability can be described by the cubic law, which is expressed as(18)k=k0φφ031−φ01−φ2+a312s.According to the literature [31], the relationship between the particle-containing fluids’ dynamic viscosity and the water’s dynamic viscosity satisfies the
(19)ηC=1−C−2.5.


### 3.3. Navier-Stokes Turbulence in Roadway Face

After the collapse column water inrush, the water flows into the roadway, and the water velocity and flow rate are large and unstable. Therefore, Navier-Stokes turbulence equations (20) and (21) describe the fluid state in the roadway.
(20)ρq⋅∇q=∇⋅−pI+μ∇q+∇qT+F,(21)ρ∇⋅q=0.


### 3.4. Transition Boundary Conditions of Flow Field

During the seepage process of the aquifer and the collapse column, the flow rate and pressure of the fluid change dynamically with time. Therefore, the two adjacent flow fields simultaneously satisfy the same flow rate and pressure. Equations (22) and (23) are satisfied at the aquifer and collapse column boundary.
(22)pdl=pbr,(23)ubr=udl.Equations (24) and (25) are satisfied at the boundary between the collapse column and the roadway.
(24)pbr=pNS,(25)ubr=uNS.The simplified model of the three-flow field is shown in Figure5. Combined with the above equations and the actual project, the COMSOL Multiphysics finite element software is employed to establish a numerical model of the coupled nonlinear seepage flow in the three-flow field for numerical calculation and analysis.Figure 5
Simplified model of the three-flow field.


## 3.1. Darcy Laminar Flow in Aquifer

The groundwater moves in the confined aquifer, the general movement speed is small, and the fluid velocity and pressure gradient satisfy the linear Darcy laminar flow. The Darcy equations are represented by(7)∇⋅ρu=Qm,(8)q=−kμ∇p,wherek is the permeability, p is the fluid pressure, ρ is the fluid density, Qm is the source-sink term, q is the Darcy velocity, and μ is the dynamic viscosity.


## 3.2. Variable Mass Seepage of Collapse Column

The collapsed column has three parts: rock mass skeleton, fluid medium, and fine particles. The variable mass seepage of the collapse column is based on the following assumptions:(1)
Collapse column is regarded as the dual medium of pore and fissure(2)
Suspended particles are incompressible and move at the same speed as the fluid(3)
The fluid flow in the collapse column is not affected by filler particles(4)
The influence of stress on the deformation of the collapse column is ignoredThe collapsed column is affected by mining, resulting in random fractures. The micro-element is composed of rock matrix and fractures. The model is shown in Figure4.Figure 4
Schematic diagram of microelement.The fracture opening isa, the matrix block length is s, s>>a, the microelement fracture rate is 3·a/s, and the matrix porosity is set as φ. The microporosity is obtained by
(9)nf=φ+3as.


### 3.2.1. Forchheimer Equation

The fluid movement in the fractured rock with different gradations conforms to the Forchheimer equation [27, 28]. Therefore, the Brinkman equation is used as the fluid motion equation in the collapse column. The finite element software COMSOL Multiphysics is described by Forchheimer’s improved Brinkman equation, and the expression is shown in
(10)qμk+βfq+Qbrn=∇−pI+μn∇q+∇qT+F,(11)ρ∇⋅q=Qbr,whereβf is the non-Darcy factor, n is the porosity, F is thevolume force affecting the fluid, and I is the unit matrix.According to the literature [29], the non-Darcy coefficient is obtained by empirical
(12)βf=0.005n5.5k.The relationship between the Forchheimer and non-Darcy coefficient satisfies(13)cf=βfk.


### 3.2.2. Mass Conservation Equation

For dual-medium variable mass seepage, the diffusion of particles is neglected, and the mass conservation equation of particles in the collapse column satisfies equation (14).Conservation of fluid mass satisfies:(14)3s∂ac∂t+∂cφ∂t+∇⋅cq¯=3s∂a∂t+∂φ∂t,(15)3s∂a1−c∂t+∂1−cφ∂t+∇⋅1−cq¯=0,whereC is the particle volume concentration, q¯ is the seepage velocity, and t is the seepage time.


### 3.2.3. Evolution Equation of Porosity and Fracture Rate

According to the literature [30], the evolution of unit matrix porosity and fracture rate under erosion satisfies
(16)∂φ∂t=λ1ρsφmax−φCq,(17)∂a∂t=λ2ρsamax−aCq,whereλ1 is the pore dissolution coefficient and λ2 is the fracture dissolution coefficient.


### 3.2.4. Auxiliary Equations

In the process of variable mass seepage, particle seepage migration causes the porosity and permeability of the collapse column to change. After the formation of the aqueduct in the collapse column, water inrush occurs, and the porosity and permeability are stable. The change in permeability can be described by the cubic law, which is expressed as(18)k=k0φφ031−φ01−φ2+a312s.According to the literature [31], the relationship between the particle-containing fluids’ dynamic viscosity and the water’s dynamic viscosity satisfies the
(19)ηC=1−C−2.5.


## 3.2.1. Forchheimer Equation

The fluid movement in the fractured rock with different gradations conforms to the Forchheimer equation [27, 28]. Therefore, the Brinkman equation is used as the fluid motion equation in the collapse column. The finite element software COMSOL Multiphysics is described by Forchheimer’s improved Brinkman equation, and the expression is shown in
(10)qμk+βfq+Qbrn=∇−pI+μn∇q+∇qT+F,(11)ρ∇⋅q=Qbr,whereβf is the non-Darcy factor, n is the porosity, F is thevolume force affecting the fluid, and I is the unit matrix.According to the literature [29], the non-Darcy coefficient is obtained by empirical
(12)βf=0.005n5.5k.The relationship between the Forchheimer and non-Darcy coefficient satisfies(13)cf=βfk.


## 3.2.2. Mass Conservation Equation

For dual-medium variable mass seepage, the diffusion of particles is neglected, and the mass conservation equation of particles in the collapse column satisfies equation (14).Conservation of fluid mass satisfies:(14)3s∂ac∂t+∂cφ∂t+∇⋅cq¯=3s∂a∂t+∂φ∂t,(15)3s∂a1−c∂t+∂1−cφ∂t+∇⋅1−cq¯=0,whereC is the particle volume concentration, q¯ is the seepage velocity, and t is the seepage time.


## 3.2.3. Evolution Equation of Porosity and Fracture Rate

According to the literature [30], the evolution of unit matrix porosity and fracture rate under erosion satisfies
(16)∂φ∂t=λ1ρsφmax−φCq,(17)∂a∂t=λ2ρsamax−aCq,whereλ1 is the pore dissolution coefficient and λ2 is the fracture dissolution coefficient.


## 3.2.4. Auxiliary Equations

In the process of variable mass seepage, particle seepage migration causes the porosity and permeability of the collapse column to change. After the formation of the aqueduct in the collapse column, water inrush occurs, and the porosity and permeability are stable. The change in permeability can be described by the cubic law, which is expressed as(18)k=k0φφ031−φ01−φ2+a312s.According to the literature [31], the relationship between the particle-containing fluids’ dynamic viscosity and the water’s dynamic viscosity satisfies the
(19)ηC=1−C−2.5.


## 3.3. Navier-Stokes Turbulence in Roadway Face

After the collapse column water inrush, the water flows into the roadway, and the water velocity and flow rate are large and unstable. Therefore, Navier-Stokes turbulence equations (20) and (21) describe the fluid state in the roadway.
(20)ρq⋅∇q=∇⋅−pI+μ∇q+∇qT+F,(21)ρ∇⋅q=0.


## 3.4. Transition Boundary Conditions of Flow Field

During the seepage process of the aquifer and the collapse column, the flow rate and pressure of the fluid change dynamically with time. Therefore, the two adjacent flow fields simultaneously satisfy the same flow rate and pressure. Equations (22) and (23) are satisfied at the aquifer and collapse column boundary.
(22)pdl=pbr,(23)ubr=udl.Equations (24) and (25) are satisfied at the boundary between the collapse column and the roadway.
(24)pbr=pNS,(25)ubr=uNS.The simplified model of the three-flow field is shown in Figure5. Combined with the above equations and the actual project, the COMSOL Multiphysics finite element software is employed to establish a numerical model of the coupled nonlinear seepage flow in the three-flow field for numerical calculation and analysis.Figure 5
Simplified model of the three-flow field.


## 4. Numerical Simulation

### 4.1. Project Overview

At 7:29 on March 1, 2010, the water inrush at the excavation face of the 16 coal +870 horizontal return air roadway suddenly increased. By 8:40, the water inrush amounted to7.6×104 m3, and the average water inflow during this period was 6.5×104 m3/h. Because the water inrush is far greater than the drainage capacity of the mine, the mine is flooded. According to the research, the net sectional area of 16 coal return air roadways is 17.8 m2, the net width is 5 m, and the net height is 3.56 m. The columnar water-conducting channel is a developed small austenitic lime water-conducting karst collapse column.


### 4.2. Numerical Model

According to the geological structure of No. 16 coal-water inrush and the characteristics of the confined aquifer, the model is simplified, as shown in Figure6. According to the water inrush flow path, the model can be divided into three parts, 40m×10m floor limestone aquifer, 8m×40m cylindrical collapse column, and the 20m×3.66m roadway. The boundary conditions of the model are set. The two sides of the aquifer are fixed pressure boundaries, p=2.1×106 Pa, and the roadway outlet is set as a fixed pressure boundary, p=0.1×106 Pa. The boundary pressure and velocity of adjacent flow fields are equal. The effect of mining on the collapse column was quantified using the Wei bull uniformity index m=3 and the crack width 1.25×10−4 m. Set up monitoring lines in aquifers, collapse columns, and roadways. Set the monitoring line M-N to analyze the changes in flow velocity and pressure at the boundary between the collapse column and the roadway. The mechanical parameters of the collapse column are shown in Table 2, and the physical and mechanical parameters of the aquifer, collapse column, and roadway are shown in Table 3.Figure 6
Schematic diagram of the numerical model.Table 2
Mechanical parameters of collapse column.
Physical parametersFall columnUnitDensityρ2260kg/m3Elastic modulusE1.0GPaPoisson’s ratioμ0.3—Matrix block length0.02mFissure width1.25e−4mPore erosion coefficientλ10.01m-1Fissure erosion coefficientλ20.01m-1Initial suspended particle concentration in pores and fissures0.01—Initial permeabilityk05×10−13m2Table 3
Hydrodynamic parameters of the three stages.
Physical parametersAquiferFall columnCoal layer laneUnitDensityρ100010001000Kg/m3Dynamic viscosity1.0e−3Formula (24)1.0e-3Pa·sAcceleration factor—1.0——Porosity0.14φ—mPenetration2.4e−11k—mNon-Darcy factor—Formula (12)—m2


### 4.3. Analysis of Simulation Results

Selecting six-time points (0 s, 4000 s, 8000 s, 12000 s, 16000 s, and 20000 s) as monitoring time points, three points from bottom to top in the collapse column (21.5, 10.5), (21.5, 30), (21.5, 44.5) for monitoring were used.


#### 4.3.1. Pressure Cloud Map

It is evident that with the increase of time, the pressure of the aquifer on the monitoring line and the pressure of the collapse column gradually decrease in Figures7 and 8. The pressure is stable at about 2.0 MPa, while the pressure in the collapse column drops rapidly. When t=1.6×104 s, the pressure changes significantly compared with the previous one and drops to 1.92 MPa, and the pressure drop rate in the collapse column somewhat decreases, indicating that the water-conducting channel is being formed. When t=2.0×104 s, the pressure at the junction of the aquifer and the collapse column drops to 1.74 MPa, which is 12.25 times that at t=6×103 s. The results show that the seepage migration of mass particles can relieve the pressure of the aquifer. The pressure drops slowly and then drops rapidly.Figure 7
Fluid pressure distribution at different time points.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sFigure 8
Monitoring pressure on line A-B-C-D-E-F.


#### 4.3.2. Speed Graph

The fluid velocity on the monitoring line changes dynamically with time in Figures9 and 10. Before t=1.2×104 s, the velocity of the aquifer and the collapse column grows slowly. When t=4.0×104 s, the maximum speed of the roadway is 0.46 m/s. When t=8×103 s, the speed is reduced and relatively, stable at 9.8×10−4 m/s. When t=1.2×104 s, the velocity increases again. The reason is that aquifer pressure drops gradually, fluid velocity increases, and particle seepage-migration speeds up. When t=1.6×104 s, the velocity in the aquifer and the collapse column increase. When t=2.0×104 s, the velocity of the aquifer and the collapse column increased significantly. The results show that the mass loss of the collapse column causes the increase of the fluid velocity, and the fluid velocity in the roadway fluctuates wildly, which satisfies the turbulent flow characteristics. The heterogeneity of the collapse column causes the fluid velocity to fluctuate on the detection line, and the overall trend is increasing.Figure 9
Monitoring speed on line A-B-C-D.Figure 10
Monitoring speed on line D-E-F.


#### 4.3.3. The Formation Process of the Water Channel

The cloud diagram of the change of the crack opening at different times is shown in Figure11. The maximum crack width does not change significantly before t=1.2×104 s, while the minimum crack is opening increases from 9.98×10−6 m to 1.25e−5 m. When t=1.6×104 s, two discontinuous seepage channels are formed in the upper part and the top of the fissure cloud. There is no connection between the channels and no connection between the channel and the top boundary. The maximum fissure width is 4.06×10−4 m, 1.58 times that at t=1.2×104 s. When the time is 2.0×104 s, the seepage channels are connected to form a water channel, and a water inrush occurs in the collapse column. The maximum width of the crack is 6.83×10−4 m.Figure 11
Cloud map of fracture opening changes at different times.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sThe boundary lineM-N velocity curve for the boundary velocity is shown in Figure 12. When time is 1.2×104 s, the aqueduct is not formed, and the boundary velocity increases slowly. When time is 1.6×104 s, the boundary velocity increases significantly, and the reason is that several discontinuous seepage channels are formed inside the collapse column, which causes the increase of fluid velocity. When time is 2.0×104 s, the fluid velocity increases suddenly at the 2.5 m position of the boundary, and the results show that the water inrush occurs in the collapse column.Figure 12
Speed on boundary monitoring lineM-N.The graph of outflow water volume at different times of the collapse column is shown in Figure13. Before the 1.2×104 s, the fluid flow increased slowly, reaching 59.4 m3/h. When time is 1.6×104 s, the seepage channel gradually forms, and the water inflow reaches 89.7 m3/h. When the time was 2.0×104 s, the water channel was formed, and the water inflow increased sharply to 176.8 m3/h. The research results show that after the formation of the aqueduct in the collapse column, water inrush occurs, and the fluid velocity and flow show a slow to rapid increase.Figure 13
The amount of water inflow from the collapse column at different times.


#### 4.3.4. Nonlinear Seepage of Collapse Column

For the judgment of the degree of nonlinearity, the Forchheimer number is commonly used, which is the ratio of the fluid inertia term to the first-order term in the Forchheimer equation [32].
(26)F0=βfkρvfμ,whereβf non-Darcy factor, calculated by equation (12), k is permeability, ρ is density, and μ is dynamic viscosity, calculated by equation (13).The permeability curves of three points in the collapse column at different times are shown in Figure14. Before the 1.2×104 s, the permeability of domain points 1, 2, and 3 increased slowly. When time is 1.6×104 s, the permeability of domains 1, 2, and 3 continues to increase. There is a sudden increase in point 3, and the permeability changed from 1.88×10−11 m2 to 7.88×10−11 m2. At time 2.0×104 s, the changing trend of permeability is the same as that at time 1.6×104 s, and the permeability at point 3 increases to 6.64×10−10 m2. It shows that the seepage migration of mass particles gradually increases the permeability of the collapse column. After the seepage channel is formed, mass loss occurs, and the permeability gradually increases. The permeability increases first at the surface position. The variation of porosity shows a similar law to that of permeability. The change in void ratio is shown in Figure 15, which is the same as the changing permeability trend.Figure 14
Domain penetration rate at different times.Figure 15
Void ratio of different time domain points.The change curves of hydrodynamic viscosity of three domain points in the model at different times are shown in Figure16. Before the time of 1.6×104 s, the hydrodynamic viscosity increases from the initial 1.03 MPa·s to 1.12 MPa·s at domains 1 and 2. The viscosity of domain 3 increased from 1.03 MPa·s to 1.28 MPa·s, 2.8 times that of domains 1 and 2. The domain point 3 increases from 1.03 MPa·s to 1.28 MPa·s. The reason is that the fluid carries certain particles, which cause the change of dynamic viscosity. On the other hand, the continuous movement and accumulation result in a large volume concentration of particles at the top and a rapid increase in the hydrodynamic viscosity. At time 2.0×104 s, the hydrodynamic viscosity of domain point 3 increased to 1.4 MPa·s, while the hydrodynamic viscosity of domain point 1 and domain point 2 decreased. The reason is that the particles at the bottom and middle of the collapse column continued to flow water. The upward movement causes a decrease in the hydrodynamic viscosity, while the top carries a large number of particles, causing an increase in the dynamic viscosity.Figure 16
Dynamic viscosity at different time domains.The Forchheimer number at different times for the three-domain points in the collapsed column is shown in Figure17. Before the time of 1.2×104 , the bottom of the collapse column is F0<1, the viscous resistance plays a dominant role in the seepage process of water, the flow state tends to be linear, and the degree of nonlinearity is weak. With time, F0 increases, and the degree of nonlinearity increases. For the viscous resistance F0>1 in the middle position, the inertial resistance plays a leading role and the degree of nonlinearity increases. It is stronger than the bottom of the collapse column, and F0 shows an increasing trend, indicating that with the loss of mass, the nonlinearity of the flow state of the water in the middle increases. At time is 1.6×104 s, the position F0 at the bottom of the collapse column is around 1 and drops to 1 at t=20000 s. In contrast, the middle position drops and rises abruptly at t=1.6×104 s and t=2.0×104 s, respectively, which indicates that the mass loss in the middle blocks the channel and leads to nonlinear weakening. The nonlinearity shows the opposite law with the bottom and middle of the collapse column.Figure 17
Forchheimer numbers of domain points at different times.


## 4.1. Project Overview

At 7:29 on March 1, 2010, the water inrush at the excavation face of the 16 coal +870 horizontal return air roadway suddenly increased. By 8:40, the water inrush amounted to7.6×104 m3, and the average water inflow during this period was 6.5×104 m3/h. Because the water inrush is far greater than the drainage capacity of the mine, the mine is flooded. According to the research, the net sectional area of 16 coal return air roadways is 17.8 m2, the net width is 5 m, and the net height is 3.56 m. The columnar water-conducting channel is a developed small austenitic lime water-conducting karst collapse column.


## 4.2. Numerical Model

According to the geological structure of No. 16 coal-water inrush and the characteristics of the confined aquifer, the model is simplified, as shown in Figure6. According to the water inrush flow path, the model can be divided into three parts, 40m×10m floor limestone aquifer, 8m×40m cylindrical collapse column, and the 20m×3.66m roadway. The boundary conditions of the model are set. The two sides of the aquifer are fixed pressure boundaries, p=2.1×106 Pa, and the roadway outlet is set as a fixed pressure boundary, p=0.1×106 Pa. The boundary pressure and velocity of adjacent flow fields are equal. The effect of mining on the collapse column was quantified using the Wei bull uniformity index m=3 and the crack width 1.25×10−4 m. Set up monitoring lines in aquifers, collapse columns, and roadways. Set the monitoring line M-N to analyze the changes in flow velocity and pressure at the boundary between the collapse column and the roadway. The mechanical parameters of the collapse column are shown in Table 2, and the physical and mechanical parameters of the aquifer, collapse column, and roadway are shown in Table 3.Figure 6
Schematic diagram of the numerical model.Table 2
Mechanical parameters of collapse column.
Physical parametersFall columnUnitDensityρ2260kg/m3Elastic modulusE1.0GPaPoisson’s ratioμ0.3—Matrix block length0.02mFissure width1.25e−4mPore erosion coefficientλ10.01m-1Fissure erosion coefficientλ20.01m-1Initial suspended particle concentration in pores and fissures0.01—Initial permeabilityk05×10−13m2Table 3
Hydrodynamic parameters of the three stages.
Physical parametersAquiferFall columnCoal layer laneUnitDensityρ100010001000Kg/m3Dynamic viscosity1.0e−3Formula (24)1.0e-3Pa·sAcceleration factor—1.0——Porosity0.14φ—mPenetration2.4e−11k—mNon-Darcy factor—Formula (12)—m2


## 4.3. Analysis of Simulation Results

Selecting six-time points (0 s, 4000 s, 8000 s, 12000 s, 16000 s, and 20000 s) as monitoring time points, three points from bottom to top in the collapse column (21.5, 10.5), (21.5, 30), (21.5, 44.5) for monitoring were used.


### 4.3.1. Pressure Cloud Map

It is evident that with the increase of time, the pressure of the aquifer on the monitoring line and the pressure of the collapse column gradually decrease in Figures7 and 8. The pressure is stable at about 2.0 MPa, while the pressure in the collapse column drops rapidly. When t=1.6×104 s, the pressure changes significantly compared with the previous one and drops to 1.92 MPa, and the pressure drop rate in the collapse column somewhat decreases, indicating that the water-conducting channel is being formed. When t=2.0×104 s, the pressure at the junction of the aquifer and the collapse column drops to 1.74 MPa, which is 12.25 times that at t=6×103 s. The results show that the seepage migration of mass particles can relieve the pressure of the aquifer. The pressure drops slowly and then drops rapidly.Figure 7
Fluid pressure distribution at different time points.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sFigure 8
Monitoring pressure on line A-B-C-D-E-F.


### 4.3.2. Speed Graph

The fluid velocity on the monitoring line changes dynamically with time in Figures9 and 10. Before t=1.2×104 s, the velocity of the aquifer and the collapse column grows slowly. When t=4.0×104 s, the maximum speed of the roadway is 0.46 m/s. When t=8×103 s, the speed is reduced and relatively, stable at 9.8×10−4 m/s. When t=1.2×104 s, the velocity increases again. The reason is that aquifer pressure drops gradually, fluid velocity increases, and particle seepage-migration speeds up. When t=1.6×104 s, the velocity in the aquifer and the collapse column increase. When t=2.0×104 s, the velocity of the aquifer and the collapse column increased significantly. The results show that the mass loss of the collapse column causes the increase of the fluid velocity, and the fluid velocity in the roadway fluctuates wildly, which satisfies the turbulent flow characteristics. The heterogeneity of the collapse column causes the fluid velocity to fluctuate on the detection line, and the overall trend is increasing.Figure 9
Monitoring speed on line A-B-C-D.Figure 10
Monitoring speed on line D-E-F.


### 4.3.3. The Formation Process of the Water Channel

The cloud diagram of the change of the crack opening at different times is shown in Figure11. The maximum crack width does not change significantly before t=1.2×104 s, while the minimum crack is opening increases from 9.98×10−6 m to 1.25e−5 m. When t=1.6×104 s, two discontinuous seepage channels are formed in the upper part and the top of the fissure cloud. There is no connection between the channels and no connection between the channel and the top boundary. The maximum fissure width is 4.06×10−4 m, 1.58 times that at t=1.2×104 s. When the time is 2.0×104 s, the seepage channels are connected to form a water channel, and a water inrush occurs in the collapse column. The maximum width of the crack is 6.83×10−4 m.Figure 11
Cloud map of fracture opening changes at different times.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sThe boundary lineM-N velocity curve for the boundary velocity is shown in Figure 12. When time is 1.2×104 s, the aqueduct is not formed, and the boundary velocity increases slowly. When time is 1.6×104 s, the boundary velocity increases significantly, and the reason is that several discontinuous seepage channels are formed inside the collapse column, which causes the increase of fluid velocity. When time is 2.0×104 s, the fluid velocity increases suddenly at the 2.5 m position of the boundary, and the results show that the water inrush occurs in the collapse column.Figure 12
Speed on boundary monitoring lineM-N.The graph of outflow water volume at different times of the collapse column is shown in Figure13. Before the 1.2×104 s, the fluid flow increased slowly, reaching 59.4 m3/h. When time is 1.6×104 s, the seepage channel gradually forms, and the water inflow reaches 89.7 m3/h. When the time was 2.0×104 s, the water channel was formed, and the water inflow increased sharply to 176.8 m3/h. The research results show that after the formation of the aqueduct in the collapse column, water inrush occurs, and the fluid velocity and flow show a slow to rapid increase.Figure 13
The amount of water inflow from the collapse column at different times.


### 4.3.4. Nonlinear Seepage of Collapse Column

For the judgment of the degree of nonlinearity, the Forchheimer number is commonly used, which is the ratio of the fluid inertia term to the first-order term in the Forchheimer equation [32].
(26)F0=βfkρvfμ,whereβf non-Darcy factor, calculated by equation (12), k is permeability, ρ is density, and μ is dynamic viscosity, calculated by equation (13).The permeability curves of three points in the collapse column at different times are shown in Figure14. Before the 1.2×104 s, the permeability of domain points 1, 2, and 3 increased slowly. When time is 1.6×104 s, the permeability of domains 1, 2, and 3 continues to increase. There is a sudden increase in point 3, and the permeability changed from 1.88×10−11 m2 to 7.88×10−11 m2. At time 2.0×104 s, the changing trend of permeability is the same as that at time 1.6×104 s, and the permeability at point 3 increases to 6.64×10−10 m2. It shows that the seepage migration of mass particles gradually increases the permeability of the collapse column. After the seepage channel is formed, mass loss occurs, and the permeability gradually increases. The permeability increases first at the surface position. The variation of porosity shows a similar law to that of permeability. The change in void ratio is shown in Figure 15, which is the same as the changing permeability trend.Figure 14
Domain penetration rate at different times.Figure 15
Void ratio of different time domain points.The change curves of hydrodynamic viscosity of three domain points in the model at different times are shown in Figure16. Before the time of 1.6×104 s, the hydrodynamic viscosity increases from the initial 1.03 MPa·s to 1.12 MPa·s at domains 1 and 2. The viscosity of domain 3 increased from 1.03 MPa·s to 1.28 MPa·s, 2.8 times that of domains 1 and 2. The domain point 3 increases from 1.03 MPa·s to 1.28 MPa·s. The reason is that the fluid carries certain particles, which cause the change of dynamic viscosity. On the other hand, the continuous movement and accumulation result in a large volume concentration of particles at the top and a rapid increase in the hydrodynamic viscosity. At time 2.0×104 s, the hydrodynamic viscosity of domain point 3 increased to 1.4 MPa·s, while the hydrodynamic viscosity of domain point 1 and domain point 2 decreased. The reason is that the particles at the bottom and middle of the collapse column continued to flow water. The upward movement causes a decrease in the hydrodynamic viscosity, while the top carries a large number of particles, causing an increase in the dynamic viscosity.Figure 16
Dynamic viscosity at different time domains.The Forchheimer number at different times for the three-domain points in the collapsed column is shown in Figure17. Before the time of 1.2×104 , the bottom of the collapse column is F0<1, the viscous resistance plays a dominant role in the seepage process of water, the flow state tends to be linear, and the degree of nonlinearity is weak. With time, F0 increases, and the degree of nonlinearity increases. For the viscous resistance F0>1 in the middle position, the inertial resistance plays a leading role and the degree of nonlinearity increases. It is stronger than the bottom of the collapse column, and F0 shows an increasing trend, indicating that with the loss of mass, the nonlinearity of the flow state of the water in the middle increases. At time is 1.6×104 s, the position F0 at the bottom of the collapse column is around 1 and drops to 1 at t=20000 s. In contrast, the middle position drops and rises abruptly at t=1.6×104 s and t=2.0×104 s, respectively, which indicates that the mass loss in the middle blocks the channel and leads to nonlinear weakening. The nonlinearity shows the opposite law with the bottom and middle of the collapse column.Figure 17
Forchheimer numbers of domain points at different times.


## 4.3.1. Pressure Cloud Map

It is evident that with the increase of time, the pressure of the aquifer on the monitoring line and the pressure of the collapse column gradually decrease in Figures7 and 8. The pressure is stable at about 2.0 MPa, while the pressure in the collapse column drops rapidly. When t=1.6×104 s, the pressure changes significantly compared with the previous one and drops to 1.92 MPa, and the pressure drop rate in the collapse column somewhat decreases, indicating that the water-conducting channel is being formed. When t=2.0×104 s, the pressure at the junction of the aquifer and the collapse column drops to 1.74 MPa, which is 12.25 times that at t=6×103 s. The results show that the seepage migration of mass particles can relieve the pressure of the aquifer. The pressure drops slowly and then drops rapidly.Figure 7
Fluid pressure distribution at different time points.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sFigure 8
Monitoring pressure on line A-B-C-D-E-F.


## 4.3.2. Speed Graph

The fluid velocity on the monitoring line changes dynamically with time in Figures9 and 10. Before t=1.2×104 s, the velocity of the aquifer and the collapse column grows slowly. When t=4.0×104 s, the maximum speed of the roadway is 0.46 m/s. When t=8×103 s, the speed is reduced and relatively, stable at 9.8×10−4 m/s. When t=1.2×104 s, the velocity increases again. The reason is that aquifer pressure drops gradually, fluid velocity increases, and particle seepage-migration speeds up. When t=1.6×104 s, the velocity in the aquifer and the collapse column increase. When t=2.0×104 s, the velocity of the aquifer and the collapse column increased significantly. The results show that the mass loss of the collapse column causes the increase of the fluid velocity, and the fluid velocity in the roadway fluctuates wildly, which satisfies the turbulent flow characteristics. The heterogeneity of the collapse column causes the fluid velocity to fluctuate on the detection line, and the overall trend is increasing.Figure 9
Monitoring speed on line A-B-C-D.Figure 10
Monitoring speed on line D-E-F.


## 4.3.3. The Formation Process of the Water Channel

The cloud diagram of the change of the crack opening at different times is shown in Figure11. The maximum crack width does not change significantly before t=1.2×104 s, while the minimum crack is opening increases from 9.98×10−6 m to 1.25e−5 m. When t=1.6×104 s, two discontinuous seepage channels are formed in the upper part and the top of the fissure cloud. There is no connection between the channels and no connection between the channel and the top boundary. The maximum fissure width is 4.06×10−4 m, 1.58 times that at t=1.2×104 s. When the time is 2.0×104 s, the seepage channels are connected to form a water channel, and a water inrush occurs in the collapse column. The maximum width of the crack is 6.83×10−4 m.Figure 11
Cloud map of fracture opening changes at different times.
(a)
0 s(b)
4000 s(c)
8000 s(d)
12000 s(e)
16000 s(f)
2000 sThe boundary lineM-N velocity curve for the boundary velocity is shown in Figure 12. When time is 1.2×104 s, the aqueduct is not formed, and the boundary velocity increases slowly. When time is 1.6×104 s, the boundary velocity increases significantly, and the reason is that several discontinuous seepage channels are formed inside the collapse column, which causes the increase of fluid velocity. When time is 2.0×104 s, the fluid velocity increases suddenly at the 2.5 m position of the boundary, and the results show that the water inrush occurs in the collapse column.Figure 12
Speed on boundary monitoring lineM-N.The graph of outflow water volume at different times of the collapse column is shown in Figure13. Before the 1.2×104 s, the fluid flow increased slowly, reaching 59.4 m3/h. When time is 1.6×104 s, the seepage channel gradually forms, and the water inflow reaches 89.7 m3/h. When the time was 2.0×104 s, the water channel was formed, and the water inflow increased sharply to 176.8 m3/h. The research results show that after the formation of the aqueduct in the collapse column, water inrush occurs, and the fluid velocity and flow show a slow to rapid increase.Figure 13
The amount of water inflow from the collapse column at different times.


## 4.3.4. Nonlinear Seepage of Collapse Column

For the judgment of the degree of nonlinearity, the Forchheimer number is commonly used, which is the ratio of the fluid inertia term to the first-order term in the Forchheimer equation [32].
(26)F0=βfkρvfμ,whereβf non-Darcy factor, calculated by equation (12), k is permeability, ρ is density, and μ is dynamic viscosity, calculated by equation (13).The permeability curves of three points in the collapse column at different times are shown in Figure14. Before the 1.2×104 s, the permeability of domain points 1, 2, and 3 increased slowly. When time is 1.6×104 s, the permeability of domains 1, 2, and 3 continues to increase. There is a sudden increase in point 3, and the permeability changed from 1.88×10−11 m2 to 7.88×10−11 m2. At time 2.0×104 s, the changing trend of permeability is the same as that at time 1.6×104 s, and the permeability at point 3 increases to 6.64×10−10 m2. It shows that the seepage migration of mass particles gradually increases the permeability of the collapse column. After the seepage channel is formed, mass loss occurs, and the permeability gradually increases. The permeability increases first at the surface position. The variation of porosity shows a similar law to that of permeability. The change in void ratio is shown in Figure 15, which is the same as the changing permeability trend.Figure 14
Domain penetration rate at different times.Figure 15
Void ratio of different time domain points.The change curves of hydrodynamic viscosity of three domain points in the model at different times are shown in Figure16. Before the time of 1.6×104 s, the hydrodynamic viscosity increases from the initial 1.03 MPa·s to 1.12 MPa·s at domains 1 and 2. The viscosity of domain 3 increased from 1.03 MPa·s to 1.28 MPa·s, 2.8 times that of domains 1 and 2. The domain point 3 increases from 1.03 MPa·s to 1.28 MPa·s. The reason is that the fluid carries certain particles, which cause the change of dynamic viscosity. On the other hand, the continuous movement and accumulation result in a large volume concentration of particles at the top and a rapid increase in the hydrodynamic viscosity. At time 2.0×104 s, the hydrodynamic viscosity of domain point 3 increased to 1.4 MPa·s, while the hydrodynamic viscosity of domain point 1 and domain point 2 decreased. The reason is that the particles at the bottom and middle of the collapse column continued to flow water. The upward movement causes a decrease in the hydrodynamic viscosity, while the top carries a large number of particles, causing an increase in the dynamic viscosity.Figure 16
Dynamic viscosity at different time domains.The Forchheimer number at different times for the three-domain points in the collapsed column is shown in Figure17. Before the time of 1.2×104 , the bottom of the collapse column is F0<1, the viscous resistance plays a dominant role in the seepage process of water, the flow state tends to be linear, and the degree of nonlinearity is weak. With time, F0 increases, and the degree of nonlinearity increases. For the viscous resistance F0>1 in the middle position, the inertial resistance plays a leading role and the degree of nonlinearity increases. It is stronger than the bottom of the collapse column, and F0 shows an increasing trend, indicating that with the loss of mass, the nonlinearity of the flow state of the water in the middle increases. At time is 1.6×104 s, the position F0 at the bottom of the collapse column is around 1 and drops to 1 at t=20000 s. In contrast, the middle position drops and rises abruptly at t=1.6×104 s and t=2.0×104 s, respectively, which indicates that the mass loss in the middle blocks the channel and leads to nonlinear weakening. The nonlinearity shows the opposite law with the bottom and middle of the collapse column.Figure 17
Forchheimer numbers of domain points at different times.


## 5. Conclusions

The seepage tests are carried out on the broken rock mass with different gradations under triaxial stress, and the nonlinear seepage coupling equations of the three-flow field are established in this paper. Besides, finite element software COMSOL Multiphysics was employed to simulate the seepage process of the three-flow field.(1)
The reasonable degree of fluid flow fitting curve of fractured rock mass and Forchheimer equation under triaxial stress with differentn value gradations is above 0.96, and Darcy’s law is reasonable 0.65-0.75(2)
From 0 s to 20000 s, the pressure from the aquifer to the collapse column boundary decreases fromp=2.02 MPa to p=1.74 MPa, and the velocity increases from 0.015 mm/s to 0.064 mm/s. The porosity of domain point 3 changed from 0.113 to 0.331, and the permeability increased from 1.03e−12 to 6.4e−10. The water volume increased from 43.2 m3/h to 176.8 m3/h(3)
Before the 12000 s, the tiny fissures developed and expanded first. The minimum fissure increased from9.98×10−6 m to 1.25e−5 m. When the time is 16000 s, the maximum crack width starts at 4.06×10−4 m, which is 1.58 times that of t=12000 s. Within 16000 s to 20000 s, the maximum width of the crack is 6.83×10−4 m. The velocity cusp appeared at the boundary position, reaching 3.96 mm/s, and the water inflow nearly doubled, reaching 176.8 m3/h(4)
Particles enter the fluid under the action of water seepage-erosion, which causes the increase of the porosity and permeability of the collapse column, the increase of the hydrodynamic viscosity, the change of the flow state, and the Forchheimer number increases from 0.72 to 1 when it is close to the aquifer. The Forchheimer number in the middle of the collapse column fluctuates between 7.5, and the degree of nonlinearity is small from laminar flow to nonlinear seepage



---

*Source: 8148203-2022-07-25.xml*