浦 海，郭世儒，劉德俊，許軍策，王 健

（1.中國礦業(yè)大學(xué) 深部巖土力學(xué)與地下工程國家重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116；2. 新疆工程學(xué)院 采礦工程與地質(zhì)學(xué)院，新疆 烏魯木齊 830091；3.中國礦業(yè)大學(xué) 力學(xué)與土木工程學(xué)院，江蘇 徐州 221116）

0 引 言

隨著我國煤炭開發(fā)基地的戰(zhàn)略性西移，西部礦區(qū)煤層開采時(shí)因其獨(dú)特的地質(zhì)條件和高強(qiáng)度開采方式涌現(xiàn)出許多工程問題，其中在采掘工作面頻發(fā)的突水潰砂災(zāi)害已成為影響煤礦安全生產(chǎn)的關(guān)鍵問題之一。 高強(qiáng)度開采過程中的擾動(dòng)促使上覆含水層底部產(chǎn)生裂隙，在裂隙發(fā)育形成導(dǎo)水通道后進(jìn)一步貫通厚松散砂層使水砂混合物涌入工作面是突水潰砂災(zāi)害形成的主要方式之一。 在形成導(dǎo)通裂隙后，水砂運(yùn)移的方式將決定礦井的危害程度，迫切需要研究水砂運(yùn)移規(guī)律，尤其是在形成裂隙通道時(shí)，水砂在其內(nèi)部的流動(dòng)規(guī)律及對(duì)通道的擴(kuò)張影響［1-3］。

近年來，眾多學(xué)者針對(duì)突水潰砂問題在理論分析、相似試驗(yàn)?zāi)M、基礎(chǔ)試驗(yàn)驗(yàn)證、數(shù)值計(jì)算模擬等方面做出大量的研究成果，為深入了解突水潰砂發(fā)生機(jī)理與災(zāi)害防控提供了有意義的參考［4-12］。 基于工程分析角度，張桂民等［13］分析探討了煤礦鉆孔誘發(fā)的突水潰砂機(jī)理，提出了防治和治理突水潰砂的措施。 許延春等［14］研究了工作面突水潰砂機(jī)制，并提出相應(yīng)的預(yù)防措施。 張勃陽等［15］采用流體力學(xué)方法，分析煤礦潰砂事故，提出了鉆孔潰砂災(zāi)害的通用計(jì)算方法。 基于試驗(yàn)研究角度，張凱等［16］根據(jù)實(shí)際工程情況進(jìn)行了涌砂判據(jù)的試驗(yàn)研究，通過滲流試驗(yàn)研究了砂巖的水力特性，得出涌砂發(fā)生的臨界條件。 李宏杰等［17］進(jìn)行了突水防砂的實(shí)例研究，根據(jù)試驗(yàn)結(jié)果提出了分層開采、降低上下分層開采厚度等優(yōu)化設(shè)計(jì)和技術(shù)措施。 袁奇［18］分析了陜北礦區(qū)潰砂主要來源，并設(shè)計(jì)整套試驗(yàn)裝置分析突水潰砂現(xiàn)象，攝錄了不同砂粒粒徑、砂厚和水頭高度情況下潰砂剖面特征，直觀展示了突水潰砂過程。 楊鑫等［19］采用自主研發(fā)的水砂兩相高速滲流試驗(yàn)設(shè)備，對(duì)西部榆橫礦區(qū)風(fēng)積砂含水層的起動(dòng)、潰砂規(guī)律及水砂兩相高速滲流進(jìn)行室內(nèi)試驗(yàn)研究。 提出西部典型礦區(qū)風(fēng)積砂含水層突水潰砂災(zāi)害的臨界判據(jù)，并分析潰砂過程中的能量傳遞機(jī)制。 楊偉峰等［20］設(shè)計(jì)了水砂混合涌流裝置，模擬了覆巖裂隙通道水砂混合突水啟動(dòng)過程及遷移過程，分析了裂縫通道中水砂混合流的突水特征。 高炳倫［21］針對(duì)工程中裂隙潰砂問題研制了裂隙潰砂試驗(yàn)系統(tǒng)，實(shí)現(xiàn)了單一平直裂隙條件下裂隙潰砂過程的試驗(yàn)?zāi)M。 基于模擬研究角度，王迎超等［22］采用DEM 數(shù)值模擬方法建立了固-水兩相流動(dòng)的三維模型，直觀地展現(xiàn)了突水潰砂時(shí)水砂運(yùn)移過程。 梁艷坤等［23-24］采用分形理論分析了垮落帶破碎巖體的塊度分布特征，并利用離散元方法（DEM） 根據(jù)破碎巖體的塊度分布生成不規(guī)則碎石，建立了垮落帶破碎巖體潰砂數(shù)值模型。

