基于COMSOL Multiphysical的蓄水運營期堤防工程安全穩定性研究

林俊屹

(雙江鎮水利局，湖南 懷化 418000)

1 概 述

堤防工程作為防洪樞紐重要的組成部分，研究蓄水運營期堤防安全穩定性具有重要作用，特別是堤防內滲流場演化特征以及邊坡失穩特性，均一定程度反映了堤防的安全穩定性[1-3]。一些專家與學者采用物理模型手段，設計模型試驗，研究在潰壩、地震荷載破壞、失穩破壞等過程中水工建筑滲流場、應力場、變形場演化特征，為水利設計提供重要參考[4-6]。利用微震或聲發射以及原位儀器，對現場運營狀態下的水工結構開展細觀特征分析，研究各個微觀數據之間關系，為評價水工建筑安全運營狀態提供依據[7-9]。以數值仿真手段，分析壩體全過程滲流及應力變形演化特征，在流固耦合、動力響應等研究中提供模型解，進而提升水工設計水平[10-12]。本文利用COMSOL Multiphysical仿真軟件[13-14]研究堤防工程滲流場演化特征以及邊坡安全穩定性，重點探討蓄水運營期下堤防安全穩定演化特征，為堤防壩體安全運營設計提供參考。

2 工程介紹及模型建立

2.1 工程概況

某防洪樞紐工程是地區內調水、控水、蓄水等復合型水資源調度水利設施，為地區供水以及水資源凈化修復提供重要作用，年可供水量超過300×104m3，按照50年一遇洪水位設計攔水大壩，確保水資源安全高效利用。該樞紐工程包括泄洪閘、抽水泵站、攔污柵以及堤防工程等水工結構，其中泄洪閘設計為多孔式水閘，單孔尺寸為5 m×3 m，閘門采用弧型鋼閘門，半徑為1.8 m，以液壓式啟閉機作為開度控制設施，確保通行流量對水工結構以及下游用水工程不產生水力威脅。攔污柵與堤防工程一起構成上游水位重要的水利控制結構，堤防長度超過1 000 m，頂寬度為3 m，主要為砂土堆料及細礫石土，含水量較低，部分細礫石土沉降變形較大，室內初步測試彈性模量遠高于普通堆石料體；堤防表面鋪設土工布，坡腳處設置厚70 cm防滲墻，共同構成堤防防滲系統。現由于該堤防工程運營年限較長，部分堤段出現滲流活躍，下游河坡的混凝土灌漿設施出現局部失穩現象，分析認為該堤防運營穩定性欠佳，局部堤身應進行加固維護。為此，工程設計部門討論先利用COMSOL Multiphysical多物理場軟件建立堤防工程數值模型，并根據實際工程水位運營高程及相關荷載條件，研究有關堤防安全的邊坡、滲流場狀態等特征參數，為評價堤防安全運營及加固除險提供重要參考。

2.2 模型建立

本文重點分析的堤段為K2+120～K2+128斷面處，上下游坡度比均為1/3.5，邊坡軸線長度為15 m，堤腳水平延伸至其他坡腳處，此設定為模型工程的邊界；上覆土層為砂土，該區段內礫石料沉降變形較大，局部還存在下臥粉土，對堤防安全具有重大威脅，因而模型試驗也以該斷面為基礎開展分析研究，并逐步擴展至其他區段內，該區段堤防斷面圖見圖1。利用COMSOL Multiphysical多物理場分析軟件建立堤防研究區段的數值模型圖，以四面體單元作為網格劃分基本單元，共劃分出網格單元15 625個，節點數14 258個，所建立的數值模型及網格模型見圖2。根據工程現場實際勘測地質資料得知，堤防斷面主要土體涉及細礫石、砂土以及黏性土等，因而相關物理參數均以室內土工試驗所測定，并寫入COMSOL軟件中計算；本文堤防邊坡安全系數計算方法采用簡化畢肖普法。為比較不同水位對堤防穩定性影響，筆者在進行仿真分析時設計高水位的不同運營時間方案，分別為0、5、10、15、20和25 d；由于堤防堆筑材料差異性會對滲流穩定性產生較大影響，特別是砂土的存在具有顯著差異性，因而本文堤防工程斷面分析方案分為多砂斷面(砂土占比50%)與全砂斷面(全斷面均為砂土)；上游水位設定高程為56.5 m；砂土、粉土、礫石土的滲透系數分別取值為3×10-5、2.6×10-3和1.85×10-2cm/s，在上述模型建立基礎上開展堤防邊坡安全穩定性以及滲流場分析。

圖1 堤防斷面圖

圖2 數值模型及網格模型圖

