基于LBM-DEM 細觀數(shù)值模擬的水力誘導(dǎo)覆蓋型巖溶地面塌陷發(fā)育過程分析

陶小虎，葉 明，龔建師，王赫生，胡曉雨

（1.中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，江蘇 南京 210016；2.自然資源部流域生態(tài)地質(zhì)過程重點實驗室，江蘇 南京 210016；3.佛羅里達州立大學(xué)地球、海洋與大氣科學(xué)學(xué)院，佛羅里達州塔拉哈西 32306；4.江蘇省水資源服務(wù)中心，江蘇 南京 210029）

0 引言

覆蓋型巖溶地面塌陷的形成是自然環(huán)境變遷的一個組成部分，而人類活動對巖溶發(fā)育地區(qū)自然平衡狀態(tài)的干擾增多，大大加快了這一地質(zhì)災(zāi)害發(fā)生的進程，由此引發(fā)的地面塌陷問題日趨頻繁、后果嚴重，已成為巖溶地區(qū)可持續(xù)發(fā)展的一大障礙[1-2]。在眾多覆蓋型巖溶塌陷形成過程中，由地下水位變化引起的巖溶塌陷居多，地下水流提供了最主要的動力。地下水過量開采、水庫蓄水等引起的地下水位變化為主要誘發(fā)因素[3-6]。

長期以來，國內(nèi)外學(xué)者對地下水位變化引起的覆蓋型巖溶地面塌陷問題一直十分關(guān)注并開展了卓有成效的研究。在數(shù)學(xué)解析方面，基于滲透力學(xué)、極限平衡理論、牛頓摩擦定律等方法建立簡單力學(xué)模型，獲得宏觀規(guī)律，確定巖溶塌陷的機制及其臨界條件的研究[7-9]。在物理實驗方面，覆蓋型巖溶塌陷模型試驗研究從單一致塌因素研究到多因素組合、從簡單的二維定性研究到三維定量分析以及從真實土到透明土可視化分析都體現(xiàn)了試驗研究的巨大進步[10-12]。在數(shù)值模擬方面，基于連續(xù)介質(zhì)理論，在巖溶塌陷臨界條件、主控因素、塌陷模式等內(nèi)容的研究上取得了較好的成果[13-14]。然而由于受到當前科學(xué)技術(shù)的局限、土洞發(fā)育以及塌陷過程的非連續(xù)、非均勻、流固耦合等特性，現(xiàn)有的物理實驗和數(shù)值模擬方法很難對水力誘導(dǎo)覆蓋型巖溶塌陷進行系統(tǒng)的巖溶塌陷發(fā)育和塌陷過程的細觀研究，對各影響因子之間的定量關(guān)系研究也受到了限制[15]。

離散元方法（discrete element method,DEM）是根據(jù)離散物質(zhì)本身的離散特性建立數(shù)值模型，與連續(xù)介質(zhì)理論方法相比，DEM 可以從細觀角度完成大變形分析，更好地描述結(jié)構(gòu)破壞過程[16-17]。格子Boltzmann 方法（lattice Boltzmann method，LBM）是一種新興的計算流體學(xué)方法，能夠從介觀尺度上描述流體的運動和模擬多相多尺度問題。與傳統(tǒng)流體力學(xué)計算方法相比，LBM在程序?qū)崿F(xiàn)上更為簡單，能夠?qū)崿F(xiàn)與大量移動顆粒的完全耦合，很好地模擬多相多尺度耦合問題，尤其是在模擬復(fù)雜邊界條件的多孔介質(zhì)水流中，更能體現(xiàn)出其相比于傳統(tǒng)流體力學(xué)計算方法的優(yōu)勢[18]。因此，為了能夠從細觀上對覆蓋性巖溶塌陷模擬研究，本文將采用DEM模擬土體顆粒運動，選擇庫侖定律來判別土體顆粒之間的黏性接觸是否受到破壞，采用LBM 中D2Q9 模型[19]模擬土體孔隙間水流運動，應(yīng)用Bouzidi 等[20]反彈邊界格式和動量交換法實現(xiàn)LBM 與DEM 雙向耦合，建立覆蓋型巖溶塌陷二維流固耦合細觀數(shù)值模型，并對承壓水下降工況下的覆蓋性巖溶塌陷過程進行探索研究。

1 數(shù)值方法

本節(jié)主要介紹LBM-DEM 雙向流固耦合的方法，主要包括固體對流體的作用、流體對固體的作用、耦合模型時間步長同步。