目前針對(duì)突水潰砂問題的研究并不透徹，主要原因是地下空間為封閉不可觀測系統(tǒng)，地層水砂分布、承壓含水層壓力、巖層裂隙網(wǎng)絡(luò)分布等因素難以觀測。 加之水砂性質(zhì)的多樣，進(jìn)行相似模擬試驗(yàn)時(shí)選取相似材料過程繁瑣，通過試驗(yàn)手段對(duì)突水潰砂危害性評(píng)估有很大的難度。 數(shù)值模擬因其可重復(fù)性和直觀性等特點(diǎn)，在研究突水潰砂災(zāi)害時(shí)有明顯的優(yōu)勢。

突水潰砂數(shù)值模擬中，重點(diǎn)在于針對(duì)砂粒及流體選取合適的計(jì)算模型，能夠?qū)烧哌M(jìn)行耦合以得出符合實(shí)際的運(yùn)動(dòng)過程。 傳統(tǒng)的研究方法中將砂看成另一種形式的流相，水砂混合流動(dòng)看成兩相流問題，或者將砂簡化為質(zhì)點(diǎn)，忽略顆粒之間的碰撞及翻轉(zhuǎn)，在分析突水潰砂下水砂在裂隙網(wǎng)絡(luò)中的運(yùn)移規(guī)律時(shí)達(dá)不到真實(shí)效果。 為使模擬更符合實(shí)際問題，筆者采用格子Boltzmann 方法（Lattice Boltzmann Method，LBM） 模 擬 液 相 流 動(dòng)， 用 離 散 單 元 法（Discrete Element Method，DEM）構(gòu)建球形砂粒，并采用浸入移動(dòng)邊界法（Immersed Moving Boundary Method， IMB）方法實(shí)現(xiàn)流體和砂粒間的耦合，由此模擬砂粒在流場中的流動(dòng)情況，有效地分析突水潰砂下水砂兩相流動(dòng)運(yùn)移規(guī)律。

1 突水潰砂分析模型

西部礦區(qū)煤層開采時(shí)突水潰砂災(zāi)害的演變過程如圖1 所示。 西部礦區(qū)煤層開采時(shí)劇烈的擾動(dòng)使上層覆巖形成大量裂隙網(wǎng)絡(luò)，當(dāng)擾動(dòng)加劇時(shí)將促使裂隙網(wǎng)絡(luò)貫通巷道和上層含砂水層，水砂混合流將通過裂隙大量涌入生產(chǎn)工作面，形成突水潰砂災(zāi)害，帶來巨大的經(jīng)濟(jì)損失和人員傷亡。 由于西部礦區(qū)地表存在風(fēng)積砂，若擾動(dòng)加劇上層含砂水層形成與地表的貫通裂隙，將有源源不斷的砂供給，進(jìn)一步加劇突水潰砂的災(zāi)害性。 為明確突水潰砂的成災(zāi)原因及危害程度，重點(diǎn)研究圖1 中紅圈位置水砂混合流在突發(fā)裂隙開口情況下的運(yùn)移情況。 將實(shí)際突水潰砂時(shí)含砂水層中水砂進(jìn)入破裂巖體這一情況簡化為恒定壓力下的單裂隙開口通道水砂運(yùn)移數(shù)值計(jì)算模型如圖2 所示。

圖1 西部礦區(qū)突水潰砂示意Fig.1 Schematic of water and sand inrush in western mining area

圖2 突水潰砂簡化計(jì)算模型Fig.2 Simplified analysis model of water and send inrush

矩形通道中寬度為d的喉道口表征突水潰砂發(fā)生時(shí)突然增大的裂隙開口，邊界壓力P表征地下含砂水層中的水力壓力情況。 水分布在開口上方的矩形區(qū)域，密度為ρw、動(dòng)力黏度為υ。 砂粒用半徑為r、密度為ρs、厚度為h、數(shù)量為N的球形顆粒群表征。 由于砂粒分布區(qū)域固定，且分析時(shí)以自然堆積方式排布，砂層厚度和砂粒數(shù)量呈正比，僅取1 項(xiàng)作為變量即可。 以P、d、h作為變量，從數(shù)值模擬角度分析水砂混合流在單裂隙下的運(yùn)移規(guī)律。

2 LBM-DEM 耦合方法及模擬方案

LBM-DEM 耦合方法在模擬中對(duì)顆粒進(jìn)行直接建模，顆粒擁有物理邊界和確切的體積，流體假定為連續(xù)相并構(gòu)建小于流場中最小渦長尺寸的網(wǎng)格進(jìn)行直接數(shù)值模擬。 隨著越來越多的學(xué)者不斷完善基于LBM-DEM 方法下的液固耦合理論，此方法能夠很好的分析液固耦合時(shí)顆粒流動(dòng)和流場分布情況，并開始廣泛應(yīng)用于生物、醫(yī)療等多領(lǐng)域中［25］。