3 高水位浸泡下堤防邊坡穩定性

3.1 物理力學特征參數變化

堤防邊坡穩定性與土體物理力學特征參數息息相關，研究物理力學特征參數受高水位浸泡作用影響效應對認識邊坡穩定性具有重要意義。圖3為不同運營時間下多砂斷面中黏聚力、內摩擦角兩個特征參數隨運營時間變化的關系曲線。從圖3中變化可知，隨運營時間增大，堤防邊坡土質黏聚力顯著降低，浸泡天數0 d時黏聚力為38.5 kPa，而運營時間增大至15、20和25 d后，黏聚力相比前者分別降低73.8%、75.1%和80%。從降低幅度變化可看出，在前5 d內，平均5 d黏聚力降低52.2%，而后續平均每5 d僅下降18.7%，表明運營時間愈長，黏聚力降低幅度減小，黏聚力的大幅降低僅發生在初始浸泡狀態下；內摩擦角在運營時間為0 d時為21.5°，而運營時間15、20和25 d時的內摩擦角相比前者分別降低13.9%、16.7%和18.5%，相比黏聚力特征參數，內摩擦角受運營時間影響敏感度顯著低于前者。筆者認為，在高水位浸泡下，初期由于滲流作用蔓延至堤防邊坡內部，滲流作用顯著較活躍，且與初始邊坡狀態具有顯著差異，局部土體孔隙水壓力顯著增大，因而造成土體黏聚力快速降低，而隨著運營時間延長，滲流作用趨于穩定，黏聚力受滲流作用影響逐漸減弱，且由于堤防內部存在一部分粉土等黏性土體，終究會具有部分黏性特征[15-16]，因而黏聚力呈先快后慢的降低態勢。相比于黏聚力，內摩擦角主要與顆粒咬合度有關，而滲流作用主要活躍于顆粒孔隙內部，對內摩擦角影響較小，因而高水位浸泡下內摩擦角平均每5d降低幅度僅為4%。

圖3 不同運營時間下土體特征參數與運營時間關系

3.2 堤防邊坡穩定性變化

圖4為堤防邊坡在不同運營時間下安全穩定系數變化曲線。從圖4中可看出，邊坡安全系數與運營時間為負相關關系，在運營時間0 d時的安全系數為5.94，而運營時間增大至20和25 d后，安全系數相比前者降低54.9%和60.4%。從降低速率來看，安全系數在前10 d內，平均每5 d安全系數降低27.3%；而在運營時間10～25 d時，安全系數平均每5 d降低8.7%。從安全系數與運營時間數學關系來看，兩者具有冪函數關系，此關系可為估算不同運營時間堤防工程安全穩定性提供參考。

圖4 不同運營時間下邊坡安全穩定系數變化曲線

圖5為不同運營時間下堤防邊坡內部危險滑動面分布狀態。從圖5中可看出，各運營時間下的危險滑動面均位于迎水側，且各危險滑動面的半徑隨運營時間顯著減小，運營時間5 d時危險滑動面半徑為3.833 m，而運營時間為20、25 d時滑動面半徑相比前者分別降低29.7%、37.9%，即運營時間愈長，愈加大了滑動面發生向外側的滑動趨勢。分析認為，當運營時間延長后，堤壩邊坡內部顆粒間有效應力顯著降低，對土顆粒骨架結構穩定性具有重要影響，減弱土顆粒間長期承載穩定性，進而增大邊坡滑動面滑移趨勢，從而產生邊坡危險滑動面半徑減小的現象。

圖5 不同運營時間下邊坡內部危險滑動面分布(橫坐標為堤防模型軸線長度、縱坐標為模型所處高度，零點為計算模型邊界點，下同)

4 高水位浸泡下堤防滲流場特征

4.1 多砂斷面

在高水位浸泡下，堤防滲流場安全穩定性亦會受到較大影響。為此，建立不同類型土質的堤防斷面有限元模型，其中多砂斷面堤防指堆筑材料砂土占比超過50%，但低于80%的堤防工程。按照網格劃分，共獲得20 167個滲流網格單元、17 658個節點，見圖6。

圖6 多砂斷面滲流場分析模型

經仿真計算，獲得多砂斷面堤防模型中流速矢量場以及浸潤線分布特征，見圖7。從圖7中流速場分布可知，不同土層中流速分布具有顯著差異，在第二層砂土層中存在隔水、阻水效應，該土層中無顯著流速矢量分布；而在首層粉土層中，流速分布較密集，滲流活躍性較大；由于砂土層在中間的隔水作用，造成底層細礫石層雖有流速分布，但分布狀態較平靜，無顯著滲流活躍性。分析指出，當堤防堆筑材料包括有除砂土層以外的其他土層時，內部滲流作用實質上仍處于較活躍狀態，特別是在高水位浸泡后。圖7中已標出浸潤線分布，其滲流逸出點位于背水面，表明滲流作用對背水面安全穩定性具有較大影響。