1.1 固體對流體的作用

為了計算固體顆粒與流體之間的相互作用的準確性，特別是流體對固體顆粒的作用力的計算，本文采用具有二階精度且更為簡單的線性插值的Bouzidi 插值反彈邊界格式（圖1），假設(shè)每個時步末節(jié)點A處-α方向的分布函數(shù)是由虛擬流體節(jié)點D處的粒子向右遷移，經(jīng)過點C反彈運動至節(jié)點A所得，其中一個時間步長所遷移長度為一個空間步長，即 |DC|+|AC|=|AB|。取q=|AC|/|AB|，fα(r,t)為在r位置離散速度在 α方向上粒子密度分布函數(shù)。當q<時，虛擬節(jié)點D處于流體節(jié)點E、A之間，可由流體節(jié)點E、A處的分布函數(shù)插值得到虛擬節(jié)點D的分布函數(shù)；當q≥時，虛擬節(jié)點D處于流體節(jié)點A與固體節(jié)點B之間，可以通過節(jié)點A處的碰撞后的分布函數(shù)進行線性插值得到點A 處下一時步的-α方向的分布函數(shù)。

圖1 Bouzidi 插值反彈邊界示意圖Fig.1 Schematic view of Bouzidi interpolation bounce-back diagram

1.2 流體對固體的作用

動量交換法能夠精確有效且便捷的處理流體流動對固體的力學(xué)作用的計算[21]。因此本文將采用動量交換法評估二維流體對固體的力和力矩。在評估流體對固體顆粒的作用力和力矩時，固體顆粒中的分布函數(shù)不參與計算。

式中：rf——顆粒邊界流體節(jié)點位置；

Fσ——流體節(jié)點處流體對固體顆粒的作用力；

Tσ——流體節(jié)點處流體對固體顆粒的作用的力矩；

δx——LBM 中格子步長；

eα——速度方向；

δt——時間步長；

fα(r,t)——在r位置離散速度在 α方向上粒子密度的分布函數(shù)；

rσ——作用力點到顆粒圓心的矢量。

1.3 耦合模型時間步長同步

在耦合模型中，涉及到兩套方法的時間步長，作如下考慮：

當LBM 時間步長 δt小于或等于DEM 預(yù)估時間步長 δt時，DEM 時間步長取LBM 時間步長，δt=δt；

當LBM 時間步長 δt大于或等于DEM 預(yù)估時間步長 ?t時，令n=，將DEM 時間步長更新為 ?t=。

2 模型優(yōu)化

2.1 二維顆粒的水力半徑

實際覆蓋層土體中，流體能夠穿過由多個顆粒形成的多孔介質(zhì)，即使相互緊密接觸的顆粒之間也存在滲透路徑。然而在二維LBM-DEM 耦合模型中，三維顆粒會被簡化成二維圓形，相互緊密接觸的圓形之間不存在任何通道，阻止了流體的流動，這與實際情況相悖。此外，三維顆粒與二維圓形的受力也不相同。因此在二維LBM-DEM 耦合模型中需要對顆粒進行處理。

Cui[22]從流體等效拖曳力出發(fā)，將三維球體顆粒等效為圓柱體，從而得到二維圓形顆粒的水力半徑rh≈0.8R，其中R為顆粒的實際半徑，該式適用于雷諾數(shù)較大（Re>1 000）的情況。本文將推導(dǎo)較小雷諾數(shù)的情況。

當雷諾數(shù)較小時，由斯托克斯阻力公式[23]可知三維中流體作用于圓球的力為：

同樣，二維中，流體作用于圓盤的力為：

式中：U——實體顆粒相對于流體的速度；

F——三維實體顆粒收到流體作用力；

U2D——二維圓形顆粒相對于流體的速度；

F2D——二維圓形顆粒收到流體的作用力。

因此當球體顆粒在靜水中自由沉降時，阻力應(yīng)與速度成正比增大，逐漸達到勻速運動，作用于球體的重力、浮力以及阻力保持平衡，將二維和三維計算中的平衡速度進行等效，則有：

本文建議LBM-DEM 耦合模型模擬多孔介質(zhì)時水力半徑因子?。?/p>

2.2 顆粒直徑與格子步長比優(yōu)化

