陳小茜，曾 斌，王春暉，鄭束寧，彭丁茂

(1.中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074；2.杭州交通投資建設(shè)管理有限公司，浙江 杭州 310005；3.浙江省交通規(guī)劃設(shè)計研究院，浙江 杭州 310006)

滲透破壞是導(dǎo)致多種地質(zhì)災(zāi)害發(fā)生的重要原因之一，如由隧道滲水引起的地表沉降、雨水入滲導(dǎo)致的滑坡、土壩滲漏引起的壩體失穩(wěn)等[1～3]。為了減少滲流引起的災(zāi)害，復(fù)雜地層結(jié)構(gòu)中滲流問題的研究則變得十分重要。滲透系數(shù)是滲流問題中一個非常重要的參數(shù)，它反應(yīng)了地層的透水性能[4～5]。對于裂隙介質(zhì)，前人多用等效多孔介質(zhì)、雙重介質(zhì)等方法求取地層滲透系數(shù)[2～3]。對于非均質(zhì)孔隙介質(zhì)，現(xiàn)階段多采用現(xiàn)場試驗或經(jīng)驗法求取滲透系數(shù)。但現(xiàn)場試驗費用高、耗時長，且獲取的滲透系數(shù)通常是局部的、近似的；經(jīng)驗法求取的參數(shù)則偏差較大，難以滿足計算要求的精度。另外，等效滲透理論也常用于計算滲透系數(shù)[2]。對于相對簡單的地質(zhì)條件，基于等效滲透系數(shù)理論的計算方法簡潔、精度較高；但對于復(fù)雜的地質(zhì)條件，其計算方法則很難獲取精確的結(jié)果。

顆粒流程序(PFC，Particle Flow code)通過離散元法模擬規(guī)則顆粒介質(zhì)之間的運動及其相互作用，適用于固體力學(xué)大變形等問題[6～8]。顆粒流程序不但可以表征宏觀物質(zhì)的物理特性，而且可以反映其他方法無法實現(xiàn)的細觀特性。前人已用顆粒流程序在工程地質(zhì)和地質(zhì)災(zāi)害等領(lǐng)域進行過諸多探索，包括對壩基滲流、砂土滲透、堤防滲透變形等工況的模擬研究[9～13]。但尚無利用顆粒流程序模擬非均質(zhì)層狀含水介質(zhì)的等效滲透系數(shù)等方面的研究。本文采用顆粒流離散元法探討地層等效滲透系數(shù)的求取，從滲流過程的細觀角度出發(fā)，分析顆粒的運動規(guī)律及應(yīng)力分布特征，重點討論介質(zhì)結(jié)構(gòu)對滲透系數(shù)的影響；并將模擬結(jié)果與經(jīng)典理論公式計算結(jié)果對比，討論數(shù)值模擬方法的可行性與有效性。

1 顆粒流數(shù)值模擬方法

PFC中用圓形顆粒代替土顆粒，引入“域”和“管道”的概念來表征和模擬真實的流體?！坝颉笔且幌盗蟹忾]的顆粒鏈，鏈上的每個鏈接都是一個連接接觸?！坝颉贝嫠畨哼M行模擬，將壓力等效成體力施加在顆粒上，在模擬中不斷更新。假想固體中流體的通道在顆粒接觸處，并與顆粒接觸處相切，稱該通道為“管道”。利用“管道”和“域”的連接模擬顆粒與流體的耦合作用。

PFC中流固耦合計算的流動方程、壓力方程和求解條件如下[5～7]：

(1)流動方程

流體管道相當(dāng)于一個平行通道，長度為L′、孔徑為a，在垂直平面方向上為單位厚度，管道內(nèi)的流速q(單位時間內(nèi)的體積)為：

(1)

式中：K′——傳導(dǎo)系數(shù)/(cm·s-1)；

P2-P1——相鄰域的壓力差/ Pa。

(2)壓力方程

周圍管道流入每個域的流量和為∑q，在單位時間步長Δt下，流體壓力增量ΔP(流入為正)為：

(2)

