基于Abaqus的某堤防工程防滲加固結構水-力耦合場分析研究

溫青山

(水利部新疆維吾爾自治區水利水電勘測設計研究院，烏魯木齊 830001)

1 概 述

我國是遭受洪澇災害全球最嚴重國家之一[1]，建立堤防工程有助于鞏固人民生產生活安全，但不可忽視的是當前許多地區堤防工程修建質量較差，導致面臨較為嚴重防洪形勢[2]。水利工程技術人員為了提升堤防工程建設質量，采用過新型堤身材料加固防滲[3-4]；或采用人工灌漿固結措施，強化堤防工程滲流穩定性[5-6]；或基于水利-巖土工程交叉學科特點，引入巖土工程分析理論，多相耦合下研究堤防工程防滲加固[7-8]。本文針對堤防工程面臨水力與自身應力變形穩定性，利用數值軟件建立分析模型，引入多相耦合理論，探討堤防工程防滲加固結構滲流與應力特性，為評判或設計堤防工程防滲結構提供重要參考。

2 水-力耦合分析理論

2.1 滲流方程

水流在巖土層中運動稱之為滲流作用，發生滲流主要是由于巖土體中孔隙發育較多，同時水流中存在水頭壓差，造成水流有一定內在推動力[9]?；谒|點運動方程確定流體宏觀效應，獲得流體滲流模型，直觀反映滲流特征參數與介質實際狀態，并得到滲流場中介質流動規律[10]。

基于質量守恒定律與水質點連續運動假定，可得到滲流連續方程為[11-12]

(1)

式中：ρ為密度；vi為滲流速度；n為含水量；S為物理參數。

當巖土體處于非飽和狀態時，含水量n替換成θ，按照上式推導獲得非飽和滲流方程為：

(2)

進而聯立上述兩式，獲得飽和-非飽和連續性方程為：

(3)

代入達西公式，得：

(4)

為求解上述微分方程，需要給定初始邊界條件，以邊界條件作為方程解評判準則，按照Galerkin余量法與邊界條件，可得到矩陣表達式滲流方程為：

(5)

求解上述矩陣方程式，得：

[K]{hc}={P}

(6)

2.2 應力變形理論

基于材料本構模型方程，獲得巖土體應力變形解，依據摩爾庫倫破壞準則，鄧肯張模型可表述為：

(7)

式中：(σ1-σ3)f為峰值偏應力；c、φ為黏聚力與內摩擦角。

引入參量S根據材料強度殘余水平，衡量巖土體內是否存在極限平衡狀態，其表述為：

(8)

聯立上述兩式有：

(9)

水利工程中常常需要澆筑大面積混凝土，而以線彈性模型表述混凝土材料為：

{σ}=[D]{ε}

(10)

式中：[D]為彈性矩陣。

(11)

其中：d11=d22=d33=λ+2G；d21=d31=d32=λ；d44=d55=d66=G，λ，G為與巖土材料參數相關。

作為人工材料(混凝土)與巖土地基硬接觸，勢必會在接觸面產生一定滑移錯動，造成混凝土材料性能降低或巖土體產生次生裂縫。故而針對接觸面，選擇Goodman單元模型[13]，表述兩種材料間接觸面之間的應力位移關系為：

[σ]=[K0][ω]

(12)

式中：[σ]=[τyxτyyτyz]T為三方向應力矩陣；[W]=[ΔUΔVΔW]T為位移矩陣；[Ko]為本構模型參量。

2.3 滲流-應力耦合理論

借助應力剛度矩陣表達式， 平衡方程組可表述為[14]：

[K]{δ}={F}+{Fs}

(13)

式中：[K]為剛度矩陣；{δ}為位移矩陣；{F}為邊界約束矩陣；{Fs}為滲透體積力的等效荷載約束矩陣。