由于流體是在孔隙介質(zhì)中流動，并對細小的孔隙介質(zhì)產(chǎn)生力的作用，為保證流體運動和流體作用力達到精度要求，需要選擇細小顆粒影響下的格子網(wǎng)絡(luò)最優(yōu)步長[24]。定義顆粒直徑與格子步長比 β=D/δx，其中D為顆粒直徑。本文通過平板間流體對圓形阻礙物作用力的收斂性分析，以選擇計算合理的 β。數(shù)值模擬取平板長度為0.15 m，寬度為0.01 m，流體為20 °C 下的水，其密度為1×103kg/m3，動力黏滯系數(shù)μ為1.010×10-6kPa·s。左端流入壓強3.4 Pa，右端流出壓強為3.3 Pa。顆粒直徑為1 mm，位于x=0.070 m、y=0.005 m 處。對3 種不同直徑D（1.0 mm、1.2 mm、1.5 mm）的顆粒，計算分析不同β 下流體對顆粒的作用力，計算結(jié)果如圖2 所示。

圖2 顆粒所受到水流的水平作用力隨β 的變化Fig.2 The variation of horizontal force on solid particle introduced by water flow with different β

顆粒受力計算值隨著β 的增大而逐漸收斂。當β≥8 時，流體對顆粒作用力與參考值（最小步長0.05 mm工況）比較，相對誤差可控制在0.5%以內(nèi)。因此，對覆蓋層多孔介質(zhì)中進行流固耦合計算時，為確保計算的準確性和收斂性，在進行流體計算時，模型中經(jīng)過水力半徑因子縮放后的最小顆粒直徑與格子步長比需滿足β≥8。

3 覆蓋型巖溶塌陷的建模

如圖3 所示，覆蓋型巖溶塌陷的流固耦合仿真計算共包括4 個部分：DEM 初始設(shè)置、LBM 初始設(shè)置、塌陷模型初始設(shè)置、塌陷模型計算。

圖3 覆蓋型巖溶塌陷的建模流程Fig.3 Modeling flowchart of covered karst ground collapse formation

（1）DEM 初始設(shè)置：鎖定LBM 程序；確定模型的范圍和邊界，形成模擬空間；模型的左右邊界、下邊界為Wall 邊界，模型內(nèi)顆粒不能穿越此類邊界。巖溶開口在前期設(shè)置中假設(shè)為Wall 邊界，在后期修改為開放邊界。上邊界為開放邊界，即顆?？梢宰杂沙鋈朐撨吔?。在空間內(nèi)按照四邊形規(guī)則排布或者隨機生成顆粒，并賦予顆粒力學(xué)參數(shù)（楊氏模量和內(nèi)摩擦角，為了顆粒間緊密，此時不設(shè)黏性鍵）；引入局部阻尼[25-26]，在重力作用下將顆粒壓實，經(jīng)過足夠長時間的循環(huán)，消除土體顆粒內(nèi)部的不平衡力，讓其達到相對穩(wěn)定狀態(tài)。

（2）LBM 初始設(shè)置：在DEM 初始設(shè)置完成后，激活LBM 計算程序，鎖定DEM 程序（固體顆粒位置固定）；按照模擬范圍生成離散網(wǎng)格；施加水力邊界條件，計算固體顆粒間的水的流動情況，直至達到穩(wěn)定；打開DEM 程序，修改DEM 程序中的時間步長，保證DEM和LBM 之間能夠?qū)崟r的同步交換，在水力和重力共同作用下，運行足夠步數(shù)循環(huán)，使得土體達到相對穩(wěn)定狀態(tài)。

（3）塌陷模型初始設(shè)置：對已經(jīng)生成的模型中的顆粒之間存在的接觸按照參數(shù)黏性力學(xué)參數(shù)配置黏性鍵，經(jīng)過足夠長時間的循環(huán)，消除土體顆粒內(nèi)部的不平衡力，完成顆粒流模型的黏性設(shè)置（黏性鍵只配置一次，后續(xù)計算中新生成的接觸或者顆粒黏性鍵破壞后的接觸不再擁有黏性屬性）；去除模型中巖溶開口處的墻體，改為開放邊界。開口處的顆粒由于失去墻體的頂托力在重力作用和應(yīng)力釋放影響下產(chǎn)生細微變形，少部分顆粒向下脫落，并形成平衡拱，運行足夠步數(shù)，使得顆粒達到相對穩(wěn)定狀態(tài)。至此，完成對覆蓋型巖溶塌陷模型的初始狀態(tài)設(shè)置。一般情況下，拆除巖溶開口處的墻體（Wall）后的計算屬于動力計算，離散元系統(tǒng)中的阻尼應(yīng)設(shè)置為0，但由于巖溶開口的形成一般是非常緩慢的擴展過程，而模型中為了方便，將巖溶開口處的墻體直接拆除，容易造成模型系統(tǒng)內(nèi)顆粒的流失，與實際狀態(tài)不符。因此，在此步驟下可以仍然考慮繼續(xù)采用阻尼來避免這一問題。