2.1 格子Boltzmann 方法（LBM）

格子Boltzmann 方法是應(yīng)用非連續(xù)介質(zhì)思想求解流體力學(xué)問題的方法。 其基本思想是把流體看成由許多只有質(zhì)量沒有體積的微小粒子組成，在時(shí)間、空間上離散，按照符合特定基本物理規(guī)律的演化規(guī)則碰撞，并沿網(wǎng)格線在節(jié)點(diǎn)之間遷移，利用分布函數(shù)的統(tǒng)計(jì)平均得到系統(tǒng)的宏觀物理參數(shù)［26］。

借助于離散的單松弛時(shí)間格子Boltzmann -BGK 方程描述流體相互作用，即

其中；x為粒子空間位置矢量；ci為粒子速度矢量；t為時(shí)間；δt為時(shí)間步長；fi x，t( ) 為粒子分布函數(shù)；Fi x，t( ) 為離散的外力項(xiàng)；Ωi為離散的碰撞算子，可表達(dá)式為

其中：τc為2 次碰撞間隔時(shí)間；u為流體宏觀速度矢量；fi，equ，ρ( ) 為離散后的平衡態(tài)速度分布函數(shù)，表達(dá)式為

其中：ωi為權(quán)重系數(shù)；cs為格子聲速，表達(dá)式為

式中：R為氣體常數(shù)；T為熱力學(xué)溫度；D為空間維數(shù)。

考慮到實(shí)際應(yīng)用，選取單松弛時(shí)間的D3Q19（3是空間維數(shù)，19 是離散速度數(shù)）三維模型，如圖3 所示。 并根據(jù)三維模型各點(diǎn)坐標(biāo)表征ci，i＝0，1，2，…，18，指可能遷移方向數(shù)。

圖3 D3Q19 模型Fig.3 D3Q19 model

若知道粒子分布函數(shù)，可計(jì)算模型的宏觀密度、速度，其表達(dá)式為

2.2 離散單元法（DEM）

離散單元法被廣泛應(yīng)用于巖石力學(xué)和顆粒流等離散系統(tǒng)模擬中，是基于非連續(xù)介質(zhì)力學(xué)的數(shù)值方法［27］。 計(jì)算各個(gè)時(shí)間步上顆粒間所受的接觸力及顆粒運(yùn)動(dòng)，分析時(shí)將運(yùn)動(dòng)分為平動(dòng)和轉(zhuǎn)動(dòng)，平動(dòng)在全局坐標(biāo)系下計(jì)算，由牛頓方程給出，即

轉(zhuǎn)動(dòng)在局部坐標(biāo)系下計(jì)算，需要計(jì)算出繞顆粒質(zhì)心的合力矩，顆粒轉(zhuǎn)動(dòng)的控制方程為

式中：Ip為主慣性矩；ω為角速度；rk為顆粒上第k個(gè)接觸處的矢徑；Md為阻尼力矩。

2.3 浸入移動(dòng)邊界法（IMB）

浸入移動(dòng)邊界法是LBM 和DEM 的耦合方法，其具備碰撞算子的局部性和遷移操作的簡單性，并在此基礎(chǔ)上克服了傳統(tǒng)耦合方法中的動(dòng)量不連續(xù)性問題，具有足夠代表性的網(wǎng)格非一致邊界［28-30］。

所有離散的格子被分為純流體格子、邊界格子和純固體格子3 類。 引入格子固含率εs定義每個(gè)格子內(nèi)流固相互作用的影響，并通過格子內(nèi)顆粒占據(jù)面積與格子占據(jù)總面積之比定量表征格子內(nèi)流場和顆粒分布情況。 流體格子的固含率εs＝0，固體格子固含率εs＝1，邊界格子的固含率εs為0～1，計(jì)算時(shí)用砂粒所占面積除以格子總面積，如圖4 所示。

圖4 格子分類示意Fig.4 Schematic of grid classification

為表征顆粒對(duì)流體的影響，在忽略外力項(xiàng)時(shí)將Boltzmann - BGK 方程添加源項(xiàng)變?yōu)?/p>

β（εs，τ） 是與格子固含率有關(guān)的附加碰撞項(xiàng)的加權(quán)函數(shù)。 表示為格子固含率εs和松弛因子τ的函數(shù)為

為考慮流體對(duì)顆粒的影響，通過加入所有被當(dāng)前顆粒覆蓋的格子上對(duì)顆粒的力和力矩表示周圍流體作用于該顆粒的力和力矩，該力即顆粒受到的曳力，曳力F和力矩T的表達(dá)式分別為