在多場耦合中，滲流場與應力場會改變彼此運動特征。例如，滲流場可能會造成應力場應力偏轉，或造成巖土材料孔隙率等參數發生變化；而應力場同樣可能會造成滲流方向在某一特定方向上重分布，故而確定耦合方程很有必要。

滲流場對應力場影響方程表述為：

(14)

式中：f為滲流體積力，fx、fy、fz為三方向分解力；H為水頭值；γw為重度；Jx、Jy、Jz為三方向上水力坡降。

而應力場對滲流場影響方程可表述為Carman-Kozeny模型，其表達式如下[15]：

(15)

式中：kz為Kozeny參數；Sp為表面積。

其中滲透系數與孔隙率存在如下關系式：

(16)

即可歸納出滲透系數與應力狀態存在間接函數關系式，表現應力場中滲流場分布狀態，如下式：

k=k(σij)

(17)

綜上所述，滲流場耦合應力場運動方程表述為：

[K]{δ}={F}+{Fs}+{F′}

[K′]{H}=[A]

(18)

利用有限元數值軟件，基于運動方程多次迭代計算，直至達到精度要求，獲得水-力耦合下滲流應力特征參數。

3 工程概況

某堤防工程為河流流經區域內重要的防洪設施，所在水庫正常蓄水位達45 m，乃是當地大型控制性水利樞紐工程，為該城市所轄5個縣區提供工農業用電，并調節農業用水資源。堤防工程沿江修建，研究區段長度約為3.5 km，該堤防工程所用材料為黏土與砂壤土材料，層次性顯著，堆筑厚度不高，密實性較差。由于土體屬中等滲透性土，含水量亦較高，達62%，防滲性較差，工程安全穩定性受到較大挑戰，故有必要對該堤防工程采取防滲加固、加高加寬坡頂等工程措施。

根據現場地質踏勘表明，該堤防工程地基上覆蓋較厚土層。研究土層性質得知，為第四系更新統堆積土，厚度約為4～6 m。另外，覆蓋層中還夾有部分防滲效果較好的黏土材料，部分地段出露有透水性較高的半風化砂礫石，易形成滲流通道，整個提防工程中砂土材料防滲效果較差，因而堤基礎防滲應考慮灌漿固化砂礫石層。地基層與上覆蓋土層之間還存在著厚3.1～4.2 m的砂壤土，屬晚更新統河流沖擊作用形成砂土，透水性較高，多層交織在基巖層與覆蓋層中。鉆孔資料表明，基巖層分別包括有泥質粉砂巖、灰巖與白云巖，樁號2+175～2+650段內基巖層為灰巖，塊狀構造，強度與完整性較高。但在樁號3+65～3+450段內包括有泥質粉砂巖，部分夾層內有砂礫石碎屑，且還發育巖溶，內充填材料含碳酸巖，連通性中等，對堤防地基防滲性是較大考驗。綜上所述地質分析，研究區段內堤防工程防滲性能受到巖溶與透水性強材料威脅，巖溶內材料在水-力作用下可能會造成滲透失穩破壞，影響堤防地基承載能力，故而考慮對堤防工程設置防滲墻，強化堤防滲流穩定性。

4 水-力耦合有限元分析

4.1 建立模型

針對該堤防工程滲流穩定性，采用防滲墻工程措施，滿足堤防穩定性與靜力穩定性要求。選取樁號1+125～1+150區段作為有限元分析計算模型，設計防滲墻厚度為0.5 m，深度達到基巖層，并穿過巖溶層，防滲墻下端采用固結灌漿方式，墻體材料抗壓強度均超過10 MPa，設計堤防迎水測水位為45.2 m，建立模型計算分析滲流場與應力變形特征。

利用Abaqus數值軟件按照選取區段剖面圖建立數值模型，設置堤防工程地基影響深度為75.8 m，堤防高度為6.9 m，簡化相鄰巖土層，采用耦合特性的六面體C3D10P微單元體，實現與堤防工程實際狀態相一致，并特別對防滲墻研究體單元網格加密，精度提高，整體模型共劃分出12 682個單元數，節點13 246個，見圖1。模型中各巖土層材料采用前文分析的本構模型，材料各物理參數均選取室內土工試驗報告，與邊界約束條件一起輸入數值軟件中。