（4）塌陷模型計算：模型完成初始狀態(tài)設(shè)置后，確保去掉阻尼（設(shè)為0）；根據(jù)計算需求修改流體邊界條件；進一步進行耦合計算至穩(wěn)定或者塌陷。（圖4）

圖4 覆蓋型巖溶塌陷LBM-DEM 模型Fig.4 Schematic view of LBM-DEM model of covered karst ground collapse

在實際中黏土顆粒非常小，模擬數(shù)量巨大，普通計算機能力有限而不能有效模擬。顆粒放大法是目前較為可行的處理方式，已廣泛應(yīng)用于工程研究中，該方法將顆粒放大，降低模型中的顆粒數(shù)目，使原物理模型問題能在合理有效的時間內(nèi)解決[27-30]。因此本文假設(shè)模型中的顆粒為多個黏性土顆粒結(jié)合的土體，以增加顆粒的半徑。如圖5 所示，模型空間范圍取15.9 cm × 8.0 cm，巖溶開口大小為0.6 cm，顆粒直徑取0.15 cm，對應(yīng)LBM 中顆粒縮放后的直徑為0.12 cm（水力半徑因子取0.8），LBM 空間步長取 δx=0.000 1 m，即β 為12，生成127.2 萬個網(wǎng)格，1 274 391 個格點。固體顆粒與流體具體參數(shù)主要參照文獻[17]、[22]、[25]（表1）。在LBM中水頭的計算可以用等效水壓代替，令隔水層上方水壓為Pt，下方巖溶開口處水壓力為Pb，其余邊界為水量平衡邊界或不透水邊界。

表1 計算模型參數(shù)Table 1 Summary of simulation model parameters

圖5 承壓水位下降工況覆蓋型巖溶塌陷過程Fig.5 The process of covered karst ground collapse formation with declining confined water level

4 覆蓋型巖溶地面塌陷的討論

為了方便模擬，本次模擬工況為初始時刻模型中各處水頭相等，等效成水壓為933.33 Pa，Pb以降速50 Pa/s下降500 Pa 的情況，模擬結(jié)果如圖5 所示。

（1）覆蓋型巖溶塌陷形成過程

在初始時刻，模型中各處的水頭相等，水流處于相對靜止狀態(tài)，對土體僅有向上的浮力作用。圖5（a）為模型完成初始狀態(tài)設(shè)置后的顆粒分布情況。

隨著承壓水的水位下降，巖溶開口處的水頭低于土層中的水頭，土層中孔隙水向下流動，對土體顆粒造成向下的力。由于水位下降初期水流對巖溶開口上方土顆粒的作用力較小，對顆粒未造成較大影響，隨著承壓水位的繼續(xù)下降，孔隙水對土顆粒的作用力逐漸增加，巖溶開口上方土體顆粒逐漸向洞口處移動，并脫離內(nèi)側(cè)相鄰?fù)馏w顆粒，見圖5（b）。如圖5（c）所示，隨著水力的作用，巖溶上方失穩(wěn)土體顆粒增多，土洞不斷擴展，但是此時黏性土層表面并不會出現(xiàn)明顯的下沉。土洞不斷擴大，在自身重力以及水力作用下，表層土體受到破壞并突然塌陷進下方的洞穴，見圖5（d），而形成覆蓋型巖溶塌陷，此后墜入土洞的土體在水力作用下逐漸通過巖溶開口流失，見圖5（e），并最終形成落水洞。

比較圖5（b）（c）（d）（e），土洞擴展由初期土體顆粒運動呈現(xiàn)單顆粒向下剝落，逐漸演變成塊狀失穩(wěn)剝落，塊體內(nèi)的土體顆粒之間僅發(fā)生極小的相對位移甚至不發(fā)生相對運動，發(fā)生失穩(wěn)破壞的塊狀構(gòu)造不斷變大，土洞擴展加速，最終導(dǎo)致土洞上方土層失穩(wěn)。隔水層中的土洞從開始擴展到土洞頂部達到極限厚度是一個逐漸加速破壞的過程，頂板失穩(wěn)塌陷是土洞擴展的特例。

（2）塌陷過程中的水力特性