式中：Kf——流體的體積模數(shù)/ kPa；

Vd——域的表觀體積/mm3。

(3)耦合方式

在處理耦合過程時，采用以下方法：①通過接觸的張開與閉合或接觸力的變化實現(xiàn)通道間隙的變化；②通過改變研究區(qū)內(nèi)域的力學(xué)特征來改變壓力；③區(qū)域孔隙壓力對其內(nèi)部顆粒有推移作用。

(4)求解方法

應(yīng)用顯示求解方法，將流量方程應(yīng)用于所有的管道，并將壓力方程應(yīng)用于所有的域之間交替求解。假設(shè)某個域內(nèi)存在擾動壓力ΔPp，由于擾動流域里的流量可以從式(1)計算得：

(3)

由水流流入引起的響應(yīng)壓力變化ΔPr為：

(4)

式中：N——域所連接的管道數(shù)；

a——管道孔徑/mm；

R——域周圍顆粒的平均半徑/mm。

(5)求解條件

保證模型運行穩(wěn)定的條件是水流入引起的壓力變化必須小于擾動壓力，當(dāng)兩者相等時可求出臨界時間步長：

(5)

(6)滲透系數(shù)K

(6)

式中：A——流體通過的橫截面積/m2;

Δh——上下游水頭差/m;

ΔL——流程長度/m。

取長度L=10 mm范圍內(nèi)所有截面的流量總和，在Δt時間內(nèi)：

(7)

則滲透系數(shù)K的表達式為：

(8)

2 滲流模型構(gòu)建

2.1 模型結(jié)構(gòu)的設(shè)計

利用PFC軟件生成圓形顆粒，模擬實際情況中的細土顆粒，滲流模型的尺寸為20 mm×20 mm。經(jīng)計算分析，當(dāng)模型尺寸較小、其他細觀參數(shù)相同時，由顆粒分布不同引起滲透系數(shù)的變化為同一數(shù)量級內(nèi)的微小變化[11]，可直接忽視。故本文中顆粒半徑采用Rmax到Rmin的均勻分布，Rmax=0.58 mm，Rmin=0.35 mm。對上、下邊界的顆粒進行固定，標(biāo)記為不排水邊界(如圖1(a)中上、下側(cè)的綠色顆粒)。通過PFC生成“域”和“管道”所組成的網(wǎng)絡(luò)。模型結(jié)構(gòu)示意圖見圖1，圖1(a)中圓形顆粒為土顆粒，圖1(b)中黑色圓點為“域”，“域”之間的淺綠色線段為“管道”，圖1(c)中藍色線段為力鏈。

模型左側(cè)水壓固定為2×105Pa，右側(cè)水壓固定為0。為減少邊界對滲流模型的影響，對文中所有模型的左、右側(cè)分別預(yù)留1.5 mm進行加壓，則實際的滲流途徑為17 mm。因無法測量某一過水?dāng)嗝娴牧髁?，故選取模型L=10 mm處的中間范圍測量總流量。

圖1 模型結(jié)構(gòu)示意圖Fig.1 Schematic diagram of the model structure

2.2 模擬參數(shù)的選取

用PFC模擬地層1、地層2，兩地層中的細觀參數(shù)——體積模量、域的管道數(shù)量等保持一致。為使模型符合實際滲流情況，經(jīng)大量數(shù)據(jù)調(diào)試可得最佳數(shù)值，參數(shù)取值如表1所示。

表1 數(shù)值模擬基本參數(shù)

PFC中的傳導(dǎo)系數(shù)是指流體在滲流通道內(nèi)的流動速度，是流體細觀運動參數(shù)。白若虛等[11]利用PFC探討滲透系數(shù)的細觀影響因素，認為滲透系數(shù)與傳導(dǎo)系數(shù)K′呈線性相關(guān)。本文通過改變傳導(dǎo)系數(shù)模擬地層1、地層2的滲透系數(shù)。結(jié)果見表2。

表2 不同傳導(dǎo)系數(shù)下的流量和滲透系數(shù)

