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

林俊屹

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

1 概 述

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

2 工程介紹及模型建立

2.1 工程概況

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

2.2 模型建立

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

圖1 堤防斷面圖

圖2 數(shù)值模型及網(wǎng)格模型圖

3 高水位浸泡下堤防邊坡穩(wěn)定性

3.1 物理力學特征參數(shù)變化

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

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

3.2 堤防邊坡穩(wěn)定性變化

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

圖4 不同運營時間下邊坡安全穩(wěn)定系數(shù)變化曲線

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

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

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

4.1 多砂斷面

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

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

經(jīng)仿真計算，獲得多砂斷面堤防模型中流速矢量場以及浸潤線分布特征，見圖7。從圖7中流速場分布可知，不同土層中流速分布具有顯著差異，在第二層砂土層中存在隔水、阻水效應，該土層中無顯著流速矢量分布；而在首層粉土層中，流速分布較密集，滲流活躍性較大；由于砂土層在中間的隔水作用，造成底層細礫石層雖有流速分布，但分布狀態(tài)較平靜，無顯著滲流活躍性。分析指出，當?shù)谭蓝阎牧习ㄓ谐巴翆右酝獾钠渌翆訒r，內(nèi)部滲流作用實質(zhì)上仍處于較活躍狀態(tài)，特別是在高水位浸泡后。圖7中已標出浸潤線分布，其滲流逸出點位于背水面，表明滲流作用對背水面安全穩(wěn)定性具有較大影響。

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

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

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

4.2 全砂斷面

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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