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

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

0 引言

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

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

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

1 數(shù)值方法

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

1.1 固體對流體的作用

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

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

1.2 流體對固體的作用

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

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

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

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

δx——LBM 中格子步長；

eα——速度方向；

δt——時(shí)間步長；

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

rσ——作用力點(diǎn)到顆粒圓心的矢量。

1.3 耦合模型時(shí)間步長同步

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

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

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

2 模型優(yōu)化

2.1 二維顆粒的水力半徑

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

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

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

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

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

F——三維實(shí)體顆粒收到流體作用力；

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

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

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

本文建議LBM-DEM 耦合模型模擬多孔介質(zhì)時(shí)水力半徑因子取：

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

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

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

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

3 覆蓋型巖溶塌陷的建模

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

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

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

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

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

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

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

在實(shí)際中黏土顆粒非常小，模擬數(shù)量巨大，普通計(jì)算機(jī)能力有限而不能有效模擬。顆粒放大法是目前較為可行的處理方式，已廣泛應(yīng)用于工程研究中，該方法將顆粒放大，降低模型中的顆粒數(shù)目，使原物理模型問題能在合理有效的時(shí)間內(nèi)解決[27-30]。因此本文假設(shè)模型中的顆粒為多個(gè)黏性土顆粒結(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 萬個(gè)網(wǎng)格，1 274 391 個(gè)格點(diǎn)。固體顆粒與流體具體參數(shù)主要參照文獻(xiàn)[17]、[22]、[25]（表1）。在LBM中水頭的計(jì)算可以用等效水壓代替，令隔水層上方水壓為Pt，下方巖溶開口處水壓力為Pb，其余邊界為水量平衡邊界或不透水邊界。

表1 計(jì)算模型參數(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 覆蓋型巖溶地面塌陷的討論

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

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

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

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

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

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

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

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

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

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

圖7 塌陷過程中巖溶開口中心線顆粒受地下水垂向作用力隨時(shí)間變化曲線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（負(fù)值表示方向向下）時(shí)失穩(wěn)剝落，導(dǎo)致顆粒0588 受到的地下水垂向作用力陡增至-0.144 N 從而失穩(wěn)，進(jìn)一步導(dǎo)致顆粒0799 上地下水垂向作用力突變?yōu)?0.129 2 N。顆粒一旦開始下降后，受到向下的水流作用力逐步減小。模型表層的顆粒因水流實(shí)現(xiàn)貫通，造成所受到的水流作用力紊動(dòng)。

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

（4）位移變化分析

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

圖8 塌陷過程中巖溶開口中心線顆粒的垂向位移隨時(shí)間變化曲線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é)作用后，巖溶開口處土洞附近的顆粒在重力和水力作用下剝落。土體顆粒的剝落容易造成土顆粒原位置和上方位置處水壓的陡降，從而造成較強(qiáng)的水力坡降，使得地下水對內(nèi)部顆粒作用力陡增，容易進(jìn)一步引起上方顆粒在地下水作用力和重力作用下失穩(wěn)，從而使得承壓水位下降引起的塌陷從土體顆粒失穩(wěn)至土層塌陷過程加速。

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