基于darcy方程的突水流動數(shù)值模擬

基于darcy方程的突水流動數(shù)值模擬

1采動巖體滲流數(shù)值模型研究采井運(yùn)動極大地改變了圍巖破壞的滲透性，導(dǎo)致頂板、斷裂帶和底板突然事故，這是礦山生產(chǎn)經(jīng)濟(jì)損失最大的事故。開展采動條件下巖體突水機(jī)制及滲流數(shù)值模型研究，對于采動巖體突水預(yù)測與防治、開采方法的改進(jìn)以及安全度的評價(jià)具有重大理論意義和實(shí)際價(jià)值。開采過程中巖體破壞突水機(jī)制十分復(fù)雜，基于滲流力學(xué)理論建立適合采動巖體滲流突水力學(xué)模型與數(shù)值計(jì)算方法，是解決問題的技術(shù)關(guān)鍵。根據(jù)采動巖體破壞非線性滲流的特點(diǎn)，本文提出采用Brinkman方程，研究水在破碎巖體中的流動規(guī)律，探索含水層不同條件對破碎巖體水滲流的作用機(jī)制，為正確預(yù)測突水水量和壓力提供科學(xué)依據(jù)。2darcy方程的應(yīng)用針對礦山生產(chǎn)過程中斷層突水、頂?shù)装逋凰蛶r溶陷落柱突水預(yù)測與防治問題，很多學(xué)者[1,2,3,4,5,6,7,8,9,10]開展了大量的研究工作。對于巖層突水問題，采動應(yīng)力和水壓力作用下巖層破斷滲透性及其滲透壓力作用機(jī)制是突水機(jī)制研究的核心問題。巖體采動破壞(尤其是斷層破碎帶和陷落柱)后裂隙尺度和流速急劇增大，Reynolds數(shù)遠(yuǎn)遠(yuǎn)大于10，呈現(xiàn)非Darcy滲流特性；同時(shí)水流在導(dǎo)水裂隙中快速運(yùn)動使得水壓作用力呈現(xiàn)慣性力和滲透動壓兩種方式。W.F.Brace在關(guān)于巖石破壞滲透性演化規(guī)律的綜述中也著重強(qiáng)調(diào)，隨著變形的增加，巖性和孔隙、微裂隙結(jié)構(gòu)等因素對應(yīng)力–滲透性關(guān)系的影響比較復(fù)雜，表現(xiàn)在對巖石峰后滲透增大幅度的量化描述上。采動條件下，無論是陷落柱、斷層破碎帶、圍巖破壞區(qū)域，都由破碎巖體組成，屬于大空隙的多孔介質(zhì)，滲流通道系統(tǒng)比較復(fù)雜。目前，針對碎裂巖體的滲流場研究比較少，一般來講，描述流體運(yùn)動的流場方程包括3種：(1)以線性層流為主、忽略流體慣性力的Darcy方程，水流在含水層中的流動符合該方程；(2)Navier-Stokes方程，忽略流體滲流阻力，突水后水流在巷道內(nèi)流動符合該方程；(3)介于Darcy方程和Navier-Stokes方程之間考慮非線性滲流性態(tài)的方程。理論上講，破碎帶就是聯(lián)系含水層(Darcy層流)和巷道(自由管流)之間的過渡區(qū)域。Darcy方程以流體壓力驅(qū)動為主，適合低滲透多孔介質(zhì)，地下水在含水層中滲流采用該方法計(jì)算。針對采動巖體破碎帶，作者在Darcy方程的基礎(chǔ)上引入滲透率突跳系數(shù)的概念，認(rèn)為滲流模型中單元破壞后滲透系數(shù)在原來的基礎(chǔ)上增大若干倍，具體增大的數(shù)值可依據(jù)室內(nèi)實(shí)驗(yàn)和現(xiàn)場水文實(shí)驗(yàn)結(jié)果。該方法的優(yōu)點(diǎn)是便于數(shù)值求解，但還是基于Darcy方程求解，不能描述非Darcy效應(yīng)。S.C.Yuan和J.P.Harrison對峰后破碎巖石采用立方定律和strainpartition技術(shù)建立滲流–體應(yīng)變關(guān)系方程和數(shù)值模型，該模型能夠定量描述峰后破碎巖石滲透性的急劇增大，但歸根到底，該模型也屬于Darcy方程。Darcy方程可表示為式中：u為流體流速(m/s)；k為滲透率(m2)；η為動黏系數(shù)(Pa·s)；p為流體壓力(Pa)；Z為位置高度(m)；ρ為流體密度(kg/m3)；g為重力加速度，一般取9.81m/s2。Navier-Stokes方程基于牛頓第二定律，刻畫流體在重力、黏性阻力和壓力的作用下運(yùn)動規(guī)律，考慮了流體靜壓能、動能和勢能平衡，以流體動能為主，不考慮滲透阻力的作用，主要研究管流，適合河道、管道流場，在巷道通風(fēng)、流體管流計(jì)算中得到廣泛應(yīng)用。顯然不適合描述水在破碎巖體中滲流。Navier-Stokes方程可表示為式中：I為單位矩陣。非Darcy方程在土石壩或堆積體水流滲流過程中得到應(yīng)用，對于破碎巖體，服從AhmedSunada(Forcheimer)關(guān)系，滲流系統(tǒng)是非線性的。Ahmed-Sunada型非Darcy孔隙滲流系統(tǒng)的控制方程系統(tǒng)的行為由Reynolds數(shù)和Darcy數(shù)2個(gè)參數(shù)調(diào)節(jié)，E.Choi等的研究結(jié)果表明：Brinkman黏性項(xiàng)影響流場的分布?？妳f(xié)興等提出了破碎巖體非Darcy方程，考慮了流體的慣性作用。由于非Darcy方程很難求解，即使一維模型，參數(shù)變化引起方程求解很不穩(wěn)定，出現(xiàn)分叉現(xiàn)象，所以繆協(xié)興等對于突水計(jì)算沒有能夠進(jìn)行數(shù)值求解。流體在介質(zhì)中流動，當(dāng)速度足夠快，以至于由于剪切作用引起的能量耗散不能忽略，需要考慮黏性流體的剪切應(yīng)力。H.C.Brinkman在Darcy方程的基礎(chǔ)上考慮Navier-Stokes方程中流體黏性剪切應(yīng)力項(xiàng)，提出了Brinkman方程，該方程基于牛頓第二定律，刻畫流體在孔隙介質(zhì)中快速運(yùn)動形成的剪切力、滲透壓力作用下的運(yùn)動規(guī)律，適合于表述孔隙介質(zhì)中的Darcy流與流體管流之間的過渡區(qū)域，該方程在多孔介質(zhì)內(nèi)部流動的非Darcy效應(yīng)理論分析和實(shí)際工程問題模擬中具有較好的效果，比較適合采動破碎巖體或陷落柱非Darcy快速滲流特點(diǎn)。Brinkman方程可表示為3巷道物理過程模擬圖1所示為陷落柱突水概念模型，由圖可知，當(dāng)巷道掘進(jìn)到陷落柱時(shí)誘發(fā)突水：含水層中的高承壓(可達(dá)幾個(gè)兆帕)水通過陷落柱快速流動進(jìn)入巷道之中，在高滲流壓力作用下，突水過程中水流經(jīng)歷了在含水層中的Darcy層流、陷落柱或斷層帶中的Brinkman快速流以及在巷道中的Navier-Stokes紊流3個(gè)物理過程。針對上述概念模型，在本文的計(jì)算模型中，聯(lián)立式(1)～(3)，可以解出突水過程中水流從含水層經(jīng)過導(dǎo)水破壞帶進(jìn)入到巷道發(fā)生突水的全過程?？梢姳疚牟捎肂rinkman方程針對破碎帶水流要兼顧流體壓力梯度和運(yùn)動作用的特點(diǎn)，比較適合破碎帶滲流運(yùn)動場，同時(shí)可以把含水層和巷道整個(gè)水流路徑連接在一起。應(yīng)用上述計(jì)算方程，采用FEMLAB系統(tǒng)進(jìn)行求解，研究突水動態(tài)過程。FEMLAB系統(tǒng)(也稱為COMSOLMultiphysics系統(tǒng))是基于偏微分方程組(PDEs)而開發(fā)的多物理場耦合過程分析工具，偏微分方程是描述科學(xué)規(guī)律的基礎(chǔ)，應(yīng)用該工具可將任意耦合偏微分方程轉(zhuǎn)化為適當(dāng)?shù)男问揭员阌跀?shù)值分析，并運(yùn)用基于有限元方法的高效求解器進(jìn)行求解。本文在該系統(tǒng)平臺上進(jìn)行二次開發(fā)，基于質(zhì)量守恒和壓力平衡，把Brinkman方程與Navier-Stokes方程和Darcy方程耦合在一起，確定適合現(xiàn)場實(shí)際的邊界條件和初始條件，求解突水問題。4計(jì)算4.1含水層流場的邊界條件如圖1所示，流體在巖溶含水層區(qū)域?yàn)镈arcy流動，進(jìn)入陷落柱破碎帶后漸變?yōu)锽rinkman流動，進(jìn)入巷道中完成Navier-Stokes流動過程。假設(shè)流體密度和黏度都是常數(shù)，流動是穩(wěn)態(tài)過程，流體進(jìn)入含水層的流量和流出巷道的壓力已知，則流體在3種流動區(qū)域過渡時(shí)的邊界條件如下：(1)含水層的Darcy流動：在含水層中的水流滿足線性Darcy定律，設(shè)pesdl為Darcy區(qū)域的流體壓力，對于穩(wěn)定流，式(1)可表示為在滿足Darcy定律的含水層和滿足Brinkman方程的陷落柱交界面上，要滿足壓力連續(xù)條件，即式中：pchns2為Brinkman區(qū)域的流體壓力。設(shè)流體進(jìn)入含水層邊界的單寬法線流速uesdln=nuesdl(n為含水層邊界處的法線單位矢量)，則可求得uesdln=1×10-3m3/(s·m)，在含水層上下兩側(cè)為隔水邊界，流體的Darcy速度為0，則整個(gè)含水層的邊界條件為(2)破碎帶Brinkman流動：Brinkman方程描述了介于Darcy流動和Navier-Stokes流動之間的一種過渡流動狀態(tài)。破碎帶中流體流動速度比較大，剪應(yīng)力引起的動量傳遞比較大，剪應(yīng)力效應(yīng)不能忽略，需要采用Brinkman方程來描述(見式(3))。在本模型中，忽略F項(xiàng)的影響，邊界條件定義如下：在Brinkman方程中，速度是獨(dú)立變量，所以在含水層Darcy流動和破碎帶Brinkman流動的交界面上，要滿足速度連續(xù)條件，即式中：uesdl,uchns2分別為流體在含水層和破碎帶的流速。在破碎帶Brinkman流動和巷道工作面N-S流動交界面上，要滿足壓力連續(xù)條件，即在破碎帶其他邊界為隔水邊界，即(3)水進(jìn)入巷道的Navier-Stokes流動：突水水流進(jìn)入巷道后為自由流動，由不可壓縮流體的NavierStokes方程描述(見式(2))，邊界條件定義如下：在Navier-Stokes和Brinkman流動的交界面上，即在巷道與陷落柱交界處，要滿足速度連續(xù)條件，即式中：uchns為流體在巷道中的流速。對于巷道出口處的壓力p0是已知的，一般等于大氣壓力0.1MPa，即在巷道其他邊界為隔水邊界，即聯(lián)立式(4)～(12)的邊界條件等式，就可以保證流體質(zhì)量守恒和壓力平衡，通過數(shù)值計(jì)算，耦合求解含水層中的Darcy層流、陷落柱或斷層帶中的Brinkman快速流以及在巷道中的Navier-Stokes紊流3個(gè)物理過程。4.2采動誘導(dǎo)坍塌柱突水模擬4.2.1底部含水層邊界根據(jù)圖1的陷落柱突水概念模型，建立了如圖2所示的陷落柱突水計(jì)算模型。設(shè)模型長40m，高22m，其中底部含水層長40m，高6m，中部陷落柱平均高14m，平均寬5m，右上側(cè)巷道平均長18m，高2m。底部含水層的邊界為進(jìn)水流量邊界，單寬流速為：uesdln=1×10-3m3/(s·m)，相當(dāng)于2.6MPa的承壓水壓。上部巷道出口的邊界為流體壓力邊界，等于大氣壓力。含水層中流體的黏滯系數(shù)為1×10-3Pa·s，先假設(shè)含水層和陷落柱的滲透率相同，為1×10-12m2，分析Brinkman區(qū)域和Darcy區(qū)域同樣滲透系數(shù)但不同滲流方程情況下Brinkman方程引起流速的增大幅度，然后把陷落柱Brinkman區(qū)域的滲透率提高10倍，分別研究當(dāng)含水層補(bǔ)給水量不變情況下高滲透率的陷落柱引起的水壓卸壓幅度和當(dāng)保持含水層水壓不變情況下的流速急劇增大的幅度，進(jìn)一步揭示陷落柱非Darcy滲流對突水過程的作用機(jī)制。4.2.2brifman區(qū)域滲透率計(jì)算結(jié)果如圖3～6所示。圖3,4分別顯示了流體從Darcy區(qū)域到Brinkman區(qū)域再到巷道中Navier-Stoke流動的水流速度和壓力的變化。圖5,6分別為圖4中沿折線A1—A2—A3—A4方向切線流速和壓力的分布曲線(圖中：工況1表示Darcy區(qū)域和Brinkman區(qū)域滲透率一致；工況2表示Brinkman區(qū)域滲透率提高10倍，但邊界流量和工況1一致；工況3表示Brinkman區(qū)域滲透率提高10倍，且保持邊界恒定高壓力)。陷落柱突水滲流模擬結(jié)果可見，速度變化(表面和箭頭)和壓力分布(表面)沿3個(gè)流域是連續(xù)變化的，流速梯度的大幅度變化主要集中在Brinkman區(qū)域，壓力梯度變化則主要在Darcy區(qū)域，表明水流在Darcy區(qū)域內(nèi)的滲流阻力最大。圖5,6表明，當(dāng)Darcy區(qū)域和Brinkman區(qū)域采用相同滲透率時(shí)，雖然在Darcy-Brinkman交界面和BrinkmanNavier-Stokes交界面上流速急劇增大，但過渡比較連續(xù)。在Darcy區(qū)域，流速為2.0×10-5m/s，接近到Brinkman區(qū)域，流速逐漸增大，在邊界面流速為7.0×10-5m/s，即使2個(gè)區(qū)域的滲透系數(shù)完全相等，由于Brinkman考慮了流體動能的耗散作用，在Brinkman區(qū)域流速也急劇增大，進(jìn)入到Navier-Stokes區(qū)域的邊界面，流速達(dá)到2.3×10-4m/s，增大了10倍，進(jìn)入到Navier-Stokes區(qū)域后，流速進(jìn)一步增大，達(dá)到3.5×10-4m/s，增大了17倍。而基于線性滲流理論，整個(gè)區(qū)域只用Darcy方程計(jì)算時(shí)，流速不會有這么急劇的階躍變化，無法描述突水時(shí)流量的非線性變化過程。而當(dāng)將Brinkman區(qū)域滲透率提高10倍時(shí)，該條件下若保持補(bǔ)給水量不變，將造成邊界水壓力從2.6MPa降低到0.7MPa，降低了2.7倍，表明當(dāng)含水層補(bǔ)給水量一定時(shí)，Brinkman區(qū)域滲透性的提高會起到較大幅度的卸壓作用，從而降低了壓力梯度，Brinkman區(qū)域滲透性的提高對水流速度變化沒有影響，進(jìn)一步研究表明，揭穿陷落柱時(shí)，若含水層補(bǔ)給水量有限，即使陷落柱滲透性很高，也不會發(fā)生突水災(zāi)害。當(dāng)含水層高水壓保持不變時(shí)，同樣的水壓梯度(水流進(jìn)入巷道中降為大氣壓)下，Brinkman區(qū)域滲透率的增大引起進(jìn)入Brinkman和Navier-Stokes區(qū)域水流流速更大幅度地增高。在Darcy區(qū)域，流速為1.4×10-5m/s，接近到Brinkman區(qū)域，流速急劇增大，在邊界面流速為5.1×10-4m/s，流速增大幅度出現(xiàn)“階越”式變化，進(jìn)入到Navier-Stokes區(qū)域的邊界面，流速達(dá)到1.2×10-3m/s，增大了85倍，進(jìn)入到Navier-Stokes區(qū)域后，流速進(jìn)一步增大，達(dá)到1.6×10-3m/s，增大了114倍。由此可見，陷落柱作為含水層滲流和巷道突水自由流動的過渡區(qū)域，其滲透性變化對于突水壓力和流速演變作用十分顯著，這進(jìn)一步表明，現(xiàn)場實(shí)際的突水現(xiàn)象就是在含水層充足的補(bǔ)給水量和保持恒定的高水壓作用下，當(dāng)陷落柱或?qū)扑閹?lián)通了含水層和巷道情況下形成的。5巷道滲流場模擬結(jié)果分析本文提出突水過程中含水層、破碎導(dǎo)水帶和巷道3種水流運(yùn)移耦合模型，采用Brinkman方程計(jì)算破碎帶水流場，比較符合破碎巖體介質(zhì)的流場特性，同時(shí)把含水層Darcy滲流和巷道Navier-Stokes流動有機(jī)聯(lián)系在統(tǒng)一流動場中，得到整個(gè)水流流場分布，更符合突水實(shí)際。采用非線性耦合滲流