基于顆粒流理論的層狀非均質(zhì)介質(zhì)等效滲透系數(shù)研究

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

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

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

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

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

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

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

(1)流動(dòng)方程

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

(1)

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

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

(2)壓力方程

周?chē)艿懒魅朊總€(gè)域的流量和為∑q，在單位時(shí)間步長(zhǎng)Δt下，流體壓力增量ΔP(流入為正)為：

(2)

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

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

(3)耦合方式

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

(4)求解方法

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

(3)

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

(4)

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

a——管道孔徑/mm；

R——域周?chē)w粒的平均半徑/mm。

(5)求解條件

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

(5)

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

(6)

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

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

ΔL——流程長(zhǎng)度/m。

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

(7)

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

(8)

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

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

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

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

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

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

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

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

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

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

3 模擬合理性驗(yàn)證

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

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

在巖土體滲流問(wèn)題求解中，往往應(yīng)用流線和地層成特殊角度下(0°、90°)的等效滲透系數(shù)公式較多。對(duì)于流線和地層成任意角度的等效滲透公式，何勇等[14]提出公式(9)，并采用IGW軟件進(jìn)行了相關(guān)實(shí)例的分析。通過(guò)IGW軟件計(jì)算證實(shí)了此公式在非均質(zhì)含水層滲流計(jì)算中的可行性及準(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ù)見(jiàn)表2，水流流經(jīng)地層1，穿過(guò)兩地層界面時(shí)發(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 模擬驗(yàn)證結(jié)果及其分析

由PFC軟件模擬計(jì)算可得：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軟件計(jì)算得出模擬流量Q3=4.55×105mm3/s，此時(shí)Q1=Q2≈Q3，說(shuō)明利用PFC軟件模擬層狀非均質(zhì)介質(zhì)的等效滲透系數(shù)的方法是合理的。

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

對(duì)地下水流線斜交、垂直和平行層狀非均質(zhì)介質(zhì)層面等三種工況進(jìn)行模擬。圖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。紅色箭頭為流速矢量場(chǎng)，此模型邊界條件為定水頭邊界，雖然邊界水壓分布對(duì)水流方向稍有影響，但流速矢量場(chǎng)大體上垂直等水頭邊界，符合地下水動(dòng)力學(xué)理論。

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

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

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

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

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

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

圖4 滲流與地層平行情況下不同時(shí)段水壓分布圖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計(jì)算非均質(zhì)層狀介質(zhì)的等效滲透系數(shù)得到的結(jié)果與使用理論公式計(jì)算得到的結(jié)果基本一致，Kb與Ka關(guān)系式的斜率約為1，最大誤差為-9.61%；表明基于顆粒流理論的數(shù)值模擬方法適用于小尺度條件下等效滲流模型的計(jì)算及等效滲透系數(shù)的求取。

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

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

[1] 劉印,張冬梅,黃宏偉.盾構(gòu)隧道局部長(zhǎng)期滲水對(duì)隧道變形及地表沉降的影響分析[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] 張奇華,李玉婕,袁東,等. 地下水封洞庫(kù)水幕孔注水試驗(yàn)及巖體等效滲透參數(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] 簡(jiǎn)文星,丁志雷,李世金.巴東縣城白土坡裂隙巖體等效滲透系數(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] 李識(shí)博,王常明,王剛城,等.松散堆積物壩基滲透淤堵試驗(yàn)及顆粒流模擬[J].水利學(xué)報(bào),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)試驗(yàn)與顆粒流細(xì)觀模擬研究[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ù)顆粒流細(xì)觀參數(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é)報(bào), 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] 薛禹群.地下水動(dòng)力學(xué)[M].北京:地質(zhì)出版社,1997.[XUE Y Q. Groundwater dynamics[M]. Beijing: Geological Publishing House, 1997. (in Chinese)]

[14] 何勇,袁明,陳旭,等.等效滲透系數(shù)在非均質(zhì)含水層滲流計(jì)算中的應(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ù)升尺度對(duì)非均質(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)]