3 模擬合理性驗證

3.1 理論等效滲透系數(shù)公式

自然界中，地下水流線通過具有不同滲透系數(shù)的地層邊界時，會像光線通過一種介質(zhì)進入另一種介質(zhì)一樣發(fā)生折射，這種現(xiàn)象是為了保持通過每個過水?dāng)嗝娴牧髁肯嗟取Ｋ鲝臐B透性較好的地層進入滲透性較差的地層，流線會變得更加稀疏，等水線間隔變小。對于層狀非均質(zhì)介質(zhì)，往往假想為均質(zhì)介質(zhì)，該介質(zhì)的水力梯度及含水層厚度和原含水層相等，Kb作為假想介質(zhì)——層狀非均質(zhì)介質(zhì)的滲透系數(shù)[13]。

在巖土體滲流問題求解中，往往應(yīng)用流線和地層成特殊角度下(0°、90°)的等效滲透系數(shù)公式較多。對于流線和地層成任意角度的等效滲透公式，何勇等[14]提出公式(9)，并采用IGW軟件進行了相關(guān)實例的分析。通過IGW軟件計算證實了此公式在非均質(zhì)含水層滲流計算中的可行性及準(zhǔn)確性。

(9)

式中：M1M2——地層1、2的含水層厚度/m；

α——等效流線和層面的角度；

α1、α2——流線和層面1、2的角度。

圖2 流線與含水層層面斜交示意圖Fig.2 Schematic diagram of the scenario with the streamline oblique to the aquifer

流線與地層1和地層2斜交，地層滲透系數(shù)見表2，水流流經(jīng)地層1，穿過兩地層界面時發(fā)生折射，角度由α1變?yōu)棣?。地層1、2厚度相同，M1=M2(圖2)。

(10)

(11)

cos2α1=0.94

(12)

cos2α2=1

(13)

由式(9)得Kb=2.39×10-6cm/s。

3.2 模擬驗證結(jié)果及其分析

由PFC軟件模擬計算可得：1區(qū)的流量Q1=4.27×105mm3/s，2區(qū)的流量Q2=4.27×105mm3/s，將Kb=2.39×10-6cm/s(傳導(dǎo)系數(shù)為1.81×10-6m/s)代

入PFC軟件計算得出模擬流量Q3=4.55×105mm3/s，此時Q1=Q2≈Q3，說明利用PFC軟件模擬層狀非均質(zhì)介質(zhì)的等效滲透系數(shù)的方法是合理的。

4 特殊工況等效滲透系數(shù)滲流模擬

對地下水流線斜交、垂直和平行層狀非均質(zhì)介質(zhì)層面等三種工況進行模擬。圖3(a)中流線和地層1的夾角為73°，M1=M2=4.8 mm。圖3(b)中，流線垂直于地層層面，M1=M2=10 mm。圖3(c)中，流線平行于地層層面，M1=M2=8.5 mm。地層1的滲透系數(shù)K1=1.32×10-5cm/s，地層2的滲透系數(shù)K2=1.32×10-6cm/s。模型左側(cè)邊界加壓2×105Pa，模型右側(cè)邊界水壓為0。紅色箭頭為流速矢量場，此模型邊界條件為定水頭邊界，雖然邊界水壓分布對水流方向稍有影響，但流速矢量場大體上垂直等水頭邊界，符合地下水動力學(xué)理論。

保持地層1、2的滲透系數(shù)不變，對厚度M1、M2取不同值進行模擬。當(dāng)循環(huán)達到一定步數(shù)時，滲流穩(wěn)定，得到滲透系數(shù)模擬值Ka，其與理論滲透系數(shù)Kb具有較高的擬合度，符合等效滲透公式。M1、M2、Ka、Kb計算結(jié)果見表3。

圖3 三種工況下的模擬結(jié)果示意圖Fig.3 Schematic diagram of the simulation results in three scenarios