其中，n為該顆粒覆蓋邊界格子的數(shù)目；m為離散速度數(shù)；xj為第j個(gè)邊界格子的位置；xc為顆粒中心位置。 根據(jù)受力和力矩，顆粒對(duì)位置進(jìn)行更新，其顆粒平動(dòng)速度和轉(zhuǎn)動(dòng)角速度的更新式為

式中：ρL為流體密度；ρp為顆粒密度；mL為流體質(zhì)量；g為重力加速度；Δt為計(jì)算時(shí)間步長。 最終通過顆粒和流體作用方程的不斷迭代計(jì)算，實(shí)現(xiàn)水砂混合流的液固耦合計(jì)算。

2.4 計(jì)算流程

基于LBM-DEM 實(shí)現(xiàn)了水砂混合流的液固耦合，大量學(xué)者已通過對(duì)比單顆粒沉降及兩顆粒DKT（追擊、接觸、翻轉(zhuǎn)）試驗(yàn)過程驗(yàn)證LBM-DEM 耦合方法的合理性，可用于進(jìn)行水砂混合流的液固耦合計(jì)算［31］。 其計(jì)算流程為在初始t0時(shí)刻輸入流場和顆粒信息，判斷各格點(diǎn)上的格子固含率。 通過LBM計(jì)算下一時(shí)刻t1時(shí)流場信息，并根據(jù)DEM 循環(huán)計(jì)算顆粒上的碰撞，將其信息更新為t1時(shí)刻。 依次對(duì)時(shí)間步循環(huán)，直至判定精度達(dá)到要求結(jié)束計(jì)算，實(shí)現(xiàn)水砂混合流的液固耦合計(jì)算過程如圖5 所示。

2.5 模型構(gòu)建及模擬方案

考慮到實(shí)際突水潰砂問題，構(gòu)建單裂隙通道三維數(shù)值計(jì)算模型如圖6 所示。

圖5 LBM-DEM 耦合計(jì)算過程Fig.5 Calculation process of LBM-DEM

圖6 單裂隙通道三維數(shù)值計(jì)算模型Fig.6 Three-dimensional numerical model of single fracture passage

以裂隙開口檔板最左側(cè)為O 點(diǎn)、建立YOZ 坐標(biāo)系，三維空間尺寸為0.015 m×0.400 m×0.600 m（X×Y×Z），水流區(qū)域和初始砂粒區(qū)域均位于裂隙開口上方，上下兩面作用恒定壓力驅(qū)使水流流動(dòng)，砂粒計(jì)算前先進(jìn)行自由堆積，再通過坐標(biāo)導(dǎo)入計(jì)算模型，使其均勻分布在開口上方。 模型參數(shù)設(shè)置：ρw＝1 000 kg/m3；υ＝0.001 Pa·s；ρs＝2 500 kg/m3；r＝0.002 m；E＝4.9 GPa。

為探究突水潰砂成災(zāi)的主導(dǎo)因素，改變不同邊界壓力P、裂隙開口寬度d和砂層厚度h分析水砂混合流在模型中的流動(dòng)情況，統(tǒng)計(jì)裂隙開口處截面流量及潰砂速率的變化，方案見表1。

表1 模擬方案Table 1 Simulation schemes

3 計(jì)算結(jié)果分析

3.1 不同邊界壓力下水砂運(yùn)移規(guī)律

模擬不同邊界壓力下水砂混合流的運(yùn)移情況如圖7 所示，隨著邊界壓力的增加，潰砂量和裂隙開口處最大流速明顯增加。 為定量研究潰砂增加量及流場變化，對(duì)裂隙開口處截面的流量Sw和單位時(shí)間潰砂速率vs與時(shí)間t關(guān)系進(jìn)行分析（圖8、圖9）。

隨著邊界壓力的增加，最終的穩(wěn)定流量逐漸增加（圖8）。 邊界壓力較小時(shí)流量先增加后保持恒定，而高壓力時(shí)在達(dá)到早期穩(wěn)定區(qū)持續(xù)一段時(shí)間后流量會(huì)再次升高，以P-5 最為明顯。 對(duì)比圖9 中潰砂速率變化可發(fā)現(xiàn)同樣的規(guī)律，壓力較小時(shí)潰砂速率先增至峰值后保持恒定，而隨著壓力的增加，潰砂速率從初始到峰值的增幅變大，但在峰值平臺(tái)區(qū)持續(xù)的時(shí)間變短，之后出現(xiàn)明顯降低，且隨著壓力增加降幅逐漸增大。 考慮到流量二次增加的時(shí)間點(diǎn)和潰砂速率從峰值開始降低的時(shí)間基本一致。 以P-5為例，分析了5 個(gè)等距時(shí)間步上截面流量、潰砂速率及潰砂形貌（圖10），并通過紅色虛線標(biāo)記不同壓力時(shí)顆粒中的主要受力區(qū)。

