考慮土石界面滲透性的土石混合體滲流有限元模擬

何玲麗，田東方，王世梅，陳 勇

(1.水利與環(huán)境國家級實驗教學示范中心，湖北 宜昌 443002；2.三峽大學水利與環(huán)境學院，湖北 宜昌 443002；3.三峽庫區(qū)地質災害省部共建教育部重點實驗室，湖北 宜昌 443002；4.三峽大學土木與建筑學院，湖北 宜昌 443002)

1 研究背景

土石混合體(SRM)是指第四紀以來形成的，由具有一定工程尺度、強度較高的巖塊、細粒土體及孔隙構成且具有一定含石量的極端不均勻松散巖土介質系統[1]。這類介質在我國各地廣泛分布于邊坡、地基及地下工程中，還經常作為工程材料應用于大壩、岸堤、路基、機場等。以SRM的細觀結構(即其中塊石的形狀、尺寸、分布等信息)為基礎，通過滲流數值模擬可以確定水分在土體與塊石中的運移過程。以此為基礎，可開展?jié)B透數值實驗，揭示細觀結構對SRM滲透特征的影響，也可從細觀尺度研究土石混合體內部變形破壞過程，有助于探索土石混合體破壞機制。

目前，基于細觀結構的滲流模擬，多將SRM視為土體和塊石構成的二元介質，采用達西定律描述土體和塊石的滲流，利用有限元實現數值模擬，如徐文杰等[2]和Yan等[3]的研究。數值實驗表明，SRM滲透系數受含石率影響較大，且隨含石率的增大而單調遞減；同時，由于土體和塊石的滲流符合達西定律，故SRM的滲流也符合達西定律。

然而，實驗表明SRM的滲透系數隨含石率的增大并非單調遞減，而是先降后增[4-5]，如圖1a所示；此外，SRM的滲透速度是水力梯度的冪函數[6-7]，如圖1b所示。從圖1a可見，現有的數值模擬不能很好地反映SRM的滲透特征。這可能是因為未考慮土石界面滲透性而造成的。土石界面通常是SRM中的薄弱環(huán)節(jié)和應力集中的部位，這將使得土石界面附近土體的孔隙分布與基質土體存在一定差異而容易形成集中滲流通道。因此，土石界面的滲透性必須加以考慮。

圖1 土石混合體滲透特性

由于土體和塊石滲流符合達西定律，那么SRM滲透速度與水力梯度之間的冪函數關系，應由土石界面的滲流引起。結合文獻[4，6]的實驗規(guī)律，本文假定界面流量是水力梯度的冪函數。

綜上，基于SRM細觀結構進行滲流數值模擬時，現有方法未考慮土石界面的滲透性，不能正確反映SRM的滲透特性。為此，本文忽略巖塊的滲透性，將SRM視為由土體和土石界面構成的二元介質，假定土體滲流為二維達西流，土石界面為一維滲流且流量是水力梯度的冪函數，基于有限元法實現考慮土石界面滲透性的SRM滲流數值模擬。

2 模型的構建

2.1 控制方程

SRM滲流模擬區(qū)域如圖2所示，考慮到巖塊的滲透性遠小于土體，故忽略巖塊的滲透性；則滲流求解域為土體區(qū)域Ω，以及M條土石界面Γr，r=1，2，…，M。假定土體各向同性且為二維達西滲流；土石界面為一維滲流。

圖2 SRM滲流模擬區(qū)域示意

土體滲流方程為

(1)

式中，H為總水頭；ks為土體飽和滲透系數；x，y為坐標，位于水平面內。

在定水頭邊界ΓH上

(2)

沿任意一條土石界面，建立弧坐標l；設在第r條界面Γr上，從土體到界面的流量為qi(l)，則

(3)

式中，l為沿界面的弧長坐標；nx、ny分別為Γr外法線單位向量在x、y軸的分量(土體為內部區(qū)域)。

記q(l)為土石界面Γr的沿程流量，假定只考慮從土體流入界面的水量，根據水量守恒有

?q/?l=qi

(4)

本文假定q為水力梯度J的冪函數，則

q=aJb

(5)

式中，a、b為土石界面滲透參數。

記h(l)為土石界面Γr上的沿程水頭，則J=-?h/?l，令kq=aJb-1，結合式(4)、(5)有

(6)

2.2 方程的Garlerkin弱解形式

(7)

由格林公式可知

(8)

將式(8)代入式(7)，引入邊界條件(3)得

(9)

式(6)的加權余量格式為

(10)

根據格林公式，式(10)左端變?yōu)?/p>

(11)

式(11)右端第一項，當Γr自行封閉時(如圖2中的Γ2)，其值為0；當Γr不封閉時(如圖2中的Γ1)，其值根據邊界條件式(2)確定。這一條件在有限元求解時，通過引入定水頭邊界得以滿足。故該項略去不寫。

將式(10)引入式(11)后，可改寫為

(12)

2.3 模型的有限元求解