在承壓水位下降的影響下，巖溶開口處的顆粒優(yōu)先發(fā)生失穩(wěn)，發(fā)生較大運動，從而影響附近的水壓分布。如圖6 所示為塌陷過程中巖溶開口附近高程水壓隨時間變化曲線，圖6 中y表示為塌陷過程中巖溶開口中心線不同高度（底部y=0 mm）。發(fā)生顆粒剝落時，顆粒向下運移，原位置和上方附近水壓有陡降，原因來自兩方面：一方面是顆粒向下運動，原位置需要得到周圍地下水補償，從而引起周圍水壓釋放；另一方面，顆粒失穩(wěn)后會在重力作用下加速下降，這個速度一般會比水流速度大得多，會對顆粒附近的局部范圍的水流造成擾動，加速水流，從而導(dǎo)致顆粒上方附近的水壓降低。這使得上方顆粒附近形成較大的水壓差，從而加大上方顆粒塌陷的可能。

圖6 塌陷過程中巖溶開口中心線不同高程處水壓隨時間變化曲線Fig.6 The variation of water pressure at different elevations of karst opening centerline during covered karst ground collapse formation

（3）地下水對顆粒作用力分析

對巖溶開口中心線上7 個顆粒（0377、0588、0799、1854、3331、4808、6285，位置從下至上排序）受到地下水垂向作用力隨時間變化作出曲線，如圖7 所示，其中0377、1854、3331、4808、6285 五個顆粒在垂向上等距，0377、0588、0799 為垂向上相鄰的顆粒。

圖7 塌陷過程中巖溶開口中心線顆粒受地下水垂向作用力隨時間變化曲線Fig.7 The variation of vertical force on the particles at karst opening centerline during covered karst ground collapse formation

從圖7 可以看出，巖溶開口附近的土體顆粒受力失穩(wěn)剝落后，導(dǎo)致內(nèi)部顆粒受到地下水的作用力明顯增大。例如，顆粒0377 在地下水垂向作用力為-0.073 2 N（負值表示方向向下）時失穩(wěn)剝落，導(dǎo)致顆粒0588 受到的地下水垂向作用力陡增至-0.144 N 從而失穩(wěn)，進一步導(dǎo)致顆粒0799 上地下水垂向作用力突變?yōu)?0.129 2 N。顆粒一旦開始下降后，受到向下的水流作用力逐步減小。模型表層的顆粒因水流實現(xiàn)貫通，造成所受到的水流作用力紊動。

由此不難推斷出，當巖溶開口上方土洞周圍存在強度較弱、甚至對土洞內(nèi)部產(chǎn)生黏性拉力土體時，該部分土體可以在較小水力作用下失穩(wěn)，但其失穩(wěn)后會造成較強的瞬間水流作用力，容易使得土洞邊沿土體失穩(wěn)。土洞發(fā)展過程中，對土體顆粒的作用力會越來越大，從而導(dǎo)致剝落力呈增大趨勢。

（4）位移變化分析

圖8 所示為塌陷過程中巖溶開口中心線不同高程處顆粒的垂向位移隨時間變化曲線。從圖8 中可以看出，巖溶開口中心線上土體顆粒從下向上逐級剝落，高程越高的土體顆粒位移變化較為一致，這說明巖溶開口上方土層從下向上，以顆粒剝落逐步演變?yōu)閴K狀剝落，直至蓋層整體失穩(wěn)，整個過程是加速過程。

圖8 塌陷過程中巖溶開口中心線顆粒的垂向位移隨時間變化曲線Fig.8 The variation of vertical displacement of the particles at karst opening centerline during covered karst ground collapse formation

5 結(jié)論

（1）承壓水位下降情況下地下水主要對巖溶開口處土層的顆粒產(chǎn)生向下的作用力。在克服土顆粒之間的黏結(jié)作用后，巖溶開口處土洞附近的顆粒在重力和水力作用下剝落。土體顆粒的剝落容易造成土顆粒原位置和上方位置處水壓的陡降，從而造成較強的水力坡降，使得地下水對內(nèi)部顆粒作用力陡增，容易進一步引起上方顆粒在地下水作用力和重力作用下失穩(wěn)，從而使得承壓水位下降引起的塌陷從土體顆粒失穩(wěn)至土層塌陷過程加速。

（2）本文所建立的流固耦合模型對進一步從細觀上研究水力驅(qū)動的覆蓋型巖溶地面塌陷的發(fā)育特征以及發(fā)育的臨界條件、定量分析各影響因子之間的關(guān)系具有理論及實際意義。