圖7 不同邊界壓力下水砂運(yùn)移Fig.7 Migration of sediment at different boundary pressures

圖8 不同邊界壓力下截面流量變化Fig.8 Flow changes of cross section under different boundary pressures

圖9 不同邊界壓力下潰砂速率變化Fig.9 Changes of sand breaking rate underdifferent boundary pressures

從圖10 所知，在①—③步中截面流量和潰砂速率基本不發(fā)生變化，此時(shí)在水流驅(qū)動(dòng)作用下砂粒將在裂隙出口處形成臨時(shí)的拱狀受力結(jié)構(gòu)，顆粒間相互作用使砂粒更加密集，對(duì)流場有一定的阻礙作用。隨著砂粒的不斷流出，出口處的拱狀結(jié)構(gòu)逐漸變化，其中心處不斷凹陷。 至第④步時(shí)拱狀結(jié)構(gòu)變?yōu)轭愃啤癕”狀結(jié)構(gòu)，仍對(duì)水流及流砂有阻礙作用。 直至第⑤步時(shí)，中心砂粒完全流出喉道口，拱狀結(jié)構(gòu)完全消失后不再對(duì)流場起到阻礙作用，此時(shí)流量會(huì)快速上升，而由于砂粒已大量流出其潰砂速率會(huì)快速下降。

圖10 P-5 截面流量和潰砂速率對(duì)比Fig.10 Comparison of cross section flow rate and sediment burst rate of Orifice P-5

為表征不同邊界壓力對(duì)水砂流動(dòng)的影響程度，統(tǒng)計(jì)單位時(shí)間潰砂速率最大值及初始階段受砂粒阻礙作用的平均流量隨壓力的變化規(guī)律如圖11 所示。其中隨著邊界壓力的增加，最大潰砂速率和平均流量均單調(diào)遞增。 計(jì)算最大潰砂速率及初始平均孔口流量隨壓力的區(qū)間增幅，并通過區(qū)間增幅除以區(qū)間范圍得出單位壓力增幅以表征不同壓力數(shù)值下的影響情況，計(jì)算結(jié)果見表2。

圖11 截面流量和潰砂速率隨邊界壓力變化Fig.11 Variations of cross section flow rate and sediment burst rate with boundary pressures

表2 截面流量和潰砂速率隨邊界壓力的增幅Table 2 Increase of cross section flow rate and sediment burst rate with boundary pressure

從表3 發(fā)現(xiàn)，低壓力時(shí)隨壓力的增加最大潰砂速率和平均流量的增幅均大于高壓力時(shí)。 說明僅考慮邊界壓力情況下，隨著邊界壓力的增加，同等時(shí)間下突水潰砂災(zāi)害程度會(huì)逐漸加劇，但低壓力下由于壓力增加引發(fā)的災(zāi)害加劇程度較大，此時(shí)應(yīng)注意由壓力增加引起的劇烈波動(dòng)。 由于模型中設(shè)置恒定壓差表征地下承壓含水層作用，未考慮實(shí)際流動(dòng)中壓力降低的作用，結(jié)論僅能在一定程度上適用于工程實(shí)際問題，后續(xù)將針對(duì)復(fù)雜模型繼續(xù)研究。

針對(duì)干砂單裂隙開口通道下流出孔口時(shí)單位質(zhì)量流量計(jì)算，1961 年BEVERLOO 等［32］最早給出單位質(zhì)量流量W的計(jì)算式為

式中：C為與顆粒間摩擦性質(zhì)相關(guān)的無因次量；ρB為顆粒堆積密度，通過顆粒密度ρs與空隙度ε計(jì)算，表示為ρB＝ρs1-ε( ) ；D0為孔口直徑；dp為顆粒直徑；k為無因次系數(shù)。

但式（1） 未考慮流體與顆粒并存情況，BULSARA 等［33］引入壓力梯度概念得出含流體時(shí)顆粒流出孔口時(shí)單位質(zhì)量流量計(jì)算式為

其中：P1和P2分別為孔板下方及上方壓力。 郭帥等［34］通過試驗(yàn)分析驗(yàn)證了其合理性，模擬時(shí)設(shè)置不同恒定壓差以驅(qū)動(dòng)水流，參照式（2）給出最大潰砂量隨壓差的計(jì)算式為

其中：ΔP為模型恒定壓差；W0為干砂在孔口的單位質(zhì)量流量，能有效表征壓差為零時(shí)干砂受自重下落的情況。 圖12 給出式（3）擬合曲線與數(shù)值模擬數(shù)據(jù)點(diǎn)對(duì)比情況，曲線擬合情況較好。

