地下集水工程滲濾系統(tǒng)的離散元數(shù)值模擬

齊 陽，高 強，唐新軍，李曉慶

(新疆農業(yè)大學 水利與土木工程學院，新疆 烏魯木齊830052)

地下集水技術是利用所營造的滲濾系統(tǒng)(河床含水層+人工濾層 +集水構筑物)的濾凈功能，將地表江河水及含水層潛水誘滲至集水構筑物內并取出地表的一種取水技術[1]。該方法是新疆農村安全飲水以及工農業(yè)用水的常用取水方式之一。但是目前該技術缺乏系統(tǒng)、完善的研究，有關層間系數(shù)的合理取值范圍缺乏明確的規(guī)定，造成集水工程濾層設計的隨機性比較大，極易出現(xiàn)滲濾系統(tǒng)淤堵或者基土細粒大量穿過濾層進入集水廊道的問題[2]。

然而，目前地下集水工程滲濾系統(tǒng)的研究大多停留在宏觀角度，難以從本質上揭示其滲濾機理。滲濾系統(tǒng)的研究問題屬于典型的流固耦合問題，目前眾多學者[3－11]在解決流固耦合問題時多采用顆粒流的方法。周健等[3]采用顆粒流的方法在砂土管涌——濾層防治方面進行了模擬研究;李識博等[4]對松散堆積物壩基滲透淤堵試驗進行顆粒流模擬;倪小東[5]采用顆粒流的方法模擬了管涌的砂槽試驗。上述數(shù)值模擬試驗結果與相關物理試驗結果均擬合較好，充分說明了顆粒流方法解決流固耦合問題的可行性。因此，本文嘗試采用顆粒流的方法對均勻基土－濾層滲濾系統(tǒng)進行不同層間系數(shù)(D15/d85，其中 D15、d85分別為濾料、基土的控制粒徑)的系列數(shù)值模擬試驗，以驗證顆粒流方法在地下集水工程滲濾系統(tǒng)研究方面的可行性，以便用其對復雜級配基土做更深一步的研究。

1 基土－濾層滲濾系統(tǒng)數(shù)值模擬原理

1.1 流固方程

對于固相顆粒，通常采用力－位移定理和動量方程[12]來模擬顆粒的運動。

對于液相介質，采用 Navier－Stokes連續(xù)方程和運動方程[13]模擬孔隙中流體的運動，具體如下式所示:

式中:ρf為流體密度;ε為孔隙率;v為流速;εv為宏觀流速;t為時間;▽為梯度;▽p為壓力梯度;μ為流體的動力黏度;fb為單位體積內流體與顆粒之間的作用力，主要為流體與固體之間的拖曳力，可用下式表示:

估算孔隙介質滲透性的公式有很多，但其核心都是Kozeny－Carman關系式，這里采用滲透率公式[14]估算滲濾系統(tǒng)各層的滲透性:

式中:B為常數(shù)，其值常取1/180;ε為孔隙率;ˉd為顆粒平均直徑。

1.2 流固耦合過程

顆粒流程序中，流固耦合方程采用顯示差分方程求解。當前時步根據(jù)顆粒間位置按力－位移定理確定接觸力，并將該力與流體對顆粒的拖曳力疊加后根據(jù)運動方程確定下一時步顆粒新位置、速度及接觸力。流體采用劃分網(wǎng)格法計算，流域內顆粒運動引起孔隙率變化，從而造成滲透率、流速等改變，進而影響到流體對顆粒的拖曳力，繼而影響顆粒下一時步的位置，具體過程如圖1所示。

圖1 顆粒流流固耦合計算過程

2 基土－濾層滲濾系統(tǒng)數(shù)值模型的建立

本文采用Cu小于5的砂樣作為基土，并根據(jù)層間系數(shù)比(D15/d85)的需要更改濾料控制粒徑生成試樣。本次數(shù)值模擬試驗設定的層間系數(shù)有3、3.5、4、4.5、5、6 六個值，因此涉及到的濾料也為六種。基土和濾料的級配曲線如圖2所示。

圖2 基土及其濾料級配曲線

綜合考慮基土細粒含量、計算機運行速度限制等因素后，將數(shù)值模型的尺寸定為60 mm×60 mm×100 mm。數(shù)值模型中，X、Y軸向墻不能透水，主要作用是限制顆粒側向移動;Z軸負半軸墻，邊界條件設定為可透水，Z軸正半軸墻為網(wǎng)狀墻，其作用等效于允許水流和顆粒穿過的邊界條件。流體單元的尺寸定為12 mm×12 mm×10 mm，試樣部分占據(jù)5×5×10個單元，如圖3灰色網(wǎng)格線所示。模型中，基土和濾料分別占總空間的2/5和3/5。數(shù)值模型及模型參數(shù)分別見圖3和表1。

