基于LBM-DEM 耦合方法的突水潰砂運移規(guī)律研究

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

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

0 引 言

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

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

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

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

1 突水潰砂分析模型

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

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

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

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

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

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

2.1 格子Boltzmann 方法（LBM）

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

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

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

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

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

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

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

圖3 D3Q19 模型Fig.3 D3Q19 model

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

2.2 離散單元法（DEM）

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

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

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

2.3 浸入移動邊界法（IMB）

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

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

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

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

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

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

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

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

2.4 計算流程

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

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

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

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

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

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

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

表1 模擬方案Table 1 Simulation schemes

3 計算結(jié)果分析

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

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

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

圖7 不同邊界壓力下水砂運移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ā)生變化，此時在水流驅(qū)動作用下砂粒將在裂隙出口處形成臨時的拱狀受力結(jié)構(gòu)，顆粒間相互作用使砂粒更加密集，對流場有一定的阻礙作用。隨著砂粒的不斷流出，出口處的拱狀結(jié)構(gòu)逐漸變化，其中心處不斷凹陷。 至第④步時拱狀結(jié)構(gòu)變?yōu)轭愃啤癕”狀結(jié)構(gòu)，仍對水流及流砂有阻礙作用。 直至第⑤步時，中心砂粒完全流出喉道口，拱狀結(jié)構(gòu)完全消失后不再對流場起到阻礙作用，此時流量會快速上升，而由于砂粒已大量流出其潰砂速率會快速下降。

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

為表征不同邊界壓力對水砂流動的影響程度，統(tǒng)計單位時間潰砂速率最大值及初始階段受砂粒阻礙作用的平均流量隨壓力的變化規(guī)律如圖11 所示。其中隨著邊界壓力的增加，最大潰砂速率和平均流量均單調(diào)遞增。 計算最大潰砂速率及初始平均孔口流量隨壓力的區(qū)間增幅，并通過區(qū)間增幅除以區(qū)間范圍得出單位壓力增幅以表征不同壓力數(shù)值下的影響情況，計算結(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)，低壓力時隨壓力的增加最大潰砂速率和平均流量的增幅均大于高壓力時。 說明僅考慮邊界壓力情況下，隨著邊界壓力的增加，同等時間下突水潰砂災害程度會逐漸加劇，但低壓力下由于壓力增加引發(fā)的災害加劇程度較大，此時應注意由壓力增加引起的劇烈波動。 由于模型中設(shè)置恒定壓差表征地下承壓含水層作用，未考慮實際流動中壓力降低的作用，結(jié)論僅能在一定程度上適用于工程實際問題，后續(xù)將針對復雜模型繼續(xù)研究。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

為表征不同開口寬度對水砂流動的影響程度，圖17 中給出單位時間潰砂速率最大值及初始階段受砂粒阻礙作用的平均流量隨開口寬度的變化規(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)，開口寬度較小時隨寬度增加對最大潰砂速率和初始平均流量的影響大于開口寬度較大時。 說明固定的砂含量下，隨著裂隙開口寬度的增加，突水潰砂災害程度的增幅逐漸減小，在寬度大于一定值后，由于砂含量的限制災害逐漸趨于恒定。 突增現(xiàn)象最容易出現(xiàn)在細小裂隙開口擴張的過程中。

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

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

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

圖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ā)生變化。 而隨著砂層厚度的增加，截面流量在達到初始峰值后均出現(xiàn)一定的下降，用對應曲線顏色的虛線記錄最終時間，①表示初始流量峰值點，②③④⑤分別對應N-2、N-3、N-4、N-5 時最終下落結(jié)束的時間點，發(fā)現(xiàn)下降持續(xù)的時間隨著砂層厚度的增加逐漸增長。 說明隨著砂層厚度的增加，下落過程中對水流的阻礙作用越明顯。

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

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

從圖20 可知，對水流及整體砂粒流速起阻礙作用的原因主要為砂粒在裂隙開口處聚集形成密實區(qū)域，對比各階段可以發(fā)現(xiàn)①②③④時刻大量砂粒受力后聚集在裂隙開口，使得水砂流動減緩，在⑤時刻砂粒形成的密實結(jié)構(gòu)崩塌，對流體的阻礙消失，截面流量逐漸開始上升。 而此時由于密實結(jié)構(gòu)崩塌后大量砂粒傾瀉將造成此時的潰砂速率達到峰值。 為明確不同砂層厚度對水砂運移的影響，圖21 分析不同砂層高度下單位時間潰砂速率最大值及初始階段受砂粒阻礙作用階段平均流量的變化規(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），在砂層厚度較小時出現(xiàn)均快速降低的現(xiàn)象。 當砂層厚度較大時，隨砂層厚度增加對最大潰砂速率和平均流量的影響較小。 為定量研究砂層厚度對水砂運移的影響程度，表4 計算給出不同砂層厚度下最大潰砂速率和平均流量的區(qū)間降幅。

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

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

4 結(jié)論及展望

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

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

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

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