圖12 單位質(zhì)量流量數(shù)據(jù)點(diǎn)與擬合曲線對(duì)比Fig.12 Comparison of data points per unit mass flow and fitting curve

3.2 不同開口寬度下水砂運(yùn)移規(guī)律

模擬不同開口寬度下的水砂混合流運(yùn)移情況表明：隨開口寬度的增加，潰砂量及開口處最大流速明顯增加，且流速高速區(qū)逐漸貫穿整個(gè)流場。 為明確潰砂增加量及流場變化情況，圖13、圖14 中結(jié)合裂隙開口處相同截面的流量Sw和單位時(shí)間潰砂速率vs進(jìn)行分析。

圖13 不同開口寬度下截面流量變化Fig.13 Variation of flow rate of orifice section with different opening widths

隨著開口寬度的增加，最終的穩(wěn)定流量逐漸增加，如圖13 所示。 開口較小時(shí)出口截面流量在達(dá)到初始流量峰值后保持恒定一段時(shí)間后再逐漸增加，開口較大時(shí)，達(dá)到初始流量峰值后將出現(xiàn)明顯的下降階段，之后又由于顆粒大量流出開口使流量逐漸增至最終穩(wěn)定流量，以K-4、K-5 最為明顯。

隨著開口寬度的增加，潰砂速率從初始到峰值的增幅變大，且在峰值平臺(tái)區(qū)持續(xù)的時(shí)間變短，如圖14 所示。 開口寬度較小時(shí)潰砂速率增加到峰值將保持恒定一段時(shí)間，后由于砂粒逐漸流出裂隙出現(xiàn)緩慢降低。 開口寬度較大時(shí)，潰砂速率達(dá)到峰值保持一定時(shí)間后由于砂粒流出開口急速降低。 在開口進(jìn)一步擴(kuò)大后，甚至在峰值持續(xù)階段出現(xiàn)一段小幅降低。

圖14 不同開口寬度下潰砂速率變化Fig.14 Change of sediment discharge rate under different opening widths

圖15 K-4 截面流量和潰砂速率對(duì)比Fig.15 Comparison of cross section flow rate and sediment burst rate of K-4 orifice

因此，以K-4 為例，結(jié)合流量、潰砂速率、小幅降低時(shí)刻的水砂分布情況進(jìn)行分析，并用紅色虛線標(biāo)示出大量砂粒受擠壓聚集在開口處的現(xiàn)象，如圖15 所示。 分析可知砂粒受擠壓在開口處形成密實(shí)堆積結(jié)構(gòu)影響水砂運(yùn)移，造成潰砂速率和截面流量小幅下降，但由于結(jié)構(gòu)中的砂粒具有流動(dòng)性，形成的密實(shí)結(jié)構(gòu)會(huì)隨著砂粒流出開口逐漸失去其阻礙效果，直至砂粒中心平面降至孔口下方時(shí)該結(jié)構(gòu)消失。

圖16 中選取潰砂速率峰值時(shí)不同孔口下水砂分布形態(tài)進(jìn)行分析，如圖16 所示，以紅色虛線畫出不同開口上方的密實(shí)結(jié)構(gòu)區(qū)域，隨著開口寬度的增加，在孔口處形成的密實(shí)結(jié)構(gòu)體積不斷增加，對(duì)流場造成的影響程度將不斷增大。 開口較小時(shí)，由于密實(shí)結(jié)構(gòu)范圍小，其受擠壓砂粒相對(duì)總體砂粒占比較小，影響效果不明顯。 開口較大時(shí)密實(shí)區(qū)范圍較大，影響砂粒數(shù)量占比較多，將對(duì)整體水砂運(yùn)移產(chǎn)生巨大影響。

圖16 不同開口寬度下水砂運(yùn)移Fig.16 Migration diagram of sediment with different opening widths

為表征不同開口寬度對(duì)水砂流動(dòng)的影響程度，圖17 中給出單位時(shí)間潰砂速率最大值及初始階段受砂粒阻礙作用的平均流量隨開口寬度的變化規(guī)律。

圖17 截面流量和潰砂速率隨開口寬度的變化Fig.17 Variation of cross section flow rate and sand breaking rate of orifice vary with opening width

隨著邊界壓力的增加，最大潰砂速率和初始階段平均流量均單調(diào)遞增（圖17）。 為量化開口寬度的影響，表3 給出最大潰砂速率和初始平均流量隨開口寬度的區(qū)間增幅。