圖3 數(shù)值模型

表1 數(shù)值模型參數(shù)

模型最小尺寸 60 mm，基土最大粒徑為 3.4 mm，濾料最大粒徑是8 mm～13.65 mm，模型最小尺寸與基土、濾料最大粒徑的比值分別為17.65、4.39 ～7.5，符合相關試驗標準[15]中不得小于 4 的規(guī)定。為防止基土、濾料顆粒在生成過程中劇烈混摻，生成一墻將模型沿Z軸方向分成兩個空間，在左側空間生成基料，循環(huán)消除內部不均衡力后，再在右側空間生成濾料，再經(jīng)過循環(huán)消除不均衡力后，刪除中間墻體。最后，在Z軸向砂樣兩端施加20 kPa的壓力差(確保水力梯度遠大于發(fā)生管涌的臨界水力梯度)，并保持右端壓力為零，使砂樣沿Z軸向產(chǎn)生一定的壓力梯度，模擬實際情況下的沿重力方向的滲流。

3 模擬結果分析

利用流體單元參數(shù)及程序內置FISH語言設置監(jiān)測函數(shù)，獲取基土中心、濾層及兩者界面處的孔隙率、滲透率的變化情況，統(tǒng)計基土流失率、侵蝕率及濾層各處基土顆粒的滯留情況，以探究系統(tǒng)水力特性變化規(guī)律?；裂豘軸方向平均分成10層，將第3層和第5層分別作為基土和界面的監(jiān)測層，濾層的參數(shù)取6～10層的均值?；翆娱g系數(shù)D15/d85表示濾層設計準則，D15為濾料控制粒徑，d85為基土控制粒徑。

3.1 不同層間系數(shù)對應的滲濾系統(tǒng)中基土顆粒的移動狀態(tài)

圖4給出了不同層間系數(shù)情況下模擬試驗結束時數(shù)值模型中顆粒的分布情況。其中，每個層間系數(shù)同時給出隱去濾料顆粒和不隱去濾料顆粒時數(shù)值模型中的顆粒分布狀態(tài)。

圖4 不同D15/d85值下模擬試驗結束時滲濾系統(tǒng)中基土顆粒的分布狀態(tài)

由圖4可知，當 D15/d85不大于3.5的時候，基土顆粒無明顯整體移動，極少量的基土細粒進入并滯留在濾層靠近界面處;當 D15/d85介于4～5之間的時候，基土顆粒出現(xiàn)明顯的整體移動，主要體現(xiàn)在，界面處基土顆粒進入濾層之后，上游基土顆粒整體向下移動填充界面處位置，出現(xiàn)不同程度的沉降;當D15/d85大于5的時候，形成集中的滲流通道，基土顆粒沿滲流通道持續(xù)移動，逐漸發(fā)展成管涌破壞。

總體而言，對于均勻基料，滲流過程中顆粒的移動規(guī)律主要表現(xiàn)為沿著滲流方向整體移動，滲透變形表現(xiàn)為較大的沉降變形。這主要是由均勻基料的級配決定的，均勻基料沒有骨架顆粒和填充顆粒之分，所有顆粒都參與骨架結構的形成，基土顆粒間孔隙較細，且較細顆粒和較粗顆粒粒徑相差不大，所以不存在細顆粒沿粗顆?？紫兑苿拥默F(xiàn)象，因此當界面處顆粒進入濾層之后，上游顆粒就會整體向下移動。

3.2 不同D15/d85值下滲濾系統(tǒng)各處孔隙率的變化

圖5給出了不同層間系數(shù)下滲濾系統(tǒng)基土、濾層及兩者界面處的孔隙率隨滲流作用時間的變化曲線。六種情況下，濾層的孔隙率皆先遞減然后再保持穩(wěn)定。但是，隨著D15/d85值的增加濾層的孔隙率的減小量加大，這是因為D15/d85值越大，濾料越粗，

圖5 基土、濾層及兩者界面處的孔隙率隨滲流作用時間的變化曲線

進入濾層的基土顆粒就愈多，孔隙率就愈小?；林行暮徒缑嫣幍目紫堵识际窍仍黾雍蟊３址€(wěn)定，并且D15/d85越小，保持穩(wěn)定的速度越快。但是對于濾層的孔隙率，層間系數(shù)較大的時候，滲流初期波動較劇烈，這是因為D15/d85較大時有較多的基土顆粒穿過界面進入濾層，因此在系統(tǒng)穩(wěn)定之前孔隙率處于波動狀態(tài)。整體而言，當D15/d85≤5時，滲濾系統(tǒng)各處孔隙率最終趨于穩(wěn)定，說明上述滲濾系統(tǒng)最終都保持穩(wěn)定，不會發(fā)展成管涌。此外，模擬試驗結束時，界面處孔隙率明顯大于基土和濾層，且不隨滲流作用時間減小，充分說明該情形下均勻基料不會出現(xiàn)在界面處淤堵的情形;當層間系數(shù)大于5時，界面處的孔隙率持續(xù)減小，說明模擬試驗結束時滲濾系統(tǒng)仍處于不穩(wěn)定狀態(tài)。