用有限元求解式(9)和式(12)時，可借鑒巖體滲流問題中的兩域滲流模型，即分別求解兩場(H和h)、迭代確定qi。但土石界面數量較多，迭代將帶來較大的計算開銷以及收斂問題。為此，本文假定水頭在土體與界面接觸處是連續(xù)的，即在任意界面Γr上，h=H。將式(9)與式(12)左右對應相加消去qi，可得

(13)

如圖3所示，本文采用適應性強的三角形單元離散Ω；土石界面將被離散為兩節(jié)點線單元，節(jié)點整體編號與三角形單元共用。

圖3 土體和界面單元示意

三角形單元e內的He、線單元τ上的Hτ采用單元節(jié)點水頭Hi、Hm的線性插值，即

(14)

([D]+[R]){H}={f}

(15)

本文采用增量法[8]求解式(15)。

3 模型驗證

3.2 算例一

算例一與徐文杰等[2]的滲透數值實驗進行對比。試樣細觀結構見文獻[2]中圖6；本文計算時，不考慮塊石滲透性，土石界面滲透參數取a=b=0，其他條件同文獻[2]。對比結果見表1，其中，ku為試樣等效滲透系數與土體滲透系數之比。

表1 不同含石率相對滲透系數ku的對比

由表1可知，本文與文獻[2]結果相差均小于2%，表明本文方法和程序可以模擬不考慮界面滲透性的情形。由于塊石透水性相比土體而言通常小幾個數量級，而塊石數量將隨含石率增大而增多從而降低試樣滲透性。因此，當只考慮土體和塊石的滲透性時，SRM的滲透性必然隨塊石的增多而減小。

3.3 算例二

本算例與王宇等[4,6]的SRM滲透性實驗進行對比。實驗試樣尺寸、邊界條件如圖4所示。試樣中塊石粒徑2～5 mm。文獻[6]中試樣的含石率分別為20%、30%、40%、50%、60%、70%；由于未提供試樣細觀結構，本文利用文獻[9]方法生成相應試樣的細觀結構；每個試樣中塊石粒徑分別為2、3、4、5 mm，各粒徑塊石的體積分數均為0.25。限于篇幅，圖4只展示了含石率為20%和70%時試樣的細觀結構。

圖4 算例二的試樣尺寸、邊界及細觀結構(灰色區(qū)域代表巖塊)

圖5 平均流速-水力梯度關系對比

數值模擬時，土體ks=6.17×10-7m·s-1，土石界面滲透參數a、b通過試算確定，具體取值見表2。

表2 算例二計算參數及實驗結果的偏差

數值模擬所得試樣中的平均流速與文獻[6]的對比見圖5。本文定義err描述某一含石率時的n組模擬和實驗數據之間的偏差，按式(16)計算，即

(16)

從圖5可以看出，本文方法可再現SRM的非達西流特性；從表2所列誤差err來看，含石率較高(如60%、70%)時，模擬結果與實驗結果較接近，含石率為30%～50%時偏差較大。

文獻[6]給出了水力梯度等于50時，滲透系數隨含石率的變化規(guī)律。數值模擬所得滲透系數與文獻[6]對比見表3。

由表3可知，盡管有時數值模擬和實驗結果相差較大(最大約35%)，但本文方法能再現SRM試樣的滲透性隨含石率的增大而先降后增的現象。

圖6 試樣平均流速與水力梯度的關系

表3 滲透系數結果對比

4 土石界面滲透參數影響分析

圖6為a、b取不同數值時，試樣平均流速與水力梯度的關系。由圖6可知，b值越大，非線性越強；a值越大，滲透速度越大；含石率越高，滲透特性對參數變化越敏感，即隨a值的增加滲透速度增大得更快，隨b值的增加非線性更強。這是因為含石率越大，土石界面數量越多，對試樣滲透特性影響越大。

本節(jié)初步研究了土石界面滲透參數(a、b)對SRM滲透性的影響。其他因素的影響將在后續(xù)研究中進一步考慮，例如土石界面流量與水力梯度的其他函數形式、塊石形狀及分布的隨機性等。

5 結 論

(1)將SRM視為由土體和土石界面構成的二元介質，假定土體滲流為二維達西流，土石界面為一維滲流且流量是水力梯度的冪函數，水頭在土體與土石界面處連續(xù)；從而避免采用迭代方法確定土體和界面間的流量交換，構建了土體和土石界面滲流場同步求解有限元模型，實現了考慮土石界面滲透性的土石混合體滲流模擬。

(2)算例驗證表明，所建方法可模擬不考慮或考慮土石界面滲透性的情形。當考慮土石界面滲透性時，本文方法可較好地再現土石混合體滲透性隨含石率增大先降后增以及滲流速度與水力梯度的非線性關系等實驗現象。

(3)滲透數值實驗表明，土石界面滲透參數對SRM滲透特性影響為：a值越大，滲透速度越大；b值越大，非線性越強；含石率越高，SRM滲透特性對參數的變化越敏感。