通過表3 發(fā)現(xiàn)，開口寬度較小時(shí)隨寬度增加對(duì)最大潰砂速率和初始平均流量的影響大于開口寬度較大時(shí)。 說明固定的砂含量下，隨著裂隙開口寬度的增加，突水潰砂災(zāi)害程度的增幅逐漸減小，在寬度大于一定值后，由于砂含量的限制災(zāi)害逐漸趨于恒定。 突增現(xiàn)象最容易出現(xiàn)在細(xì)小裂隙開口擴(kuò)張的過程中。

表3 截面流量和潰砂速率隨開口寬度的增幅Table 3 Variation of cross section flow rate and sand breaking rate increase with width of opening

3.3 不同砂層厚度下水砂運(yùn)移規(guī)律

密實(shí)填充下砂層厚度與砂粒數(shù)目呈正比，不同砂層厚度可表征不同砂粒數(shù)目對(duì)水砂運(yùn)移的影響。模擬不同砂層厚度下水砂運(yùn)移結(jié)果表明：隨著砂層厚度的增加，砂粒數(shù)目增加，最終的潰砂量明顯增加，且開口兩側(cè)堆積的未流出砂粒也逐漸增加，但流速未產(chǎn)生明顯變化。 為明確潰砂增加量及流場變化情況，圖18、圖19 結(jié)合開口處的截面流量Sw和單位時(shí)間潰砂速率vs進(jìn)行分析。

圖18 不同砂層厚度下截面流量變化Fig.18 Variation of cross section flow of orifice with sand thickness

圖19 不同砂層厚度下潰砂速率變化Fig.19 Changes of sediment discharge rate under different sand thickness

從圖18 所知，最終穩(wěn)定流量不隨砂層厚度的增加發(fā)生變化。 而隨著砂層厚度的增加，截面流量在達(dá)到初始峰值后均出現(xiàn)一定的下降，用對(duì)應(yīng)曲線顏色的虛線記錄最終時(shí)間，①表示初始流量峰值點(diǎn)，②③④⑤分別對(duì)應(yīng)N-2、N-3、N-4、N-5 時(shí)最終下落結(jié)束的時(shí)間點(diǎn)，發(fā)現(xiàn)下降持續(xù)的時(shí)間隨著砂層厚度的增加逐漸增長。 說明隨著砂層厚度的增加，下落過程中對(duì)水流的阻礙作用越明顯。

如圖19 所示，隨著砂層厚度的增加，潰砂速率峰值點(diǎn)位置不斷延后。 在砂層厚度較大時(shí)，均呈現(xiàn)先急速增至初始高速率后持續(xù)降低一段時(shí)間，再升至峰值速率后快速下降的趨勢。 同樣用①表示初始速率峰值點(diǎn)，②③④⑤分別對(duì)應(yīng)N-2、N-3、N-4、N-5 時(shí)速率持續(xù)下降結(jié)束的時(shí)間點(diǎn)，發(fā)現(xiàn)隨著砂層厚度的增加，從初始高速率降至最小速率的持續(xù)時(shí)間不斷增長，說明砂層厚度的增加對(duì)砂層整體潰突阻礙作用逐漸增加。 為分析砂層厚度較大時(shí)對(duì)水砂運(yùn)移的阻礙形式，圖20 以N-5 為例，結(jié)合流量、潰砂速率及不同時(shí)刻水砂分布情況進(jìn)行分析。

圖20 N-5 截面流量和潰砂速率對(duì)比Fig.20 Comparison of cross section flow rate and sediment burst rate of Orifice N-5

從圖20 可知，對(duì)水流及整體砂粒流速起阻礙作用的原因主要為砂粒在裂隙開口處聚集形成密實(shí)區(qū)域，對(duì)比各階段可以發(fā)現(xiàn)①②③④時(shí)刻大量砂粒受力后聚集在裂隙開口，使得水砂流動(dòng)減緩，在⑤時(shí)刻砂粒形成的密實(shí)結(jié)構(gòu)崩塌，對(duì)流體的阻礙消失，截面流量逐漸開始上升。 而此時(shí)由于密實(shí)結(jié)構(gòu)崩塌后大量砂粒傾瀉將造成此時(shí)的潰砂速率達(dá)到峰值。 為明確不同砂層厚度對(duì)水砂運(yùn)移的影響，圖21 分析不同砂層高度下單位時(shí)間潰砂速率最大值及初始階段受砂粒阻礙作用階段平均流量的變化規(guī)律。

圖21 截面流量和潰砂速率隨砂層厚度的變化Fig.21 Variation of cross section flow rate of orifice and rate of sediment burst change with thickness of sand layer