3.3 不同D15/d85值下滲濾系統(tǒng)各處滲透率的變化

圖6、圖7分別給出了不同層間系數(shù)下滲濾系統(tǒng)中基土、界面處的滲透率，濾層的滲透率隨滲流作用時間的變化曲線。由滲透率的計算公式可知，滲透率的大小主要與孔隙率和監(jiān)測范圍內的顆粒平均粒徑有關。對于均勻基料，由于沒有細料遷移現(xiàn)象，

孔隙率不會有太大變動，又由于基土中心和界面處顆粒平均粒徑相差不大，所以該情況下，基土滲透率K基和界面處滲透率K界的變化曲線與孔隙率的變化曲線接近。由于界面處有顆粒較大的濾料顆粒，平均粒徑會大于基土中心處，因此對于任一D15/d85值皆有K界＞K基。

6種情況下，濾層的滲透率都是先遞減后保持穩(wěn)定。初始滲透率隨著D15/d85的增加而增大，但是穩(wěn)定后的滲透率隨著D15/d85的增加而減小，甚至會出現(xiàn)低于界面處的情況。這主要是由于D15/d85值較小的時候，基土顆粒難以或者少量進入濾層，濾層保持了較大的滲透率，但是 D15/d85值較大時，大量基土顆粒進入濾層造成濾層孔隙率和顆粒平均粒徑顯著降低，致使?jié)B透率顯著下降。

由上述分析可知，對于均勻基料，當 D15/d85的取值能夠保持滲濾系統(tǒng)穩(wěn)定的時候，始終有K基＜K界＜K濾，但是在層間系數(shù)大于5的時候，由于濾層滯留的基土顆粒過多，會出現(xiàn)滲透性低于界面處的情況，應予以避免。

圖6 基土、界面處的滲透率隨滲流作用時間的變化曲線

3.4 不同D15/d85值下滲濾系統(tǒng)濾層各處侵蝕基土顆粒量的變化

將數(shù)值模型濾層部分沿滲流方向平均分成三層(根據(jù)靠近基土的程度，依次為第一層、第二層、第三層)，統(tǒng)計出各層的侵蝕基土顆粒的總體積隨滲流作用時間的變化情況，以便分析不同時刻濾層各層侵蝕基土顆粒的分布情況，如圖8所示。由圖8可知，進入濾層的基土顆粒主要滯留在第一層，并隨著D15/d85的增加，呈現(xiàn)逐漸增多的趨勢;第二、三層滯留的基土顆粒量很少，隨 D15/d85的增加略微增大，但D15/d85＞5時，呈現(xiàn)顯著增加的趨勢。

由此可知，D15/d85≤5時，進入濾層的基土顆粒通常會滯留在濾層中，但這時由于進入濾層的顆?？偭枯^少不會造成濾層淤堵;D15/d85稍大于5時，由于濾層第一層滯留的基土顆粒過多，會造成濾層出現(xiàn)明顯淤堵的現(xiàn)象。

圖7 不同層間系數(shù)下濾層的滲透率隨滲流作用時間的變化曲線

圖8 不同D15/d85值下濾層各分層侵蝕基土顆粒量的變化曲線

3.5 不同D15/d85值下進入濾層和穿出濾層的基土百分量的變化

圖9 、圖10分別給出了不同層間系數(shù)下進入濾層的基土百分量和穿出濾層的基土百分量隨滲流作用時間的變化曲線，其中進入濾層的基土量又稱為侵蝕基土量。從圖9、圖10中可以明顯看出，對于均勻級配基土，在層間系數(shù)小于或者等于5的情況下進入濾層的侵蝕基土顆粒量曲線最終都趨于水平，并且穿出濾層的侵蝕基土百分量亦趨于平緩，且保持在一個較小的范圍內，說明在該層間系數(shù)范圍內，滲濾系統(tǒng)最終會趨于穩(wěn)定，不會發(fā)展成為管涌;層間系數(shù)大于5時，進入濾層的侵蝕基土顆粒量和穿出濾層的侵蝕基土量在滲流后期仍處于增長的趨勢，且模擬結束時侵蝕基土量已達到了16%，穿出濾層的侵蝕基土量也已經(jīng)是層間系數(shù)不大于5時的數(shù)倍，因此勢必會發(fā)生較大的沉降變形，所以認為該情況下系統(tǒng)不穩(wěn)定。