圖7 流速矢量場以及浸潤線分布特征(多砂斷面堤防模型)

經COMSOL計算出的滲流場特征還包括滲透坡降分布，見圖8。從圖8中可以看出，在首層土體背水面堤防肩部處所計算的滲透坡降為0.7，相比粉土等黏性土體的安全坡降允許值，首層土體在背水面極易發生滲透破壞，對于工程設計部門來說，應注重加固維護首層土體與中間土層在背水面的交界處。第三層礫石土層中計算獲得滲透坡降為0.05，迎水面的單寬滲流量為0.5 m3/d，從滲透坡降與單寬滲流量計算值可知，細礫石層中并不會出現管涌或流土等危險性現象。綜上分析，多砂斷面堤防工程應重點關注砂土層與粉土層交界面出現的滲透破壞。

圖8 滲透坡降分布(多砂斷面堤防模型)

4.2 全砂斷面

同理，建立全砂斷面分析模型。全砂斷面堤防指堆筑材料均為砂土，按照四面體單元劃分后獲得網格數2 458個，節點數1 962個，全砂端面堤防工程模型見圖9。

圖9 全砂斷面滲流場分析模型

圖10為全砂斷面堤防模型計算獲得流速矢量場分布及浸潤線特征。從圖10中流速分布可知，由頂部至底部，砂土層流速分布逐漸增強，特別是在底部流速矢量線分布較密集，滲流作用較活躍。分析認為，此與高水位浸泡下孔隙水壓力的發展方向有關，由堤防坡腳等部位逐漸蔓延至堤身，甚至背水側。與多砂斷面模型具有顯著差異的是，浸潤線由迎水面至背水側基本較平緩，滲流逸出點相比多砂斷面的高程更高，COMSOL所計算處的高程相比多砂斷面增大15%，表明全砂斷面堤防工程應注意水位上漲對堤防傾覆的威脅性。

圖10 流速矢量場以及浸潤線分布特征(全砂斷面堤防模型)

圖11為計算獲得的滲流場全砂斷面模型滲流坡降特征。滲透坡降反映了土體滲透破壞的危險性，從圖11中計算參數可知，全砂斷面中最大滲透坡降為0.65，位于堤防背水一側，且坡降較集中，雖相比多砂斷面有所降低，但仍高于堤防土體滲透坡降安全允許值，因而應考慮對堤防防滲加固。從堤防單寬滲流量亦可知，其最大滲流量可達3.58 m3/d，相比多砂斷面增長6.2倍，因而對于全砂斷面應加強對堤防滲漏等潰壩危險性方面監測。

圖11 滲透坡降分布(全砂斷面堤防模型)

5 結 論

利用COMSOL Multiphysical多物理場軟件，研究計算了高水位浸泡下堤防工程邊坡穩定性，主要分析了土體邊坡穩定性以及滲流場特征，主要結論如下：

1) 隨運營時間增大，邊坡土質黏聚力與內摩擦角均降低，但黏聚力受之影響更較顯著。黏聚力在前5d內，平均5d浸泡下黏聚力降低52.2%，而后續平均每5 d下降18.7%；內摩擦角在高水位浸泡下，平均每5 d降低幅度僅為4%。

2) 邊坡安全系數與運營時間為負相關關系，且兩者具有冪函數關系。前10 d運營時間內，平均每5 d安全系數降低27.3%，而在運營時間10～25 d時，安全系數平均每5 d降低8.7%；各危險滑動面的半徑隨運營時間延長而減小，加大了滑動面往外側滑移的趨勢，運營時間為20、25 d時，滑動面半徑相比運營時間5 d時分別降低29.7%、37.9%。

3) 多砂斷面堤防模型流速矢量場集中于粉土層與砂土層交界面，而浸潤線逸出點位于交界面背水側，堤防模型中最大滲透坡降為0.7，位于粉土層背水面肩部，易出現滲透破壞，迎水面的單寬滲流量為0.5 m3/d，細礫石土層中不會出現滲透破壞。

4) 全砂斷面堤防模型中，流速場由頂至底逐漸增強，底部滲透較活躍，滲流逸出點高程相比多砂斷面模型中增大15%；單寬滲流量相比多砂斷面模型增長6.2倍，且全砂斷面中最大滲透坡降為0.65，超出坡降安全允許值。