隨著砂層厚度的增加，最大潰砂速率和平均流量均單調(diào)遞減（圖21），在砂層厚度較小時(shí)出現(xiàn)均快速降低的現(xiàn)象。 當(dāng)砂層厚度較大時(shí)，隨砂層厚度增加對(duì)最大潰砂速率和平均流量的影響較小。 為定量研究砂層厚度對(duì)水砂運(yùn)移的影響程度，表4 計(jì)算給出不同砂層厚度下最大潰砂速率和平均流量的區(qū)間降幅。

表4 截面流量和潰砂速率隨砂層高度的降幅Table 4 Decrease of cross section flow rate and sediment burst rate with thickness of sand layer

由表5 知，在砂層厚度較小時(shí)隨砂層厚度增加最大潰砂速率和平均流量的降幅遠(yuǎn)大于砂層厚度較大時(shí)，在N-2 時(shí)的降幅基本為N-5 時(shí)的10 倍。 說明在邊界壓力、砂層范圍及裂隙開口恒定的情況下，隨著砂層厚度的增加，包含的砂粒數(shù)量增多，同等時(shí)間下由于砂粒阻礙作用增強(qiáng)其突水潰砂災(zāi)害程度會(huì)逐漸減小。 當(dāng)砂層厚度較小時(shí)，水驅(qū)動(dòng)砂快速流出裂隙孔口，災(zāi)害最為嚴(yán)重。 而砂層越厚砂粒越容易在裂隙口擁堵阻礙砂粒及流體運(yùn)動(dòng)，此時(shí)在砂粒密實(shí)結(jié)構(gòu)失穩(wěn)造成砂粒傾瀉時(shí)發(fā)生最嚴(yán)重的災(zāi)害。 分析僅考慮同等時(shí)間下的災(zāi)害影響，未研究不同砂層厚度砂粒全部流出情況，此情況明顯為砂層厚度越大，可流出砂粒增加造成災(zāi)害愈加嚴(yán)重。

4 結(jié)論及展望

1）隨邊界壓力增加，截面流量及最大潰砂速率均單調(diào)增加，同等時(shí)間下突水潰砂災(zāi)害程度會(huì)逐漸加劇。 由于高壓力下砂粒更加容易在喉道口堆積成臨時(shí)性密實(shí)結(jié)構(gòu)阻礙水砂運(yùn)移，升高單位壓力時(shí)引起的最大潰砂量及平均截面流量的增幅逐漸降低。在Bulsara 式基礎(chǔ)上考慮流體與顆粒并存情況下給出的計(jì)算式能較好地?cái)M合最大單位質(zhì)量流量曲線。

2）隨開口寬度增加，截面流量及最大潰砂速率均單調(diào)增加，同等時(shí)間下突水潰砂災(zāi)害程度會(huì)逐漸加劇。 砂粒在開口處形成的密實(shí)結(jié)構(gòu)體積逐漸增加單增速不斷降低，截面流量及最大潰砂速率增幅逐漸降低，小開口時(shí)由寬度增加引發(fā)的災(zāi)害加劇大于大開口情況。 發(fā)現(xiàn)砂粒在開口處形成的密實(shí)結(jié)構(gòu)體積增加速率對(duì)災(zāi)害程度起關(guān)鍵影響，細(xì)小裂隙開口擴(kuò)張的過程中更易出現(xiàn)水砂流速突增，誘發(fā)突水潰砂災(zāi)害。

3）隨砂層厚度增加，截面流量及最大潰砂速率均單調(diào)減少，且降幅在小厚度時(shí)達(dá)最大后逐級(jí)減小至趨于平穩(wěn)，小厚度區(qū)間降幅可達(dá)大厚度的10 倍，在開口處形成的砂粒密實(shí)結(jié)構(gòu)持續(xù)時(shí)間增長，對(duì)水砂運(yùn)移的阻礙效果逐漸增加，同等時(shí)間下突水潰砂災(zāi)害程度會(huì)逐漸減弱。 不同厚度下災(zāi)害特點(diǎn)不同，厚度較小時(shí)，災(zāi)害發(fā)生在水驅(qū)動(dòng)砂快速流出裂隙孔口階段，厚度較大時(shí)，災(zāi)害發(fā)生在砂粒密實(shí)結(jié)構(gòu)失穩(wěn)造成砂粒傾瀉階段。

實(shí)際突水潰砂問題不僅需要考慮裂隙開口處的水砂流動(dòng)問題，還需要分析在裂隙網(wǎng)絡(luò)中的流動(dòng)及對(duì)裂隙網(wǎng)絡(luò)的侵蝕等相關(guān)問題。 后續(xù)將針對(duì)平面裂隙網(wǎng)絡(luò)中的流動(dòng)問題、空間破裂巖體通道中的流動(dòng)問題及裂隙通道中的滲流侵蝕問題進(jìn)行研究，為突水潰砂災(zāi)害防控提供理論指導(dǎo)。