圖9 進入濾層的基土百分量變化曲線

圖10 穿出濾層的基土百分量變化曲線

綜合上述分析可知，對于均勻級配基土，在設定的水力梯度下，對系統(tǒng)的保土性分析后發(fā)現(xiàn)，D15/d85≤5時，滲濾系統(tǒng)流失的顆粒較少，濾層的保土功能良好;D15/d85＞5時，滲濾系統(tǒng)會流失較多的基土顆粒，形成較大的沉降變形，濾層的保土功能較差，滲濾系統(tǒng)瀕臨失穩(wěn)的邊緣。該結果與太沙基準則基本一致，說明采用顆粒流的方法對基土－濾層滲濾系統(tǒng)進行數(shù)值模擬研究是可行的。對系統(tǒng)的透水性分析后發(fā)現(xiàn)，D15/d85在3～5之間時，系統(tǒng)始終有K基＜K界＜K濾，滲濾系統(tǒng)透水功能良好。對濾層侵蝕基土顆粒的滯留情況分析后發(fā)現(xiàn)，隨著D15/d85的增加，侵蝕基土量遞增，但D15/d85≤5時，由于進入濾層的侵蝕基土量有限不會造成濾層淤堵。

4 結論

(1)顆粒流方法，在基土－濾層滲濾系統(tǒng)模擬方面具有獨到的優(yōu)勢，它不僅可以對各粒徑區(qū)間顆粒的移動過程進行直觀的觀測，還可以監(jiān)測滲濾系統(tǒng)各處的孔隙率、滲透率，基土流失百分量等參數(shù)的變化過程，從而實現(xiàn)對滲濾系統(tǒng)的精準分析。

(2)層間系數(shù)D15/d85的取值是決定濾層是否有效及滲濾系統(tǒng)水力特性優(yōu)劣的主要因素，當D15/d85≤5時，滲濾系統(tǒng)保土功能和透水功能都比較好;當D15/d85＞5時，滲濾系統(tǒng)保土功能不足，有集中滲漏通道形成，系統(tǒng)失穩(wěn)。由此可知，對于地下集水工程中的濾層，若含水層為均勻基土，宜控制層間系數(shù)D15/d85不大于5。

(3)數(shù)值試驗結果與太沙基準則的一致性，充分說明了采用顆粒流的方法對基土－濾層滲濾系統(tǒng)進行數(shù)值模擬研究是可行的，結果是可信的，值得用其對滲濾系統(tǒng)進行更深入的研究。

[1]吳正淮.滲渠取水[M].北京:中國建筑工業(yè)出版社，1981.

[2]李曉慶，高 強，唐新軍.細沙河道上修建地下集水廊道的工程實踐[J].人民長江，2013，44(17):29-32.

[3]周 健，周凱敏，姚志雄，等.砂土管涌 －濾層防治的離散元數(shù)值模擬[J].水利學報，2010，41(1):17-24.

[4]李識博，王常明，王剛城，等.松散堆積物壩基滲透淤堵試驗及顆粒流模擬[J].水利學報，2012，43(10):1163-1170.

[5]倪小東，王 媛，王 飛.管涌的砂槽試驗研究及顆粒流模擬[J].四川大學學報:自然科學版，2009，41(6):51-57.

[6]劉先珊，劉 洋.考慮流固耦合效應的飽和砂土滲流破壞數(shù)值模擬[J].蘭州大學學報:自然科學版，2013，49(5):633-643.

[7]游碧波，周翠英.雙層堤基條件下管涌逸出的顆粒流模擬[J].中山大學學報:自然科學版，2010，49(6):42-49.

[8]鄭 剛，戴 軒，張曉雙.地下工程漏水漏砂災害發(fā)展過程的試驗研究及數(shù)值模擬[J].巖石力學與工程學報，2014，33(12):2458-2471.

[9]周 健，姚志雄，張 剛.砂土滲流過程的細觀數(shù)值模擬[J].巖土工程學報，2007，29(7):977-981.

[11]劉 洋，周 健，付建新.飽和砂土流固耦合細觀數(shù)值模型及其在液化分析中的應用[J].水利學報，2009，40(2):250-256.

[12]Itasca Consulting Group Inc. PFC3D Version4.0[M].USA:Itasca，2008.

[13]Itasca Consulting Group Inc. PFC3D Particle Flow Code in 3 Dimensions-Optional Features[M]. USA:Minneapolis，2005.

[14]JACOB Bear. Hydraulics of Groundwater[M]. McGrao－ Hill，1979.

[15]中華人民共和國水資源部.SL237－057－1999粗顆粒土的滲透及滲透變形試驗[S].北京:中國水利水電出版社，1999.