圖1 數值模型

4.2 滲流場特征

圖2為計算出堤防工程中覆蓋層滲透系數在防滲墻施工過程中各工況下演化曲線。從圖2中可看出，在地應力平衡與溶洞固結灌漿工況下，覆蓋層土體滲透系數曲線幾乎為水平狀態；當上覆荷載增加，滲透系數大幅度降低，此主要是由于堆筑荷載增加，覆蓋層土體孔隙率降低，滲透系數故降低。在堤防堆筑工況滲透系數降低階段中，原狀土與填土滲透系數分別降低10.3%和11%；但細砂與圓礫石滲透系數分別降低3.2%和3.5%。更有甚者，砂礫石滲透系數幾乎僅僅只有一絲“抖動”變化。由此表明，砂礫石滲流特性幾乎不受填土堆筑影響，原狀土與填土滲流特性受堆筑工況影響最大。蓄水期，除填土外，其他土體材料滲透系數均呈降低態勢，受蓄水影響最大降低量為0.02；填土由于并未達到固結狀態，因而滲透系數呈上升態勢，提升了1.5%。

圖2 各工況下覆蓋層滲透系數演化曲線

圖3為未灌漿與灌漿工況下堤防中壓力水頭與總水頭等值線。從圖3可看出，在未灌漿工況下，迎水測壓力水頭隨水位上升，穩步上升至水平狀態；但在灌漿狀態下，總水頭線積聚在防滲墻周邊，表明附近區域水頭降低顯著，防滲效果較佳；灌漿導致防滲墻與之聯系成為防滲整體，滲透系數較低，與基巖相近，減少了溶洞與上覆蓋層滲漏量，總水頭下降0.33 m，浸潤面高程亦下降2.2 cm。

圖3 堤防中壓力水頭與總水頭等值線(從左至右依次為壓力水頭、總水頭)

計算出在溶洞灌漿與未灌漿狀態下防滲墻內最大水力坡降，見圖4。灌漿后防滲墻內最大水力坡降為30.14，且隨防滲墻高程增大仍維持在較高水平，高程40 m處相比高程10 m處水力坡降增大6.85倍。不可忽視的是，未灌漿狀態下隨高程增大，水力坡降穩定性愈差，在高程40～75 m范圍內，水力坡降最大變化幅度17%，不利于堤防工程防滲穩定性。

圖4 防滲墻內最大水力坡降

表1為未灌漿與灌漿狀態下各滲流特征參數。從表1中可看出，未灌漿狀態下堤防總滲漏量為4.56×10-4m3/s，當對溶洞灌漿與防滲墻成為整體防滲結構，總滲漏量減少3.18×10-4m3/s，僅為1.38×10-4m3/s。由于溶洞內灌漿，抑制滲流作用，溶洞內最大滲流速度僅為1.98×10-7m/s，相比未灌漿狀態下減少一個量級之多。但在堤防自由面上最大滲流速度灌漿與否變化不大，兩種狀態下最大滲流速度均接近7×10-6m/s，表明灌漿措施并不改變自由面上滲流狀態，但改變了溶洞內滲流狀態，引起總滲漏量降低。

表1 未灌漿與灌漿狀態下各滲流特征參數

4.3 應力變形特性

4.3.1 應力特征

圖5為施工期防滲墻應力分布云圖，從最大最小主應力特征來看，未灌漿狀態下最大拉應力為0.361 MPa，最大壓應力為1.04 MPa，拉應力隨墻內深度增加而增大，不論是灌漿亦或不灌漿，壓應力區域集中在覆蓋層與防滲墻接觸區，高程在墻體35 m上部。對溶洞灌漿后，防滲墻最大拉、壓應力分別為0.348和0.77 MPa，集中在覆蓋層與基巖層之間接觸面處。