模型運行初期，水壓從左到右分布不均，滲透系數(shù)越大，水壓越大；隨著模型運行時間的增加，滲流模型趨于穩(wěn)定，水壓分布逐漸均勻，且呈線性分布(圖4)。由此可知，在水力梯度相同的情況下，水壓在土層之間的傳遞時間和傳導(dǎo)系數(shù)成反比。

三種工況下模擬值Ka和理論值Kb擬合方程的斜率均接近為1，截距小，總體上PFC計算結(jié)果與公式計算結(jié)果基本一致，Ka和Kb的關(guān)系如圖5所示。由誤差分析可知(表3)，模擬值與理論值存在一定的差異，原因在于理論公式計算中直接忽略了地層邊界對計算結(jié)果產(chǎn)生的影響。

在水流和地層層面斜交的工況下，保持其他參數(shù)不變，改變模型大小，發(fā)現(xiàn)非均勻介質(zhì)地層等效滲透系數(shù)計算結(jié)果的可靠性與模型尺度有關(guān)：對于小尺度模型，可直接用PFC軟件求取等效滲透系數(shù)；但當(dāng)模型尺度較大時，地層的非均質(zhì)程度更高，建模難度較大，建議選取不同位置的地層分別建模計算[15]。也有學(xué)者建議保持地層結(jié)構(gòu)不變，對模型等比例縮放，但研究結(jié)果仍需通過實踐檢驗。因此在大尺度范圍內(nèi)，如何應(yīng)用PFC軟件精確求取等效滲透系數(shù)，尚需深入研究。

表3 不同地層厚度下的模擬滲透系數(shù)和理論滲透系數(shù)的結(jié)果對比

圖4 滲流與地層平行情況下不同時段水壓分布圖Fig.4 Schematic diagram of pressure distribution at different time intervals in the scenario with the streamline parallel to the layer

圖5 模擬結(jié)果Ka與理論滲透系數(shù)Kb關(guān)系圖Fig.5 Relationship between the simulation results Ka and theoretical coefficient of permeability Kb

5 結(jié)論

(1)采用PFC計算非均質(zhì)層狀介質(zhì)的等效滲透系數(shù)得到的結(jié)果與使用理論公式計算得到的結(jié)果基本一致，Kb與Ka關(guān)系式的斜率約為1，最大誤差為-9.61%；表明基于顆粒流理論的數(shù)值模擬方法適用于小尺度條件下等效滲流模型的計算及等效滲透系數(shù)的求取。

(2)采用PFC數(shù)值模擬方法，可得到滲流模型內(nèi)各質(zhì)點、任意時段的水壓、流速矢量場，且水力梯度相同時，水壓在土層之間的傳遞時間和傳導(dǎo)系數(shù)成反比。這一結(jié)果可為研究管涌、巖溶塌陷等的臨界水力梯度、臨界水壓提供參考依據(jù)。

(3)對復(fù)雜地質(zhì)條件下的非均質(zhì)含水層的滲透系數(shù)的求取，應(yīng)采用適用性更強的PFC進行模擬計算。

[1] 劉印,張冬梅,黃宏偉.盾構(gòu)隧道局部長期滲水對隧道變形及地表沉降的影響分析[J].巖土力學(xué),2013,34(1): 290-298.[LIU Y,ZHANG D M,HUANG H W. Influence of long-term partial drainage of shield tunnel on tunnel deformation and surface settlement[J].Rock and Soil Mechanics,2013,34(1): 290-298.(in Chinese)]

[2] 張奇華,李玉婕,袁東,等. 地下水封洞庫水幕孔注水試驗及巖體等效滲透參數(shù)分析[J]. 巖土力學(xué), 2015, 36(9):2648-2658. [ZHANG Q H,LI Y J,YUAN D,etal.Water injection test about water curtain borehole for underground water-sealed cavern and analysis of rock equivalent permeability parameter[J].Rock and Soil Mechanics, 2015, 36(9):2648-2658. (in Chinese)]

[3] 簡文星,丁志雷,李世金.巴東縣城白土坡裂隙巖體等效滲透系數(shù)數(shù)值分析[J].水文地質(zhì)工程地質(zhì),2015,42(3):65-68.[JIAN W X, DING Z L, LI S J. Numerical analysis of equivalent coefficient of permeability of the Baitupo slope in Badong County [J].Hydrogeology & Engineering Geology,2015,42(3):65-68. (in Chinese)]

