从岩层渗流失稳的视角认识煤矿突水机理

从岩层渗流失稳的视角认识煤矿突水机理*孙明贵1李天珍2黄先伍3，4岳建华3（1安徽理工大学淮南232000，2徐州工程学院徐州221008，3中国矿业大学徐州221008，4贵州工业大学贵阳550004）E-mail:[e-mail]摘要：将岩体描述成由多层不同厚度和不同渗透特性岩石构成的多孔介质，建立了岩体非Darcy渗流系统的非线性动力学模型，讨论了系统的平衡态的存在性，用分离变量法建立了系统的失稳条件。从渗流系统失稳的视角揭示了矿井突水的力学机理，认为突水是层状岩体渗流失稳的过程。研究表明，矿井突水是由各层岩石的厚度、渗透特性、渗流液体的力学性质及边界压力共同决定的。当渗流系统的控制参量与边界压力满足一定的条件时，系统失去稳定性，引发突水灾害；当岩层中某层岩石的渗透率比其他岩层低一个量级，而非Darcy流β因子比其他岩层高两个量级时，该岩层为渗流系统的关键层，其渗透特性对系统的稳定性起支配作用。关键词：煤矿突水非Darcy流渗透特性稳定性1引言突水是煤矿重大灾害之一。近40年来，我国共发生矿井突水事故2000余次，有2000余人丧生，直接经济损失高达40亿元。据统计，我国国有煤矿中半数煤矿存在突水危险性，并且突水危险越来越严重[1]。我国学者在突水机理研究和突水灾害防治技术攻关中取得举世瞩目的成就，代表了当前采矿技术领域的世界领先水平。为探索矿井突水机理，先后提出多种力学模型，主要有：毕贤顺、王晋平的流-固耦合模型[2]，桂和荣、龚乃勤、孙本魁的渗流-应力-温度多场耦合模型[3]，施龙青、尹增德、刘永法的损伤模型[4]，勒德武、王延福、马培智的塞子模型[5]等。而且，当代非线性科学的研究成果在突水防治中得到了应用[6-8]。然而，这些研究工作中突水与岩体的几何特性和渗流特性联系很少。缪协兴、陈占清等试图从峰后岩石非Darcy流的失稳过程解释突水机理，但其模型与工程实际相差很大[9-11]。本文将煤矿围岩视为由多岩层组成的孔隙介质，各岩层的矿物成分不同、应力-应变状态不同决定了各岩层的渗透特性（渗透率、非Darcy流β因子河和加速度系数）不同，水在岩层中的渗流失稳引发突水。2层状岩体渗流的非线性动力学模型设围岩由层岩石组成，各层的厚度为，孔隙率为，渗透率为，非Darcy流nihiφikβ因子为，加速度系数为，两端的压力为和，iβiac1−ipip),,,(niL21=，见图1。岩层渗流的动力学方程由质量方程、动量方程、孔隙压缩方程和流体压缩方程组成，对于一维无源流动，当忽略体积力时，这些方程分别为[12]*国家杰出青年基金(50225414),国家自然科学基金重点项目(50134040)和教育部博士点基金(20010290003)资助。_______________________________________________________________________________⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫+=+=+−∂∂−=∂∂=∂∂+∂∂)1()1(0)()(002pcpcVVkxptVcVxtfiiiiiiaiρρφφρβµρρρφφ(（1）),,,ihip1−ipxWnp0p岩层渗流的动力学模型图1h1−iHiHniL21=其中：V为渗流速度；为流体的相对压力；pρ和µ分别为水的质量密度和动力粘度；，φ分别为参考压力下对应的流体质量密度和第层岩石的孔隙率；为流体的等温压缩系数；为第i层岩石的孔隙压缩系数。0ρi00pifcicφ假设，，则孔隙压缩方程与流体压缩方程可以合写为11pcpcifφ,),,,(niL21=tpctitii∂∂=∂∂00φρρφ)((),,,niL21=（2）其中为第层岩石的综合压缩系数。这样，控制方程可化为ifitcccφ+=),,,(niL21=i⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧∂∂++−+−∂∂−=∂∂∂∂−=∂∂xWWpcccWWkpcxpctWxWctpfitfiifiaiti)()]()([11111002000φρβµρφρ（3）),,,),,(),,[niHHxtiiL2101=∈+∞∈−其中；W表示动量密度，即),,,(,∑===ijjinihH121LVWρ=。当忽略水的压缩性时，式（3）可以简化成⎪⎪⎩⎪⎪⎨⎧−−∂∂−=∂∂∂∂−=∂∂)]([2000111WWkxpctWxWctpiiiaitiβµρφρ（4）),,,),,(),,[niHHxtiiL2101=∈+∞∈−方程（4）即为岩层渗流系统的动力学模型。动力系统（4）的平衡态由方程（5）确定⎪⎪⎩⎪⎪⎨⎧=+−+∂∂−=∂∂=∂∂=∂∂0100020)(:)(:sisissWWkxptWxWtpβµρ),,,(niL21=（5）可见，在平衡态各层岩石的动量密度成均匀分布，从而压力梯度为常量。由式（5）直接积分，可以得到2中国科技论文在线_______________________________________________________________________________⎪⎩⎪⎨⎧−=−=−)(201sisiiiisWWkhppconstWβµρ),,,(niL21=（6）将（6）中第二式对i求和，得到2sspBWAWG−=（7）其中∑==niiihkhA10ρµ，∑==niiihhB10ρβ，hppGnp−=0，（8）∑===niinhHh1h为岩层的总厚度，为岩层的平均压力梯度，pGA和是反映岩层渗透特性的综合参量。B由（7）可以看出，层状岩体渗流系统平衡态的存在性由判别式pBGAD42−=（9）决定，即当时，系统存在两平衡态0D⎪⎪⎩⎪⎪⎨⎧=∈−−−=−−=−−−−niHHxHxhppppBBGAAWiiiiiiiss,,,),,(),(L2124111121（10）⎪⎪⎩⎪⎪⎨⎧=∈−−−=−+=−−−−niHHxHxhppppBBGAAWiiiiiiiss,,,),,(),(L2124111122（11）当时，系统只有一个平衡态0=D⎪⎪⎩⎪⎪⎨⎧=∈−−−=−=−−−−niHHxHxhppppBAWiiiiiiiss,,,),,(),(L2121111（12）当时，系统不存在平衡态。这时，不论系统从什么初始条件演化，都不会达到稳定状态，即系统是不稳定的。因此，0D0D是系统的失稳条件之一。引入无量纲的量thBAt~0ρ=，hxx=~，phGpp=~，WABW=~（13）则有tABht~0ρ=，xhx~=，phGpp~=，WBAW~=（14）3中国科技论文在线_______________________________________________________________________________于是方程组（4）可以改写成⎪⎪⎩⎪⎪⎨⎧+−∂∂−=∂∂∂∂−=∂∂23410WaWaxpatWxWatpiiii~~~~~~~~~~),,,),,(~),,[~nihHhHxtiiL2101=∈+∞∈−（15）其中hGcapitii001φ=,iapichGa01ρ=，BkchAaiiaiµ=4，223BchAaiaiiβ=（16）以无量纲变量p~，W~表示的系统平衡态为niHHxhHxhGpphGppiiiipiipis,,,),,(),~(~L211111=∈−−−=−−−−（17a）⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧−=−−⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−±=时当时当时当041none0412104141121~2222ABGABGABGABGWpppps（17b）现在考虑系统在平衡态附近的情况，令)~,~()(~~xtpxpps1+=，)~,~(~~xt=（18）将式（19）代入式（15），得到⎪⎪⎩⎪⎪⎨⎧++∂∂−=∂∂∂∂−=∂∂21312111101WaWaxpatWxWatpiiii~~~~),,,),,(~),,[~nihHhHxtiiL2101=∈+∞∈−（19）其中isiiaWaa4322−=~（20）式（19）即以无量纲量表示的层状岩体非Darcy渗流系统在平衡态附近的演化方程，当系统不存在平衡态时，可在（20）中令0=sW~。3系统稳定性分析引入函数η，使tW~~∂∂=η，xapi~~0∂∂−=η（21）4中国科技论文在线_______________________________________________________________________________这样，方程（21）的第一式自然满足，而第二式成为234221022~~~~⎟⎠⎞⎜⎝⎛∂∂+∂∂−∂∂=∂∂tataxaatiiiiηηηη，),,,),,(~),,[~nihHhHxtiiL2101=∈+∞∈−（22）设（22）具有分离变量形式的解，即)()~(2~1)()(2111thHxhGppxhGpatxiipiipiiςςξη+⎥⎥⎦⎤⎢⎢⎣⎡−−−−=+=−−−0，),,,),,(~),,[~nihHhHxtiiL2101=∈+∞∈−（23）代（23）入（22），得到23411ςςς&&&&iiipiiiaahGppa+−−=−，),,2,1),,(~),,0[~1nihHhHxtiiL=∈+∞∈−（24）由于ς在),(~1hHhHxii−∈上任意点数值相等，并且在任意瞬时t，式（24）右端函数应该连续，这样在任意瞬时t，有∑∑∑===−−+−−=+−−==+−−=+−−niniiniiipiiinnnpnnnppananhGppanaahGppaaahGppaaahGppa121314112341122324221212131411011111ςςςςςςςς&&&&L&&&&)()(（25）于是式样（24）便可简化成∑∑===−+−−=niniiniiipiiiananhGppan12131411111ςςς&&&&)()(∑（26）对方程（26）积分，得到∫∫∑∑∑=⎟⎟⎠⎞⎜⎜⎝⎛+⎟⎟⎠⎞⎜⎜⎝⎛−−===−dthhanhhanhGppandiniiiniiniipiii21314111111ςςς&&&（27）将式（27）的左端积分函数的分母配方，得到∫=+−dtd3221λλλςς)(&&∫（28）其中0213141⎟⎟⎠⎞⎜⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛=∑∑==hhahhainiiiniiλ，213141311122⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛−⎟⎟⎠⎞⎜⎜⎝⎛−=∑∑∑∑====−hhahhahhahGppainiiiniiiniiniipiiiλ，01133=∑=hhaniniiλ（29）当时，可将写成，这种情形下方程（28）的解为02λ2λ242λλ=5中国科技论文在线_______________________________________________________________________________))(tan(~13441CtW++==λλλλς&（30）这里为积分常数。可见时，系统失稳。而当时，可将写成，这种情形下方程（28）的解为1C02λ02λ2λ242λλ−=))(tanh(~23441CtW++==λλλλς&（31）这里为积分常数。可见时，系统是稳定的。这时，系统存在两个平衡态2C02λ2141lim~λλλλ−±=+==+∞→WWts当时，方程（28）的解为02=λ3311~CtW+−==λλς&（32）其中是积分常数，可见时，系统的稳定性取决于，即取决于初始条件。当时，系统是稳定的。这时，系统存在唯一的平衡态3C02=λ3C03C1lim~λ==+∞→WWts基于以上讨论，是系统的分岔条件。将式（16）代入式（29），得到02=λ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛−−⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎝⎛=∑∑∑∑====2111002142)(42niiiaininjjajjianiniik