圖5 施工期防滲墻應力分布云圖(依次為拉應力、壓應力)

圖6 蓄水期防滲墻應力分布云圖

從蓄水期防滲墻應力云圖可看出(圖6)，蓄水期迎水側拉、壓應力均顯著增長，并出現有受拉區與受壓區，未灌漿狀態下最大拉應力為4.56 MPa，最大壓應力為6.23 MPa?；鶐r與覆蓋層中應力分布呈突變性，分析是由于基巖層彈性模量與覆蓋層彈性模量相差較大，造成應力值出現急劇折變現象。灌漿后，最大拉應力與最大壓應力均有一定程度增大，灌漿后基巖透水性得到抑制，造成迎水側水頭差相應增大，故引起拉應力亦增大。綜上應力分析可知，灌漿對堤防滲漏性可極大抑制，但一定程度亦會加劇防滲墻材料應力集中破壞性。因此，防滲墻應挑選柔性材料，緩和較強拉應力。

4.3.2 變形特征

圖7為防滲墻豎向位移(沉降)與高程關系圖。從圖7中可看出，各工況下豎向位移整體均隨防滲墻高度增加而逐漸增大，最大沉降量均位于墻頂部，但各工況中以未灌漿滿蓄水位為最大沉降量，達2.49 mm，灌漿防滲墻施工期最大沉降位移最小。堤防工程在蓄水后防滲墻頂面臨較大的水頭壓力，造成滲流場與應力場多重沖擊下，沉降位移故而亦是最大。

圖7 防滲墻豎向沉降與高程關系圖

圖8為施工期與蓄水期防滲墻水平向位移與墻高程關系，施工期最大水平位移出現在墻中部，墻底部與頂部水平位移較低，施工期模擬計算時已假定防滲墻為剛性結構水平位移約束較強，無自由度變形。灌漿堆筑期與未灌漿填筑期最大水平位移均出現在高程38 m處，分別為0.59和0.77 mm，堤防填筑勢必會加劇防滲墻背水測的彎矩變形，引起局部拉應力集中，故而填筑期再防滲墻中上部地帶出現最大位移。滿蓄期不論是灌漿還是未灌漿，水平向位移曲線均呈一致性，蓄水期迎水側面臨較大水頭壓力，此時防滲墻可視為剛性構件，其水平位移受溶洞灌漿與否影響不大；相比施工期，蓄水期受水頭壓差影響，防滲墻勢必水平位移會增大，最大水平位移達6.8 mm，高出施工期最大值7.8倍。

圖8 防滲墻水平位移-墻高程關系

5 結 論

針對某堤防工程防滲加固結構，引入水-力耦合分析理論，基于Abaqus建立堤防工程防滲加固結構模型，分析滲流場與應力場特征，得到以下結論：

1) 覆蓋層滲透系數受堤防堆筑與蓄水影響較大，隨堆筑荷載增加與蓄水位增大，滲透系數逐漸降低，填土埋深最淺，蓄水導致滲透系數降低0.02；灌漿與防滲墻整體存在可強化堤防結構防滲，降低總滲漏量3.18×10-4m3/s，降低總水頭0.33 m，浸潤面下降2.2 cm。

2) 施工期、蓄水期防滲墻應力分布表明，施工期壓應力區集中在覆蓋層與防滲墻接觸面，墻體高程35 m處；蓄水期拉、壓應力增大，未灌漿下最大拉、壓應力分別為4.56和6.23 MPa。

3) 防滲墻豎向沉降、水平位移與墻高程關系顯示，豎向沉降隨墻高度增加而增大，最大沉降位于墻頂部，未灌漿滿蓄期沉降最大，達2.49 mm；施工期最大水平位移出現在墻中部，墻底部與頂部水平位移較低，接近零，蓄水期灌漿與未灌漿位移曲線一致性，最大水平位移達6.8 mm，高出施工期最大值7.8倍。