[4] 張人權(quán),梁杏,靳孟貴,等.水文地質(zhì)學(xué)基礎(chǔ)[M].北京:地質(zhì)出版社. 2010. [ZHANG R Q,LIANG X,JIN M G,etal. Fundaments of hydrogeology [M].Beijing:Geological Publishing House,2010.(in Chinese)]

[5] JACOB Bear. Hydraulics of groundwater[M]. New York: Dover Publications, 2007.

[6] Itasca Consulting Group Inc.PFC2Dparticle flow code in 2 dimensions-optional features [M].Minneapolies: Itasca Consulting Group Inc,2005.

[7] Itasca Consulting Group. Verification Problems and Example Applications[M]. Minneapolies: Itasca Consulting Group Inc,2005.

[8] Cundall P A. PFC2DUser ’s Manual (Version2. 0) [R]. Minnesota: Itasca Consulting Group Inc,1999.

[9] 李識博,王常明,王剛城,等.松散堆積物壩基滲透淤堵試驗及顆粒流模擬[J].水利學(xué)報,2012,43(10):1163-1170.[LI S B, WANG C M, WANG G C,etal. Experimental study on seepage and silting and seepage of granular pile foundation [J]. Journal of Hydraulic Engineering, 2012, 43(10): 1163-1170. (in Chinese)]

[10] 張剛,周健,姚志雄.堤壩管涌的室內(nèi)試驗與顆粒流細觀模擬研究[J].水文地質(zhì)工程地質(zhì),2007,34(6): 83-86. [ZHANG G,ZHOU J, YAO Z X. Study on mesomechanical simulation of piping with model tests and PFC2D[J]. Hydrogeology & Engineering Geology, 2007,34(6): 83-86. (in Chinese)]

[11] 白若虛,程雪松,鄭剛.關(guān)于土滲透系數(shù)顆粒流細觀參數(shù)研究 [J].低 溫 建 筑 技 術(shù),2012(1):64-67. [BAI R X, CHENG X S, ZHENG G. Study on the mesoscopic parameters of soil Hydraulic conductivity [J]. Low Temperature Architecture Technology, 2012(1): 64-67. (in Chinese)]

[12] 倪小東,朱春明,王媛. 基于三維離散-連續(xù)耦合方法的堤防工程滲透變形數(shù)值模擬方法[J].土木工程學(xué)報, 2015, 48(增刊1): 159-165. [NI X D, ZHU C M, WANG Y. Analyzing the behavior of seepage deformation by three-dimensional coupled continuum-discontinuum methods [J]. China civil Engineering Journal,2015,48(Sup1): 159-165. (in Chinese)]

[13] 薛禹群.地下水動力學(xué)[M].北京:地質(zhì)出版社,1997.[XUE Y Q. Groundwater dynamics[M]. Beijing: Geological Publishing House, 1997. (in Chinese)]

[14] 何勇,袁明,陳旭,等.等效滲透系數(shù)在非均質(zhì)含水層滲流計算中的應(yīng)用[J]. 電網(wǎng)與清潔能源, 2010, 26(4): 73-76. [HE Y, YUAN M, CHEN X,etal. Application of equivalent penetration coefficient method in seepage flow calculation of non-homogeneous aquifer [J]. Power System and Clean Energy, 2010, 26(4): 73-76. (in Chinese)]

[15] 張弛,仵彥卿,覃榮高.滲透系數(shù)升尺度對非均質(zhì)含水層溶質(zhì)遷移影響研究[J].水文地質(zhì)工程地質(zhì),2014,41(5): 20-25.[ZHANG C,WU Y Q, TAN Y G. Research on effect of hydraulic conductivity upscaling on groundwater solute transport in heterogeneous aquifer [J]. Hydrogeology & Engineering Geology, 2014,41(5):20-25. (in Chinese)]