基于darcy方程的突水流动数值模拟

基于darcy方程的突水流动数值模拟

1采动岩体渗流数值模型研究采井运动极大地改变了围岩破坏的渗透性，导致顶板、断裂带和底板突然事故，这是矿山生产经济损失最大的事故。开展采动条件下岩体突水机制及渗流数值模型研究，对于采动岩体突水预测与防治、开采方法的改进以及安全度的评价具有重大理论意义和实际价值。开采过程中岩体破坏突水机制十分复杂，基于渗流力学理论建立适合采动岩体渗流突水力学模型与数值计算方法，是解决问题的技术关键。根据采动岩体破坏非线性渗流的特点，本文提出采用Brinkman方程，研究水在破碎岩体中的流动规律，探索含水层不同条件对破碎岩体水渗流的作用机制，为正确预测突水水量和压力提供科学依据。2darcy方程的应用针对矿山生产过程中断层突水、顶底板突水和岩溶陷落柱突水预测与防治问题，很多学者[1,2,3,4,5,6,7,8,9,10]开展了大量的研究工作。对于岩层突水问题，采动应力和水压力作用下岩层破断渗透性及其渗透压力作用机制是突水机制研究的核心问题。岩体采动破坏(尤其是断层破碎带和陷落柱)后裂隙尺度和流速急剧增大，Reynolds数远远大于10，呈现非Darcy渗流特性；同时水流在导水裂隙中快速运动使得水压作用力呈现惯性力和渗透动压两种方式。W.F.Brace在关于岩石破坏渗透性演化规律的综述中也着重强调，随着变形的增加，岩性和孔隙、微裂隙结构等因素对应力–渗透性关系的影响比较复杂，表现在对岩石峰后渗透增大幅度的量化描述上。采动条件下，无论是陷落柱、断层破碎带、围岩破坏区域，都由破碎岩体组成，属于大空隙的多孔介质，渗流通道系统比较复杂。目前，针对碎裂岩体的渗流场研究比较少，一般来讲，描述流体运动的流场方程包括3种：(1)以线性层流为主、忽略流体惯性力的Darcy方程，水流在含水层中的流动符合该方程；(2)Navier-Stokes方程，忽略流体渗流阻力，突水后水流在巷道内流动符合该方程；(3)介于Darcy方程和Navier-Stokes方程之间考虑非线性渗流性态的方程。理论上讲，破碎带就是联系含水层(Darcy层流)和巷道(自由管流)之间的过渡区域。Darcy方程以流体压力驱动为主，适合低渗透多孔介质，地下水在含水层中渗流采用该方法计算。针对采动岩体破碎带，作者在Darcy方程的基础上引入渗透率突跳系数的概念，认为渗流模型中单元破坏后渗透系数在原来的基础上增大若干倍，具体增大的数值可依据室内实验和现场水文实验结果。该方法的优点是便于数值求解，但还是基于Darcy方程求解，不能描述非Darcy效应。S.C.Yuan和J.P.Harrison对峰后破碎岩石采用立方定律和strainpartition技术建立渗流–体应变关系方程和数值模型，该模型能够定量描述峰后破碎岩石渗透性的急剧增大，但归根到底，该模型也属于Darcy方程。Darcy方程可表示为式中：u为流体流速(m/s)；k为渗透率(m2)；η为动黏系数(Pa·s)；p为流体压力(Pa)；Z为位置高度(m)；ρ为流体密度(kg/m3)；g为重力加速度，一般取9.81m/s2。Navier-Stokes方程基于牛顿第二定律，刻画流体在重力、黏性阻力和压力的作用下运动规律，考虑了流体静压能、动能和势能平衡，以流体动能为主，不考虑渗透阻力的作用，主要研究管流，适合河道、管道流场，在巷道通风、流体管流计算中得到广泛应用。显然不适合描述水在破碎岩体中渗流。Navier-Stokes方程可表示为式中：I为单位矩阵。非Darcy方程在土石坝或堆积体水流渗流过程中得到应用，对于破碎岩体，服从AhmedSunada(Forcheimer)关系，渗流系统是非线性的。Ahmed-Sunada型非Darcy孔隙渗流系统的控制方程系统的行为由Reynolds数和Darcy数2个参数调节，E.Choi等的研究结果表明：Brinkman黏性项影响流场的分布。缪协兴等提出了破碎岩体非Darcy方程，考虑了流体的惯性作用。由于非Darcy方程很难求解，即使一维模型，参数变化引起方程求解很不稳定，出现分叉现象，所以缪协兴等对于突水计算没有能够进行数值求解。流体在介质中流动，当速度足够快，以至于由于剪切作用引起的能量耗散不能忽略，需要考虑黏性流体的剪切应力。H.C.Brinkman在Darcy方程的基础上考虑Navier-Stokes方程中流体黏性剪切应力项，提出了Brinkman方程，该方程基于牛顿第二定律，刻画流体在孔隙介质中快速运动形成的剪切力、渗透压力作用下的运动规律，适合于表述孔隙介质中的Darcy流与流体管流之间的过渡区域，该方程在多孔介质内部流动的非Darcy效应理论分析和实际工程问题模拟中具有较好的效果，比较适合采动破碎岩体或陷落柱非Darcy快速渗流特点。Brinkman方程可表示为3巷道物理过程模拟图1所示为陷落柱突水概念模型，由图可知，当巷道掘进到陷落柱时诱发突水：含水层中的高承压(可达几个兆帕)水通过陷落柱快速流动进入巷道之中，在高渗流压力作用下，突水过程中水流经历了在含水层中的Darcy层流、陷落柱或断层带中的Brinkman快速流以及在巷道中的Navier-Stokes紊流3个物理过程。针对上述概念模型，在本文的计算模型中，联立式(1)～(3)，可以解出突水过程中水流从含水层经过导水破坏带进入到巷道发生突水的全过程。可见本文采用Brinkman方程针对破碎带水流要兼顾流体压力梯度和运动作用的特点，比较适合破碎带渗流运动场，同时可以把含水层和巷道整个水流路径连接在一起。应用上述计算方程，采用FEMLAB系统进行求解，研究突水动态过程。FEMLAB系统(也称为COMSOLMultiphysics系统)是基于偏微分方程组(PDEs)而开发的多物理场耦合过程分析工具，偏微分方程是描述科学规律的基础，应用该工具可将任意耦合偏微分方程转化为适当的形式以便于数值分析，并运用基于有限元方法的高效求解器进行求解。本文在该系统平台上进行二次开发，基于质量守恒和压力平衡，把Brinkman方程与Navier-Stokes方程和Darcy方程耦合在一起，确定适合现场实际的边界条件和初始条件，求解突水问题。4计算4.1含水层流场的边界条件如图1所示，流体在岩溶含水层区域为Darcy流动，进入陷落柱破碎带后渐变为Brinkman流动，进入巷道中完成Navier-Stokes流动过程。假设流体密度和黏度都是常数，流动是稳态过程，流体进入含水层的流量和流出巷道的压力已知，则流体在3种流动区域过渡时的边界条件如下：(1)含水层的Darcy流动：在含水层中的水流满足线性Darcy定律，设pesdl为Darcy区域的流体压力，对于稳定流，式(1)可表示为在满足Darcy定律的含水层和满足Brinkman方程的陷落柱交界面上，要满足压力连续条件，即式中：pchns2为Brinkman区域的流体压力。设流体进入含水层边界的单宽法线流速uesdln=nuesdl(n为含水层边界处的法线单位矢量)，则可求得uesdln=1×10-3m3/(s·m)，在含水层上下两侧为隔水边界，流体的Darcy速度为0，则整个含水层的边界条件为(2)破碎带Brinkman流动：Brinkman方程描述了介于Darcy流动和Navier-Stokes流动之间的一种过渡流动状态。破碎带中流体流动速度比较大，剪应力引起的动量传递比较大，剪应力效应不能忽略，需要采用Brinkman方程来描述(见式(3))。在本模型中，忽略F项的影响，边界条件定义如下：在Brinkman方程中，速度是独立变量，所以在含水层Darcy流动和破碎带Brinkman流动的交界面上，要满足速度连续条件，即式中：uesdl,uchns2分别为流体在含水层和破碎带的流速。在破碎带Brinkman流动和巷道工作面N-S流动交界面上，要满足压力连续条件，即在破碎带其他边界为隔水边界，即(3)水进入巷道的Navier-Stokes流动：突水水流进入巷道后为自由流动，由不可压缩流体的NavierStokes方程描述(见式(2))，边界条件定义如下：在Navier-Stokes和Brinkman流动的交界面上，即在巷道与陷落柱交界处，要满足速度连续条件，即式中：uchns为流体在巷道中的流速。对于巷道出口处的压力p0是已知的，一般等于大气压力0.1MPa，即在巷道其他边界为隔水边界，即联立式(4)～(12)的边界条件等式，就可以保证流体质量守恒和压力平衡，通过数值计算，耦合求解含水层中的Darcy层流、陷落柱或断层带中的Brinkman快速流以及在巷道中的Navier-Stokes紊流3个物理过程。4.2采动诱导坍塌柱突水模拟4.2.1底部含水层边界根据图1的陷落柱突水概念模型，建立了如图2所示的陷落柱突水计算模型。设模型长40m，高22m，其中底部含水层长40m，高6m，中部陷落柱平均高14m，平均宽5m，右上侧巷道平均长18m，高2m。底部含水层的边界为进水流量边界，单宽流速为：uesdln=1×10-3m3/(s·m)，相当于2.6MPa的承压水压。上部巷道出口的边界为流体压力边界，等于大气压力。含水层中流体的黏滞系数为1×10-3Pa·s，先假设含水层和陷落柱的渗透率相同，为1×10-12m2，分析Brinkman区域和Darcy区域同样渗透系数但不同渗流方程情况下Brinkman方程引起流速的增大幅度，然后把陷落柱Brinkman区域的渗透率提高10倍，分别研究当含水层补给水量不变情况下高渗透率的陷落柱引起的水压卸压幅度和当保持含水层水压不变情况下的流速急剧增大的幅度，进一步揭示陷落柱非Darcy渗流对突水过程的作用机制。4.2.2brifman区域渗透率计算结果如图3～6所示。图3,4分别显示了流体从Darcy区域到Brinkman区域再到巷道中Navier-Stoke流动的水流速度和压力的变化。图5,6分别为图4中沿折线A1—A2—A3—A4方向切线流速和压力的分布曲线(图中：工况1表示Darcy区域和Brinkman区域渗透率一致；工况2表示Brinkman区域渗透率提高10倍，但边界流量和工况1一致；工况3表示Brinkman区域渗透率提高10倍，且保持边界恒定高压力)。陷落柱突水渗流模拟结果可见，速度变化(表面和箭头)和压力分布(表面)沿3个流域是连续变化的，流速梯度的大幅度变化主要集中在Brinkman区域，压力梯度变化则主要在Darcy区域，表明水流在Darcy区域内的渗流阻力最大。图5,6表明，当Darcy区域和Brinkman区域采用相同渗透率时，虽然在Darcy-Brinkman交界面和BrinkmanNavier-Stokes交界面上流速急剧增大，但过渡比较连续。在Darcy区域，流速为2.0×10-5m/s，接近到Brinkman区域，流速逐渐增大，在边界面流速为7.0×10-5m/s，即使2个区域的渗透系数完全相等，由于Brinkman考虑了流体动能的耗散作用，在Brinkman区域流速也急剧增大，进入到Navier-Stokes区域的边界面，流速达到2.3×10-4m/s，增大了10倍，进入到Navier-Stokes区域后，流速进一步增大，达到3.5×10-4m/s，增大了17倍。而基于线性渗流理论，整个区域只用Darcy方程计算时，流速不会有这么急剧的阶跃变化，无法描述突水时流量的非线性变化过程。而当将Brinkman区域渗透率提高10倍时，该条件下若保持补给水量不变，将造成边界水压力从2.6MPa降低到0.7MPa，降低了2.7倍，表明当含水层补给水量一定时，Brinkman区域渗透性的提高会起到较大幅度的卸压作用，从而降低了压力梯度，Brinkman区域渗透性的提高对水流速度变化没有影响，进一步研究表明，揭穿陷落柱时，若含水层补给水量有限，即使陷落柱渗透性很高，也不会发生突水灾害。当含水层高水压保持不变时，同样的水压梯度(水流进入巷道中降为大气压)下，Brinkman区域渗透率的增大引起进入Brinkman和Navier-Stokes区域水流流速更大幅度地增高。在Darcy区域，流速为1.4×10-5m/s，接近到Brinkman区域，流速急剧增大，在边界面流速为5.1×10-4m/s，流速增大幅度出现“阶越”式变化，进入到Navier-Stokes区域的边界面，流速达到1.2×10-3m/s，增大了85倍，进入到Navier-Stokes区域后，流速进一步增大，达到1.6×10-3m/s，增大了114倍。由此可见，陷落柱作为含水层渗流和巷道突水自由流动的过渡区域，其渗透性变化对于突水压力和流速演变作用十分显著，这进一步表明，现场实际的突水现象就是在含水层充足的补给水量和保持恒定的高水压作用下，当陷落柱或导水破碎带联通了含水层和巷道情况下形成的。5巷道渗流场模拟结果分析本文提出突水过程中含水层、破碎导水带和巷道3种水流运移耦合模型，采用Brinkman方程计算破碎带水流场，比较符合破碎岩体介质的流场特性，同时把含水层Darcy渗流和巷道Navier-Stokes流动有机联系在统一流动场中，得到整个水流流场分布，更符合突水实际。采用非线性耦合渗流