- Email: [email protected]
Bifurcation and catastrophe of seepage flow system in broken rock
M INING SCIENCE AND TECHNOLOGY Mining Science and Technology 19 (2009) 0001–0007
www.elsevier.com/locate/jcumt
Bifurcation and catastrophe of seepage flow system in broken rock MIAO Xie-xing1, LI Shun-cai2, CHEN Zhan-qing1 1
2
School of Sciences, China University of Mining & Technology, Xuzhou, Jiangsu 221008, China School of Mechanical and Electrical Engineering, Xuzhou Normal University, Xuzhou, Jiangsu 221116, China
Abstract: The study of dynamical behavior of water or gas flows in broken rock is a basic research topic among a series of key projects about stability control of the surrounding rocks in mines and the prevention of some disasters such as water inrush or gas outburst and the protection of the groundwater resource. It is of great theoretical and engineering importance in respect of promotion of security in mine production and sustainable development of the coal industry. According to the non-Darcy property of seepage flow in broken rock dynamic equations of non-Darcy and non-steady flows in broken rock are established. By dimensionless transformation, the solution diagram of steady-states satisfying the given boundary conditions is obtained. By numerical analysis of low relaxation iteration, the dynamic responses corresponding to the different flow parameters have been obtained. The stability analysis of the steady-states indicate that a saddle-node bifurcaton exists in the seepage flow system of broken rock. Consequently, using catastrophe theory, the fold catastrophe model of seepage flow instability has been obtained. As a result, the bifurcation curves of the seepage flow systems with different control parameters are presented and the standard potential function is also given with respect to the generalized state variable for the fold catastrophe of a dynamic system of seepage flow in broken rock. Keywords: broken rock; non-Darcy flow; stability; saddle-node bifurcation; fold catastrophe
1
Introduction
When mining, rocks (coal) undergo fractures or break and, consequently, their permeability increases drastically. At the same time, groundwater greatly affects the deformation and stability of the rock structure by hydraulic fractures and all kinds of other physical and chemical actions such as lubrication, erosion, hydroxylation and coupling. Hence, water inrush accidents are apt to take place in the broken zones of aquifers[1–3]. For example, a disaster of water inrush took place in the Fangezhuang mine, with the greatest water discharge of 2053 cubic meters per minute on June 2, 1984. The origin of this disaster was that the water filled the collapsed column connected to high-pressure karst water of the Ordovician. Some areas with seriously weathered wall rock in the railway tunnel easily developed a broken belt and caused the inrush of water and mud[4]. Much experimental research has been carried out by earlier scientists on the seepage flow in broken rock (or rock fill). Some basic formulas of non-Darcy flows have been obtained by Forchheimer, Irmay and
others, but these formulas all assume a solid phase of rigid porous media[5–7]. However, for broken rock acting as the bearing structure in underground engineering projects, the edges and corners of the grains are prone to break under the load, so its permeability changes markedly with its load or displacement. From an experiment in steady seepage on broken sandstone, tailings and limestone of different grain sizes, 5–10 mm, 10–15 mm, 15–20 mm, 20–25 mm and a compounded grain size, Li Shun-cai has drawn the following conclusions[8–9]: 1) The water seepage flow in the broken rock is a non-Darcy flow. When the flow is steady, the relation ∂p between the pore pressure gradient and the flow ∂x velocity v can be fitted with the following second ∂p μ degree polynomial: = − v − bv 2 , where the ∂x k quadratic term bv 2 is a modification of the linear μ term v . The coefficient b is called the coefficient k [8] of deviation from the Darcy-flow .
Received 13 May 2008; accepted 14 September 2008 Projects 50490273 and 50674087 supported by the National Natural Science Foundation of China and BK2007029 by the Natural Science Foundation of Jiangsu Province Corresponding author. Tel: +86-516-83899037; E-mail address: [email protected]
2
Mining Science and Technology
2) Over its entire range, the permeability k decreases, while the absolute value of coefficient b increases with a decrease in porosity φ , but both the k- φ and the b- φ curves incur some local fluctuation. The greater the grain size, the more fluctuation of the curves. 3) The permeability k of the five kinds of broken rock of different grain sizes has a magnitude of 10−16 – 10−12 m2, while the coefficient b ranges from 4 1012 – 1015 kg/m ; 4) The coefficient b exists possibly in both positive and negative forms. When the grains become seriously broken, the coefficient b is most likely to appear as a negative value. At this time, the seepage flow in the broken rock may show the condition that the superficial permeability is greater than its absolute permeability[10]. In relatively intact rock, the flow is extremely slow and can be studied as a linear and steady course. But in broken or fractured rock, the flow is nonlinear, unsteady and variable over time[11–15]. Especially, in mining engineering, the gradual variation of the seepage parameters of the flow in surrounding rock affected by mining is apt to lead to flow instability and catastrophe, thus resulting in some dynamic disasters such as water inrush or gas outbursts. Hence, dynamical investigation of nonlinear flows in broken rock is of great theoretical importance with future application.
2
The complete, dynamic equations for non-Darcy and non-steady flow in broken rock are composed of three parts[8], i.e., the continuity Eq.(1) the kinetic flow Eq.(2) and the state Eqs.(3) and (4).
(1)
∂v ∂p μ ρca = − − v − bv 2 − ρg ∂t ∂x k
(2)
ρ = ρ 0 [1 + cf ( p − p0 )]
(3)
φ = φ0 [1 + cφ ( p − p0 )]
(4)
ρ 0φ0 ct
∂p ∂ ( ρv) =− ∂t ∂x
(5)
where ct = cφ + cf is called the composite compression coefficient. Then from Eqs.(2) and (5) we can obtain the dynamic equations of the one-dimensional non-Darcy flow in broken rock as follows: ∂p ∂ ( ρv) °° ρ 0φ0 ct ∂t = − ∂x ® ° ρc ∂v = − ∂p − μ v − bv 2 − ρg °¯ a ∂t ∂x k
∂v ∂p °° ∂t = − a0 ∂x ® ° ∂v = − a ∂p − a v − a v 2 − a 1 2 3 4 °¯ ∂t ∂x
(8)
2
p §bk · ¸ , where a0 = , a1 = 0 ¨¨ p0φ0 ct ca ρ 0 © μ ¸¹ 2
§ b k · Hg bH , a4 = ¨ a2 = a3 = ¸ ca ρ0 © μ ¹ ca
(9)
3 Steady-states of the flow system Fig. 1 is a sketch map of the flow in broken rock. The boundary conditions of the flow are as follows: p x = 0 = p1 , p x = H = p2 and the flow is along the positive direction of axis x, that is to say, p1 − p2 > ρ0 gH , where H is the piling height of the broken rock.
mass density of the fluid, ca the acceleration coefficient[16], b the coefficient of deviation from the Darcy-flow, p the relative pore pressure with regard to the atmosphere pressure, μ the dynamic viscosity
φ 0 , ρ 0 are the porosity and mass density corre-
(6)
By ignoring the compressibility of the water, we can carry out a dimensionless transformation as follows: x v t p , x= , v= , t= (7) p= μ (bk ) b kH μ p0 H Then the Eq.(6) can be simplified as
where v is the flow velocity, φ the porosity, ρ the
of the fluid, k the permeability of the broken rock, cf the isothermal compression coefficient of the fluid, cφ the compression coefficient of the pore and
No.1
sponding to the reference pressure p 0 (relative to the atmosphere pressure) respectively. When multiplying Eq.(3) with Eq.(4) and neglecting the terms containing two order compression coefficients and taking this product into the Eq.(1), we obtain
1
Dynamic model
∂ ( ρφ ) ∂ ( ρv) + =0 ∂t ∂x
Vol.19
Fig. 1
Sketch map of flow in broken rock
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
According to Eq.(8), the steady-states ( ps , vs ) of the flow system satisfy ∂vs =0 °° ∂x ® °a ∂ps + a v + a v 2 + a = 0 2 s 2 s 4 °¯ 1 ∂x
presented as Fig. 2. It follows that vs is composed of three parts, i.e., solution 1 ( b > 0, vs > 0 ), solution 2 ( b < 0, vs > −0.5 ) and solution 3 ( b < 0, vs < −0.5 ).
(10)
vs 6ROXWLRQ
Making use of the boundary conditions of the pore pressure p1 x = 0 = p1 p0 , p2 x =1 = p2 p0 , we can
6ROXWLRQ
p1 ( p1 − p2 ) x − p0 p0
− a2 ± a2
2
=
4a2 + ρ 0 ca
2
According to these two equations, the pore pressure along the flow direction decreases linearly but the flow velocity of each point is equal when the flow is steady. Because the flow velocity v is along the positive direction of the axis x, i.e., v > 0 , we have, from Eq.(7), the condition that v > 0 corresponding to b > 0 and v < 0 corresponding to b < 0 . The steady-states of flow velocity corresponding to Eq.(12) can be developed as follows (a detailed development can be found in Reference [8]) 1 1 4b k 2 ( p1 − p2 − ρ 0 gH ) 1+ , °− + μ2H ° 2 2 ° 2 °− 1 # 1 1 + 4bk ( p1 − p2 − ρ0 gH ) , °° 2 2 μ2H vs = ® ° μ2H
b>0
(13)
2
μ H 2
4k ( p1 − p2 − ρ 0 gH )
Solutions map of seepage velocity for water flow system of broken rock
Stability analyses of the steady-states of flow system
Assume the initial conditions of the flow system p −p are as follows: pore pressure p0 ( x) = p01 + 02 01 x H and flow velocity v0 ( x) = v0 , where p01 , p02 are respectively the initial pore pressures corresponding to the top and the bottom boundary of the pile of broken rock and the flow velocity is along the positive direction of the x axis. By dimensionless transformation, we obtain p0 ( x ) =
bk p01 p02 − p01 v0 , x ∈ [0,1] + x , v0 ( x ) = μ p0 p0
The responses corresponding to Eq.(8) can be obtained by iteration of successive low relaxation. In Eq.(8), partitioning the nodes along the direction of the time and space and designating the subscript “i” and “j” for the time and space variables respectively, then using the forward difference formula for the partial derivatives with respect to time and using the central difference formula for the partial derivatives with respect to height x (except the boundary nodes), we obtain vi , j +1 − vi, j −1 ° pi+1, j = pi, j + Δt ⋅ (−a0 ) 2Δx ° ® p ª º i , j +1 − pi , j −1 °v − a2vi , j − a3vi, j 2 − a4 » i +1, j = vi , j + Δt ⋅ «− a1 ° 2Δx ¬ ¼ ¯
Define bs = −
Fig. 2
4
§ bk · ¨¨ ¸¸ ( p1 − p2 − ρ 0 gH ) ©μ ¹ (12) 2a2
β
6ROXWLRQ
(11)
§ p − p2 · 2 − a2 ± a2 + 4a2 ¨¨ a1 1 − a4 ¸¸ p0 © ¹ vs = 2a2
2
$ EV
solve these differential equations, after which the pore pressure ps and the flow velocity vs corresponding to the steady-states can be developed as ps ( x ) =
3
(14)
where bs is called the critical coefficient of deviation from the Darcy-flow. The solution diagram corresponding to the dimensionless flow velocity vs is
We used the following iteration formula of low relaxation to carry out the iteration for every node, pi +1, j ← (1 − ω ) pi , j + ω pi +1, j , qi +1, j ← (1 − ω ) qi , j + ω qi +1, j where w is a relaxation factor. We then obtained the
4
Mining Science and Technology
time series of the pore pressure and the flow velocity corresponding to each node by iteration. Assume p1 = 0.6 MPa and p2 = 0.2 MPa for the boundary with fixed pore pressure and the initial porosity φ0 = 0.32 corresponding to the reference pressure p0 = 0.3 MPa , the dynamic viscosity μ = 1.01×10−3 (Pa ⋅ s) , permeability k = 0.5 × 10 −13 (m 2 ) ,
the acceleration coefficient ca = 9.46 × 109 , the compression coefficient of the fluid cf = 4.75 × 10 −10 Pa–1, the compression coefficient of the pore
cφ = 2.02 ×10 −9 Pa–1, ρ0 = 1000 (kg m3 ) , height H = 5 m , Δx = 0.1 (here the piling height is equally divided into ten parts), time interval Δt =1 s, ac,
hence Δt = 1.13 × 10−3 and the relaxation factor ω = 0.6 . Some cases will be discussed as follows.
(a) Curve of dimensionless pore pressure vs. time
Vol.19
No.1
4.1 Analysis to the stability of steady-state when b>0 Arbitrarily we selected the parameter b satisfying b > 0 , for example b = 3.5847 × 1012 kg/m4 and then, according to Eqs.(11) and (12), we can develop the dimensionless steady states of pore pressure and flow velocity respectively as follows: 4 ps ( x ) = 2.0 − x , vs = 6.145 × 10 −4 3 When the initial pore pressure and flow velocity show a little deviation from their steady-states and if let p01 = 0.50 MPa , p02 = 0.0 MPa and v0 = 6.323×10−4 , then after 20000 iterations, the time series of the pore pressure and the flow velocity of the five nodes corresponding to x=0.1H, 0.3H, 0.5H, 0.7H and 0.9H are shown in Figs. 3a and b.
(b) Curve of dimensionless flow velocity vs. time
Fig. 3 Stability analysis of steady-state corresponding to b > 0
It is clear that the pore pressure and flow velocity of each node are all attracted to its steady-state in a stable fashion, so the steady-states corresponding to b > 0 are stable.
4.2
Analysis to stability of steady-state when bs < b < 0
As shown in Fig. 2, when bs < b < 0 the steady-states have two parts, i.e., solutions 2 and 3 with the critical point A (bsˈ− 0.5) as their boundary, where for each parameter b, there are two steady-states of flow velocity, respectively noted as v s 2 and v s 3 . Making use of Eq.(14) and the parameters k , ca , H, p1 , p 2 , μ and ρ 0 , mentioned before, the corresponding critical parameter b can be obtained, i.e., bs = −1.457 × 1015 (kg/m4). Thus, when b changes and satisfies bs < b < 0 , by numerical analysis, we can observe the stability conditions of the steady-states in every branch. If we arbitrarily select the parameter b satisfying bs < b < 0 , for example, b = −1.457 × 1014 (kg/m4), then we obtain the two steady-states of flow velocity corresponding to the solution diagram (Fig. 2), i.e.,
vs2 = −0.02563 and vs3 = −0.97431 . 4.2.1 Stability of steady-states corresponding to the solution 2 When we study the stability of the steady-state vs2 (0.02563), we can discuss two cases in respect of the initial flow velocity v02 , which can either be larger
or smaller than the steady-state velocity vs2 : 1) v02 > vs2 , let v02 = −0.02493 , i.e., the initial velocity deviates from the steady-state and lies on its high side. By evolving, the flow velocity eventually comes to its steady-state vs2 (–0.02563) in a stable way, as shown in Fig. 4, where only the curves of flow state variables vs time at x=0.5H are presented. In particular, it should be pointed out that when b < 0 , the time variables in all the curves below are the magnitude or the absolute value of dimensionless time t . 2) v02 < vs2 , if we let v02 = −0.02573 , by evolving, the flow velocity eventually comes to its steady-state vs2 (–0.02563) in a stable way, as shown in Fig. 5. It is obvious that the steady-states corresponding to solution 2 are stable nodes.
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
(a) Curve of dimensionless pore pressure vs. time
Fig. 4
(b) Dimensionless flow velocity vs. time
Dynamic responses corresponding to v02 < vs2 for branch 2
4.2.2 Stability of steady-states corresponding to solution 3 We continue to analyze the stability of steadystates corresponding to solution 3 and also discuss the stability according to two cases: 1) v03 > vs3 , if we let v03 = −0.9671 , then by evolv-
(a) Dimensionless pore pressure vs. time
Fig. 6
(b) Curve of dimensionless pore pressure vs. time
Dynamic responses corresponding to v02 > vs2 for branch 2
(a) Dimensionless pore pressure vs. time
Fig. 5
5
ing, the flow velocity in the end does not come to its corresponding steady-state vs3 (–0.9743); on the contrary, it is attracted to the steady-state vs2 (–0.02563) corresponding to solution 2, as shown in Fig. 6.
(b) Dimensionless flow velocity vs. time
Dynamic responses corresponding to v03 > vs3 for branch 3
2) v03 < vs3 , if we let v03 = −1.0464 and the time space is the same as above, then by evolving, the flow velocity in the end does not come either to its corresponding steady-state vs3 (–0.9743); on the contrary, it deviates from the steady-state farther and farther and moreover, is attracted to negative infinity as shown in Fig. 7. It is obvious that the steady-states corresponding to solution 3 are unstable. When the initial velocity (dimensionless) is larger than its steady-state, its phase trajectory is, in the end, attracted to the steady-states corresponding to branch 2. When the initial flow velocity is less than its steady-state, its phase trajectory deviates from its steady-state and is attracted to nega-
tive infinity. So, the steady-states corresponding to solution 3 are saddle points. 4.2.3 Analysis to the flow stability when b < bs When b < bs , there are no steady-states. If we select the parameter b satisfying b < bs , for example, b = −1.457 × 1016 < bs = −1.457 × 1015 (kg/m4), then corresponding to a freely given initial velocity v0 (no matter how small), by evolving, the dimensionless flow velocity in an absolute sense tends to negative infinity, i.e., the flow system lost its stability. The dynamic responses corresponding to the given b and v0 = 2.63 × 10−6 (m/s) are shown in Fig. 8.
6
Mining Science and Technology
Vol.19
(a) Dimensionless flow velocity vs. time
Fig. 7
(b) Phase trajectory
Dynamic responses corresponding to v03 < vs3 for branch 3
(a) Dimensionless flow velocity vs. time
Fig. 8
5
No.1
(b) Phase trajectory
Unstable motion of non-darcy flow system when b < bs
Studies on bifurcation and catastrophe of flow in broken rock
From our investigation, we know that when b < 0 , there are two branches of the steady-states. The steady-states corresponding to solution 2 are stable nodes while those corresponding to solution 3 are unstable saddles. When b decreases gradually and tends to the critical parameter bs , the saddle and the
larger). Yet, for the different flow system they all satisfy vs = −0.5 when b = bs . So, we can obtain different bifurcation curves corresponding to the flow systems with different controlling parameter as shown in Fig. 10, where the real line and the dotted line are, respectively, corresponding to stable and unstable flows.
nodes meet and annihilate each other. When b < bs the steady-states vanish, a condition where all the phase-trajectories flow to negative infinity (only one trajectory is drawn in Fig. 8). When bs < b < 0 , the phase trajectories are divided into two flow areas with 4b k 2 ( p1 − p2 − ρ0 gH ) 1 1 branch 3, i.e., vs3 = − − 1 + 2 2 μ2H serving as their boundary, that is to say, for the trajectories corresponding to v0 > vs3 , their flow velocities finally tend to the steady-state vs2 = −
Fig. 9
Saddle-node bifurcation of flow in broken rock
4b k 2 ( p1 − p2 − ρ 0 gH ) 1 1 , while for + 1+ 2 2 μ 2H
the trajectories corresponding to v0 < vs3 ; their flow velocities finally tend to negative infinity. Hence, there is a saddle-node bifurcation (fold bifurcation) at b = bs in the dynamic system of flows in broken rock. The corresponding sketch map of bifurcation is shown in Fig. 9. For different flow systems, the permeability k and critical parameter bs also vary. Moreover, the smaller k, the smaller bs (but its absolute value is
Fig. 10
Bifurcation curve for the flow systems with different controlling parameter
Saddle-node bifurcation (fold bifurcation) is a typical discontinuous bifurcation. Catastrophes of stability of the flow system take place at the bifurcation point. According to the theory of catastrophe[17], a fold catastrophe is apt to occur at the point of fold bifurcation where, at this time, any tiny disturbance can induce some dynamic disasters such as water in-
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
rush or gas bursts. Given the second equation in Eq.(8), the potential function U , with regard to the generalized variable v , should satisfy ∂U ∂p = a2v 2 + a2v + a1 + a4 (15) ∂v ∂x ∂p p2 − p1 where the pore pressure gradient = =C . ∂x p0 Integration of Eq.(15) obtains: 1 1 U = a2v 3 + a2v 2 + a4v + a1Cv (16) 3 2 Developing the first and second order derivatives with respect to v , from Eq.(16), the set S of singularity points satisfies ∂U = a2v 2 + a2v + a1C + a4 = 0 ∂v ∂ 2U = 2 a 2 v + a2 = 0 ∂v 2 Solving these two equations jointly, we obtain the singularity point
μ 2H (17) 4k 2 ( p1 − p2 − ρ0 gH ) which is precisely the bifurcation point obtained earlier. In order to acquire the standard potential function, by coordinate translation, let v = v′ + d (18) Substituting this formula into Eq.(16) and letting d = −1 2 to eliminate the term containing the square of v ′ , we have v = −0.5 , b = −
a 1 U = a2 v ′3 + ( − 2 + a4 + a1C )v ′ 3 4 1 a4 + (− a1C − + c0 ) (19) 2 2 where c0 is the constant of integration. Eq.(19) is the standard potential function with respect to the generalized state variable v ′ for the fold catastrophe of the flow system in broken rock.
6
Conclusions
Given the permeating experimental results of a non-Darcy flow in broken rock, the nonlinear dynamic equations of non-Darcy and non-steady seepage flows in broken rock have been established. By dimensionless transformation, the solution diagram of the steady-states under the given boundary conditions is obtained for the dynamic system of the seepage flow. With a change in the seepage flow parameter b, the steady-states of flow velocity are divided into three parts. Numerical analyses indicate that there is a saddle-node bifurcation (fold bifurcation) at the critical parameter bs for the seepage flow in broken rock. Bifurcation curves of the different seepage flow sys-
7
tems have also been provided. Near the bifurcation point, catastrophe of the stability of the steady-states takes place when, at the same time, the breaking phenomenon of the solid grains of the broken rock is rather serious. So for the seepage flow system in the broken rock catastrophe is apt to occur at the bifurcation point and will result in dynamic disasters such as a water inrush.
References [1]
[2]
[3] [4]
[5] [6] [7] [8] [9]
[10] [11]
[12] [13]
[14] [15]
[16] [17]
Zhu X S. The application of engineering geophysics in the remediation of water inrush in fractured zone. Express Information of Mining Industry, 2004(11): 51–53. (In Chinese) Wang L, Liu C A. Analysis and disposal to the flooded mine resulting from the water inrush in breaking zone of shaft. Mine Construction Technology, 1999, 20(2): 8–11. (In Chinese) Zhou R G, Cheng B F, Ye G J, Wu Q. The effect of water busting in fault rupture zone. Journal of Engineering Geology, 2000, 8(4): 411–415. (In Chinese) Li F G. Analysis and disposal of seepage and water gushing in tunnel fault broken zone. West-China Exploration Engineering, 2004, 132(12): 112–114. (In Chinese) Stephson D. Hydraulic Calculation in Rockfill Engineering. Beijing: Ocean Publishing Company, 1984. Wang Z F. Notes of rockfill seepage research. Journal of Nanjing Hydraulic Research Institute, 1990(4): 445–452. (In Chinese) Bear J. Dynamics of Fluid in Porous Media. New York: Elsevier, 1979. Li S C. Nonlinear Dynamical Study on Non-darcy Flow in Broken Rock [Ph.D. dissertation]. Xuzhou: China University of Mining & Technology, 2006. (In Chinese) Miao X X, Li S C, Huang X W. Experimental study of seepage properties of non-Darcy flow in granular gangues. Journal of University of Mining and Technology, 2006, 16(2): 105–109. (In Chinese) Deng Y E, Liu C Q, Huang R Q. Advanced Seepage Theory and Method. Beijing: Scientific Publishing Company, 2004. (In Chinese) Chen Z Q. Bifurcation of Dynamic System of Non-Darcy Flow in Post-Failure Rock [Ph.D. dissertation]. Xuzhou: China University of Mining & Technology, 2003. (In Chinese) Li T Z, Li Y S, Ma Z G. Testing study of non-Darcy seepage flow in fractured rocks. Engineering Mechanics, 2003, 20(4): 132–135. (In Chinese) Tang P, Huang X W, Li T Z. Test study on the penetrating properties of non-Darcy flow in rock. Journal of Xiangtan Normal University, 2002, 24(2): 34–38. (In Chinese) Miao X X, Chen Z Q, Mao X B, Chen R H. The bifurcation of non-Darcy flow in post-failure rock. Acta Mechanica Sinica, 2003, 35(6): 660–667. (In Chinese) Pradip K G N, Venkataraman P. Non-Darcy converging flow through coarse granular media. Journal of the Institution of Engineers (India), Civil Engineering Division, 1995(76): 6–11. Kong X Y. Advanced Mechanics of Fluid in Porous Media. Hefei: Press of China University of Science and Technology, 1999. (In Chinese) Liu B Z, Peng J H. Nonlinear Dynamics. Beijing: Higher Education Press, 2003. (In Chinese)
www.elsevier.com/locate/jcumt
Bifurcation and catastrophe of seepage flow system in broken rock MIAO Xie-xing1, LI Shun-cai2, CHEN Zhan-qing1 1
2
School of Sciences, China University of Mining & Technology, Xuzhou, Jiangsu 221008, China School of Mechanical and Electrical Engineering, Xuzhou Normal University, Xuzhou, Jiangsu 221116, China
Abstract: The study of dynamical behavior of water or gas flows in broken rock is a basic research topic among a series of key projects about stability control of the surrounding rocks in mines and the prevention of some disasters such as water inrush or gas outburst and the protection of the groundwater resource. It is of great theoretical and engineering importance in respect of promotion of security in mine production and sustainable development of the coal industry. According to the non-Darcy property of seepage flow in broken rock dynamic equations of non-Darcy and non-steady flows in broken rock are established. By dimensionless transformation, the solution diagram of steady-states satisfying the given boundary conditions is obtained. By numerical analysis of low relaxation iteration, the dynamic responses corresponding to the different flow parameters have been obtained. The stability analysis of the steady-states indicate that a saddle-node bifurcaton exists in the seepage flow system of broken rock. Consequently, using catastrophe theory, the fold catastrophe model of seepage flow instability has been obtained. As a result, the bifurcation curves of the seepage flow systems with different control parameters are presented and the standard potential function is also given with respect to the generalized state variable for the fold catastrophe of a dynamic system of seepage flow in broken rock. Keywords: broken rock; non-Darcy flow; stability; saddle-node bifurcation; fold catastrophe
1
Introduction
When mining, rocks (coal) undergo fractures or break and, consequently, their permeability increases drastically. At the same time, groundwater greatly affects the deformation and stability of the rock structure by hydraulic fractures and all kinds of other physical and chemical actions such as lubrication, erosion, hydroxylation and coupling. Hence, water inrush accidents are apt to take place in the broken zones of aquifers[1–3]. For example, a disaster of water inrush took place in the Fangezhuang mine, with the greatest water discharge of 2053 cubic meters per minute on June 2, 1984. The origin of this disaster was that the water filled the collapsed column connected to high-pressure karst water of the Ordovician. Some areas with seriously weathered wall rock in the railway tunnel easily developed a broken belt and caused the inrush of water and mud[4]. Much experimental research has been carried out by earlier scientists on the seepage flow in broken rock (or rock fill). Some basic formulas of non-Darcy flows have been obtained by Forchheimer, Irmay and
others, but these formulas all assume a solid phase of rigid porous media[5–7]. However, for broken rock acting as the bearing structure in underground engineering projects, the edges and corners of the grains are prone to break under the load, so its permeability changes markedly with its load or displacement. From an experiment in steady seepage on broken sandstone, tailings and limestone of different grain sizes, 5–10 mm, 10–15 mm, 15–20 mm, 20–25 mm and a compounded grain size, Li Shun-cai has drawn the following conclusions[8–9]: 1) The water seepage flow in the broken rock is a non-Darcy flow. When the flow is steady, the relation ∂p between the pore pressure gradient and the flow ∂x velocity v can be fitted with the following second ∂p μ degree polynomial: = − v − bv 2 , where the ∂x k quadratic term bv 2 is a modification of the linear μ term v . The coefficient b is called the coefficient k [8] of deviation from the Darcy-flow .
Received 13 May 2008; accepted 14 September 2008 Projects 50490273 and 50674087 supported by the National Natural Science Foundation of China and BK2007029 by the Natural Science Foundation of Jiangsu Province Corresponding author. Tel: +86-516-83899037; E-mail address: [email protected]
2
Mining Science and Technology
2) Over its entire range, the permeability k decreases, while the absolute value of coefficient b increases with a decrease in porosity φ , but both the k- φ and the b- φ curves incur some local fluctuation. The greater the grain size, the more fluctuation of the curves. 3) The permeability k of the five kinds of broken rock of different grain sizes has a magnitude of 10−16 – 10−12 m2, while the coefficient b ranges from 4 1012 – 1015 kg/m ; 4) The coefficient b exists possibly in both positive and negative forms. When the grains become seriously broken, the coefficient b is most likely to appear as a negative value. At this time, the seepage flow in the broken rock may show the condition that the superficial permeability is greater than its absolute permeability[10]. In relatively intact rock, the flow is extremely slow and can be studied as a linear and steady course. But in broken or fractured rock, the flow is nonlinear, unsteady and variable over time[11–15]. Especially, in mining engineering, the gradual variation of the seepage parameters of the flow in surrounding rock affected by mining is apt to lead to flow instability and catastrophe, thus resulting in some dynamic disasters such as water inrush or gas outbursts. Hence, dynamical investigation of nonlinear flows in broken rock is of great theoretical importance with future application.
2
The complete, dynamic equations for non-Darcy and non-steady flow in broken rock are composed of three parts[8], i.e., the continuity Eq.(1) the kinetic flow Eq.(2) and the state Eqs.(3) and (4).
(1)
∂v ∂p μ ρca = − − v − bv 2 − ρg ∂t ∂x k
(2)
ρ = ρ 0 [1 + cf ( p − p0 )]
(3)
φ = φ0 [1 + cφ ( p − p0 )]
(4)
ρ 0φ0 ct
∂p ∂ ( ρv) =− ∂t ∂x
(5)
where ct = cφ + cf is called the composite compression coefficient. Then from Eqs.(2) and (5) we can obtain the dynamic equations of the one-dimensional non-Darcy flow in broken rock as follows: ∂p ∂ ( ρv) °° ρ 0φ0 ct ∂t = − ∂x ® ° ρc ∂v = − ∂p − μ v − bv 2 − ρg °¯ a ∂t ∂x k
∂v ∂p °° ∂t = − a0 ∂x ® ° ∂v = − a ∂p − a v − a v 2 − a 1 2 3 4 °¯ ∂t ∂x
(8)
2
p §bk · ¸ , where a0 = , a1 = 0 ¨¨ p0φ0 ct ca ρ 0 © μ ¸¹ 2
§ b k · Hg bH , a4 = ¨ a2 = a3 = ¸ ca ρ0 © μ ¹ ca
(9)
3 Steady-states of the flow system Fig. 1 is a sketch map of the flow in broken rock. The boundary conditions of the flow are as follows: p x = 0 = p1 , p x = H = p2 and the flow is along the positive direction of axis x, that is to say, p1 − p2 > ρ0 gH , where H is the piling height of the broken rock.
mass density of the fluid, ca the acceleration coefficient[16], b the coefficient of deviation from the Darcy-flow, p the relative pore pressure with regard to the atmosphere pressure, μ the dynamic viscosity
φ 0 , ρ 0 are the porosity and mass density corre-
(6)
By ignoring the compressibility of the water, we can carry out a dimensionless transformation as follows: x v t p , x= , v= , t= (7) p= μ (bk ) b kH μ p0 H Then the Eq.(6) can be simplified as
where v is the flow velocity, φ the porosity, ρ the
of the fluid, k the permeability of the broken rock, cf the isothermal compression coefficient of the fluid, cφ the compression coefficient of the pore and
No.1
sponding to the reference pressure p 0 (relative to the atmosphere pressure) respectively. When multiplying Eq.(3) with Eq.(4) and neglecting the terms containing two order compression coefficients and taking this product into the Eq.(1), we obtain
1
Dynamic model
∂ ( ρφ ) ∂ ( ρv) + =0 ∂t ∂x
Vol.19
Fig. 1
Sketch map of flow in broken rock
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
According to Eq.(8), the steady-states ( ps , vs ) of the flow system satisfy ∂vs =0 °° ∂x ® °a ∂ps + a v + a v 2 + a = 0 2 s 2 s 4 °¯ 1 ∂x
presented as Fig. 2. It follows that vs is composed of three parts, i.e., solution 1 ( b > 0, vs > 0 ), solution 2 ( b < 0, vs > −0.5 ) and solution 3 ( b < 0, vs < −0.5 ).
(10)
vs 6ROXWLRQ
Making use of the boundary conditions of the pore pressure p1 x = 0 = p1 p0 , p2 x =1 = p2 p0 , we can
6ROXWLRQ
p1 ( p1 − p2 ) x − p0 p0
− a2 ± a2
2
=
4a2 + ρ 0 ca
2
According to these two equations, the pore pressure along the flow direction decreases linearly but the flow velocity of each point is equal when the flow is steady. Because the flow velocity v is along the positive direction of the axis x, i.e., v > 0 , we have, from Eq.(7), the condition that v > 0 corresponding to b > 0 and v < 0 corresponding to b < 0 . The steady-states of flow velocity corresponding to Eq.(12) can be developed as follows (a detailed development can be found in Reference [8]) 1 1 4b k 2 ( p1 − p2 − ρ 0 gH ) 1+ , °− + μ2H ° 2 2 ° 2 °− 1 # 1 1 + 4bk ( p1 − p2 − ρ0 gH ) , °° 2 2 μ2H vs = ® ° μ2H
b>0
(13)
2
μ H 2
4k ( p1 − p2 − ρ 0 gH )
Solutions map of seepage velocity for water flow system of broken rock
Stability analyses of the steady-states of flow system
Assume the initial conditions of the flow system p −p are as follows: pore pressure p0 ( x) = p01 + 02 01 x H and flow velocity v0 ( x) = v0 , where p01 , p02 are respectively the initial pore pressures corresponding to the top and the bottom boundary of the pile of broken rock and the flow velocity is along the positive direction of the x axis. By dimensionless transformation, we obtain p0 ( x ) =
bk p01 p02 − p01 v0 , x ∈ [0,1] + x , v0 ( x ) = μ p0 p0
The responses corresponding to Eq.(8) can be obtained by iteration of successive low relaxation. In Eq.(8), partitioning the nodes along the direction of the time and space and designating the subscript “i” and “j” for the time and space variables respectively, then using the forward difference formula for the partial derivatives with respect to time and using the central difference formula for the partial derivatives with respect to height x (except the boundary nodes), we obtain vi , j +1 − vi, j −1 ° pi+1, j = pi, j + Δt ⋅ (−a0 ) 2Δx ° ® p ª º i , j +1 − pi , j −1 °v − a2vi , j − a3vi, j 2 − a4 » i +1, j = vi , j + Δt ⋅ «− a1 ° 2Δx ¬ ¼ ¯
Define bs = −
Fig. 2
4
§ bk · ¨¨ ¸¸ ( p1 − p2 − ρ 0 gH ) ©μ ¹ (12) 2a2
β
6ROXWLRQ
(11)
§ p − p2 · 2 − a2 ± a2 + 4a2 ¨¨ a1 1 − a4 ¸¸ p0 © ¹ vs = 2a2
2
$ EV
solve these differential equations, after which the pore pressure ps and the flow velocity vs corresponding to the steady-states can be developed as ps ( x ) =
3
(14)
where bs is called the critical coefficient of deviation from the Darcy-flow. The solution diagram corresponding to the dimensionless flow velocity vs is
We used the following iteration formula of low relaxation to carry out the iteration for every node, pi +1, j ← (1 − ω ) pi , j + ω pi +1, j , qi +1, j ← (1 − ω ) qi , j + ω qi +1, j where w is a relaxation factor. We then obtained the
4
Mining Science and Technology
time series of the pore pressure and the flow velocity corresponding to each node by iteration. Assume p1 = 0.6 MPa and p2 = 0.2 MPa for the boundary with fixed pore pressure and the initial porosity φ0 = 0.32 corresponding to the reference pressure p0 = 0.3 MPa , the dynamic viscosity μ = 1.01×10−3 (Pa ⋅ s) , permeability k = 0.5 × 10 −13 (m 2 ) ,
the acceleration coefficient ca = 9.46 × 109 , the compression coefficient of the fluid cf = 4.75 × 10 −10 Pa–1, the compression coefficient of the pore
cφ = 2.02 ×10 −9 Pa–1, ρ0 = 1000 (kg m3 ) , height H = 5 m , Δx = 0.1 (here the piling height is equally divided into ten parts), time interval Δt =1 s, ac,
hence Δt = 1.13 × 10−3 and the relaxation factor ω = 0.6 . Some cases will be discussed as follows.
(a) Curve of dimensionless pore pressure vs. time
Vol.19
No.1
4.1 Analysis to the stability of steady-state when b>0 Arbitrarily we selected the parameter b satisfying b > 0 , for example b = 3.5847 × 1012 kg/m4 and then, according to Eqs.(11) and (12), we can develop the dimensionless steady states of pore pressure and flow velocity respectively as follows: 4 ps ( x ) = 2.0 − x , vs = 6.145 × 10 −4 3 When the initial pore pressure and flow velocity show a little deviation from their steady-states and if let p01 = 0.50 MPa , p02 = 0.0 MPa and v0 = 6.323×10−4 , then after 20000 iterations, the time series of the pore pressure and the flow velocity of the five nodes corresponding to x=0.1H, 0.3H, 0.5H, 0.7H and 0.9H are shown in Figs. 3a and b.
(b) Curve of dimensionless flow velocity vs. time
Fig. 3 Stability analysis of steady-state corresponding to b > 0
It is clear that the pore pressure and flow velocity of each node are all attracted to its steady-state in a stable fashion, so the steady-states corresponding to b > 0 are stable.
4.2
Analysis to stability of steady-state when bs < b < 0
As shown in Fig. 2, when bs < b < 0 the steady-states have two parts, i.e., solutions 2 and 3 with the critical point A (bsˈ− 0.5) as their boundary, where for each parameter b, there are two steady-states of flow velocity, respectively noted as v s 2 and v s 3 . Making use of Eq.(14) and the parameters k , ca , H, p1 , p 2 , μ and ρ 0 , mentioned before, the corresponding critical parameter b can be obtained, i.e., bs = −1.457 × 1015 (kg/m4). Thus, when b changes and satisfies bs < b < 0 , by numerical analysis, we can observe the stability conditions of the steady-states in every branch. If we arbitrarily select the parameter b satisfying bs < b < 0 , for example, b = −1.457 × 1014 (kg/m4), then we obtain the two steady-states of flow velocity corresponding to the solution diagram (Fig. 2), i.e.,
vs2 = −0.02563 and vs3 = −0.97431 . 4.2.1 Stability of steady-states corresponding to the solution 2 When we study the stability of the steady-state vs2 (0.02563), we can discuss two cases in respect of the initial flow velocity v02 , which can either be larger
or smaller than the steady-state velocity vs2 : 1) v02 > vs2 , let v02 = −0.02493 , i.e., the initial velocity deviates from the steady-state and lies on its high side. By evolving, the flow velocity eventually comes to its steady-state vs2 (–0.02563) in a stable way, as shown in Fig. 4, where only the curves of flow state variables vs time at x=0.5H are presented. In particular, it should be pointed out that when b < 0 , the time variables in all the curves below are the magnitude or the absolute value of dimensionless time t . 2) v02 < vs2 , if we let v02 = −0.02573 , by evolving, the flow velocity eventually comes to its steady-state vs2 (–0.02563) in a stable way, as shown in Fig. 5. It is obvious that the steady-states corresponding to solution 2 are stable nodes.
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
(a) Curve of dimensionless pore pressure vs. time
Fig. 4
(b) Dimensionless flow velocity vs. time
Dynamic responses corresponding to v02 < vs2 for branch 2
4.2.2 Stability of steady-states corresponding to solution 3 We continue to analyze the stability of steadystates corresponding to solution 3 and also discuss the stability according to two cases: 1) v03 > vs3 , if we let v03 = −0.9671 , then by evolv-
(a) Dimensionless pore pressure vs. time
Fig. 6
(b) Curve of dimensionless pore pressure vs. time
Dynamic responses corresponding to v02 > vs2 for branch 2
(a) Dimensionless pore pressure vs. time
Fig. 5
5
ing, the flow velocity in the end does not come to its corresponding steady-state vs3 (–0.9743); on the contrary, it is attracted to the steady-state vs2 (–0.02563) corresponding to solution 2, as shown in Fig. 6.
(b) Dimensionless flow velocity vs. time
Dynamic responses corresponding to v03 > vs3 for branch 3
2) v03 < vs3 , if we let v03 = −1.0464 and the time space is the same as above, then by evolving, the flow velocity in the end does not come either to its corresponding steady-state vs3 (–0.9743); on the contrary, it deviates from the steady-state farther and farther and moreover, is attracted to negative infinity as shown in Fig. 7. It is obvious that the steady-states corresponding to solution 3 are unstable. When the initial velocity (dimensionless) is larger than its steady-state, its phase trajectory is, in the end, attracted to the steady-states corresponding to branch 2. When the initial flow velocity is less than its steady-state, its phase trajectory deviates from its steady-state and is attracted to nega-
tive infinity. So, the steady-states corresponding to solution 3 are saddle points. 4.2.3 Analysis to the flow stability when b < bs When b < bs , there are no steady-states. If we select the parameter b satisfying b < bs , for example, b = −1.457 × 1016 < bs = −1.457 × 1015 (kg/m4), then corresponding to a freely given initial velocity v0 (no matter how small), by evolving, the dimensionless flow velocity in an absolute sense tends to negative infinity, i.e., the flow system lost its stability. The dynamic responses corresponding to the given b and v0 = 2.63 × 10−6 (m/s) are shown in Fig. 8.
6
Mining Science and Technology
Vol.19
(a) Dimensionless flow velocity vs. time
Fig. 7
(b) Phase trajectory
Dynamic responses corresponding to v03 < vs3 for branch 3
(a) Dimensionless flow velocity vs. time
Fig. 8
5
No.1
(b) Phase trajectory
Unstable motion of non-darcy flow system when b < bs
Studies on bifurcation and catastrophe of flow in broken rock
From our investigation, we know that when b < 0 , there are two branches of the steady-states. The steady-states corresponding to solution 2 are stable nodes while those corresponding to solution 3 are unstable saddles. When b decreases gradually and tends to the critical parameter bs , the saddle and the
larger). Yet, for the different flow system they all satisfy vs = −0.5 when b = bs . So, we can obtain different bifurcation curves corresponding to the flow systems with different controlling parameter as shown in Fig. 10, where the real line and the dotted line are, respectively, corresponding to stable and unstable flows.
nodes meet and annihilate each other. When b < bs the steady-states vanish, a condition where all the phase-trajectories flow to negative infinity (only one trajectory is drawn in Fig. 8). When bs < b < 0 , the phase trajectories are divided into two flow areas with 4b k 2 ( p1 − p2 − ρ0 gH ) 1 1 branch 3, i.e., vs3 = − − 1 + 2 2 μ2H serving as their boundary, that is to say, for the trajectories corresponding to v0 > vs3 , their flow velocities finally tend to the steady-state vs2 = −
Fig. 9
Saddle-node bifurcation of flow in broken rock
4b k 2 ( p1 − p2 − ρ 0 gH ) 1 1 , while for + 1+ 2 2 μ 2H
the trajectories corresponding to v0 < vs3 ; their flow velocities finally tend to negative infinity. Hence, there is a saddle-node bifurcation (fold bifurcation) at b = bs in the dynamic system of flows in broken rock. The corresponding sketch map of bifurcation is shown in Fig. 9. For different flow systems, the permeability k and critical parameter bs also vary. Moreover, the smaller k, the smaller bs (but its absolute value is
Fig. 10
Bifurcation curve for the flow systems with different controlling parameter
Saddle-node bifurcation (fold bifurcation) is a typical discontinuous bifurcation. Catastrophes of stability of the flow system take place at the bifurcation point. According to the theory of catastrophe[17], a fold catastrophe is apt to occur at the point of fold bifurcation where, at this time, any tiny disturbance can induce some dynamic disasters such as water in-
MIAO Xie-xing et al
Bifurcation and catastrophe of seepage flow system in broken rock
rush or gas bursts. Given the second equation in Eq.(8), the potential function U , with regard to the generalized variable v , should satisfy ∂U ∂p = a2v 2 + a2v + a1 + a4 (15) ∂v ∂x ∂p p2 − p1 where the pore pressure gradient = =C . ∂x p0 Integration of Eq.(15) obtains: 1 1 U = a2v 3 + a2v 2 + a4v + a1Cv (16) 3 2 Developing the first and second order derivatives with respect to v , from Eq.(16), the set S of singularity points satisfies ∂U = a2v 2 + a2v + a1C + a4 = 0 ∂v ∂ 2U = 2 a 2 v + a2 = 0 ∂v 2 Solving these two equations jointly, we obtain the singularity point
μ 2H (17) 4k 2 ( p1 − p2 − ρ0 gH ) which is precisely the bifurcation point obtained earlier. In order to acquire the standard potential function, by coordinate translation, let v = v′ + d (18) Substituting this formula into Eq.(16) and letting d = −1 2 to eliminate the term containing the square of v ′ , we have v = −0.5 , b = −
a 1 U = a2 v ′3 + ( − 2 + a4 + a1C )v ′ 3 4 1 a4 + (− a1C − + c0 ) (19) 2 2 where c0 is the constant of integration. Eq.(19) is the standard potential function with respect to the generalized state variable v ′ for the fold catastrophe of the flow system in broken rock.
6
Conclusions
Given the permeating experimental results of a non-Darcy flow in broken rock, the nonlinear dynamic equations of non-Darcy and non-steady seepage flows in broken rock have been established. By dimensionless transformation, the solution diagram of the steady-states under the given boundary conditions is obtained for the dynamic system of the seepage flow. With a change in the seepage flow parameter b, the steady-states of flow velocity are divided into three parts. Numerical analyses indicate that there is a saddle-node bifurcation (fold bifurcation) at the critical parameter bs for the seepage flow in broken rock. Bifurcation curves of the different seepage flow sys-
7
tems have also been provided. Near the bifurcation point, catastrophe of the stability of the steady-states takes place when, at the same time, the breaking phenomenon of the solid grains of the broken rock is rather serious. So for the seepage flow system in the broken rock catastrophe is apt to occur at the bifurcation point and will result in dynamic disasters such as a water inrush.
References [1]
[2]
[3] [4]
[5] [6] [7] [8] [9]
[10] [11]
[12] [13]
[14] [15]
[16] [17]
Zhu X S. The application of engineering geophysics in the remediation of water inrush in fractured zone. Express Information of Mining Industry, 2004(11): 51–53. (In Chinese) Wang L, Liu C A. Analysis and disposal to the flooded mine resulting from the water inrush in breaking zone of shaft. Mine Construction Technology, 1999, 20(2): 8–11. (In Chinese) Zhou R G, Cheng B F, Ye G J, Wu Q. The effect of water busting in fault rupture zone. Journal of Engineering Geology, 2000, 8(4): 411–415. (In Chinese) Li F G. Analysis and disposal of seepage and water gushing in tunnel fault broken zone. West-China Exploration Engineering, 2004, 132(12): 112–114. (In Chinese) Stephson D. Hydraulic Calculation in Rockfill Engineering. Beijing: Ocean Publishing Company, 1984. Wang Z F. Notes of rockfill seepage research. Journal of Nanjing Hydraulic Research Institute, 1990(4): 445–452. (In Chinese) Bear J. Dynamics of Fluid in Porous Media. New York: Elsevier, 1979. Li S C. Nonlinear Dynamical Study on Non-darcy Flow in Broken Rock [Ph.D. dissertation]. Xuzhou: China University of Mining & Technology, 2006. (In Chinese) Miao X X, Li S C, Huang X W. Experimental study of seepage properties of non-Darcy flow in granular gangues. Journal of University of Mining and Technology, 2006, 16(2): 105–109. (In Chinese) Deng Y E, Liu C Q, Huang R Q. Advanced Seepage Theory and Method. Beijing: Scientific Publishing Company, 2004. (In Chinese) Chen Z Q. Bifurcation of Dynamic System of Non-Darcy Flow in Post-Failure Rock [Ph.D. dissertation]. Xuzhou: China University of Mining & Technology, 2003. (In Chinese) Li T Z, Li Y S, Ma Z G. Testing study of non-Darcy seepage flow in fractured rocks. Engineering Mechanics, 2003, 20(4): 132–135. (In Chinese) Tang P, Huang X W, Li T Z. Test study on the penetrating properties of non-Darcy flow in rock. Journal of Xiangtan Normal University, 2002, 24(2): 34–38. (In Chinese) Miao X X, Chen Z Q, Mao X B, Chen R H. The bifurcation of non-Darcy flow in post-failure rock. Acta Mechanica Sinica, 2003, 35(6): 660–667. (In Chinese) Pradip K G N, Venkataraman P. Non-Darcy converging flow through coarse granular media. Journal of the Institution of Engineers (India), Civil Engineering Division, 1995(76): 6–11. Kong X Y. Advanced Mechanics of Fluid in Porous Media. Hefei: Press of China University of Science and Technology, 1999. (In Chinese) Liu B Z, Peng J H. Nonlinear Dynamics. Beijing: Higher Education Press, 2003. (In Chinese)