董志勇,龐鑫杰,周建芬
(1.浙江工業大學土木工程學院,浙江 杭州 310023;2.浙江水利水電學院水利與環境工程學院,浙江 杭州 310018)
蘭溪古城墻沿蘭江修建,始建于北宋,經歷了多次修復和重建,于2017年被列為浙江省省級文物保護對象。該城墻擔負著老城區的防洪安全重任,同時也是當地重要的人文景觀。在梅汛期和臺汛期,由于連綿陰雨或大暴雨的作用,古城墻存在安全穩定性問題。如2017年洪水期間,古城墻出現多處冒水、管涌等現象,且部分坍塌,洪水灌入蘭溪老城,導致城區被淹,約10萬人受災。
目前,計算流體力學和計算固體力學的快速發展實現了滲流場與應力場耦合。土體滲流與應力耦合理論最早是根據多孔介質的力學特征描述的[1-2]。太沙基(Terzaghi)和比奧(Biot)是多孔介質中耦合問題研究的理論奠基者。20世紀初,太沙基提出了有效應力原理、一維固結理論,并進一步擴展到三維,建立了太沙基-倫杜立克方程,為土力學中滲流場與應力場的耦合奠定了理論基礎。比奧推導出三維流固耦合方程,建立了土體骨架變形與孔隙水壓力之間的關系,并在求解過程中直接體現了滲流場與應力場的相互作用[3]。但是比奧固結理論設計參數較多,難以準確測得,且得出方程的解析解也較困難。因此,在一段時間內比奧固結方程沒有在工程中得到廣泛應用。隨著現代計算技術的發展,尤其是有限元法的發展,比奧固結理論才得以應用于工程實踐中。國外,Noorishad等[4-5]提出了裂隙滲流與應力的耦合分析模型;Asaolollahi等[6]在凍土地區隧道涌水問題中,考慮了滲流場、應力場以及溫度場三場之間的相互作用,并探究了在三場耦合作用下隧道的穩定性問題。國內,平揚等[7]以比奧固結理論為基礎,研究在開挖過程中滲流場與應力場的變化規律及其對基坑穩定的影響。王媛[8]根據比奧固結理論、伽遼金有限元法,聯立了滲流與應力方程,建立了以位移和孔隙水壓力為未知量的方程,可以同時計算出滲流場與應力場,并以重力壩為例驗證了該方法和計算程序的可行性。柴軍瑞等[9]從均質土壩的滲透特性出發,分析了均質土壩滲流場與應力場相互作用的力學機制,提出了均質土壩滲流場與應力場間接耦合分析的連續介質數學模型,并討論了該數學模型的有限元數值解法。
國內外學者的廣泛關注和研究,推動了滲流-應力耦合理論的應用和發展。與此同時,邊坡滲流穩定分析也得到了迅速發展。毛昶熙等[10]用土體單元所受的滲透力替代土條周邊的孔隙水壓力,從而開辟了新的分析土坡穩定性的途徑。這種方法考慮了滲流場與應力場的相互作用,從而提高了穩定性分析的精度,并在幾個實際工程中得到了理想的結果。隨后,黃俊等[11]提出用土條所受的滲透力代替孔隙水壓力的邊坡穩定計算方法,在實際工程中得到了較好的結果。汪自力等[12]在有限元計算基礎上,也用土體單元所受的滲透力代替其周邊的孔隙水壓力,采用數學規劃中的單純形法自動尋找最小安全系數;李吉慶等[13]以極限平衡理論為基礎,對比分析了畢肖普條分法和有限元單元法,指出了條分法在計算滲流作用下邊坡安全系數時存在的缺點,并論證了以土條受到的滲透力來代替土條邊界水壓力的有限元單元法的優越性,并通過模型土壩及幾個實例土壩加以驗證;徐衛亞等[14]綜述了國內外降雨型堆積體滑坡滲流穩定性研究進展,并對今后研究作了展望;劉永濤等[15]結合EEMD和RVM算法的優點,構建了土石壩滲流量時間序列的預測模型,為土石壩滲流分析提供了新的方法。在研究非飽和土方面,張信貴等[16]根據飽和與非飽和土體抗剪強度及滲流特性,推導了飽和與非飽和土體在應力場和滲流場耦合作用下的有限元方程。韋立德等[17]建立一個考慮飽和與非飽和滲流場與應力場耦合的三維強度折減有限元模型。張玉軍等[18-20]通過推導非飽和土中滲流與應力耦合有關的方程獲得解析解。田東方等[21]闡述了坡面徑流-非飽和滲流場與應力場耦合的重要意義和作用,提出了滲流場與應力場間接耦合的計算方法。
隨著有限元商業軟件的發展,許多學者利用ABAQUS、FLAC3D、GEO-SLOPE、PLAXIS、北京理正等研究邊坡穩定性問題。徐海奔[22]利用ABAQUS有限元軟件對飽和與非飽和滲流規律進行了初步探討。唐曉松等[23-25]通過PLAXIS、ADINA等有限元軟件對滲流作用下的邊坡穩定性進行了分析,并利用GEO-SLOPE程序進行驗算;Desai等[26-27]利用有限元軟件模擬庫水位驟降下邊坡的穩定性。綜上,雖然國內外學者對滲流場與應力場耦合分析研究較多,但在實際加固設計中考慮較少,尤其是耦合作用對擋土墻穩定性的影響。本文在降雨入滲條件下,基于ABAQUS軟件,對蘭溪市某防洪古城墻邊坡措施的滲流場和應力場進行數值模擬,比較分析了不同加固措施對古城墻穩定性的影響。
邊坡內滲流場與應力場的相互作用、相互影響較為復雜。一方面,滲流場會因外界條件的改變而變化,致使土體孔隙水壓力分布發生變化,從而導致土體和作用在土體上的有效應力與滲透力發生變化,即改變了應力場;另一方面,應力場的變化會導致土體發生變形,從而改變土體的孔隙率和滲透系數,即改變了滲流場。
由于巖土體的孔隙結構比較復雜,通常把流體和固體看作一個理想的連續體,用這個連續體代替復雜的孔隙結構,建立滲流與應力耦合的控制方程。ABAQUS軟件控制方程就是基于連續介質的概念將多孔介質分為固相和液相。固相應力平衡可用虛功原理表示:
(1)
式中:σ為有效應力;δε為虛變形速率;δv為虛速度場;T為單位面積的表面力;f為體積力(不包括流體質量);s為固相的飽和度;n為固相材料的孔隙率;ρw為流體密度;g為重力加速度。
應用有限元法的位移模式和拉格朗日公式可以將虛功方程離散化,得到固相的有限元網格,并使液相通過網格。流體必須滿足連續性方程,使得某時段流入或流出的流量等于流體體積的變化率,滲流連續性方程為
(2)
式中:ρe為流體的參考密度;vf為滲流速度;n為S面外法線方向。
式(2)應用了反向歐拉法得到近似積分,將孔隙水壓力看作一個變量,并對其進行有限元離散。另外,孔隙內流體遵循2個定律,分別在不同條件下采用:一是Darcy定律,用于流速較低的情況;二是Forchheimer定律,用于流速較高的情況。
降雨入滲隨時間和空間變化,通常用土體允許入滲容量fp、土體飽和滲透系數ks和降水強度q這3個因素之間的關系描述降雨入滲模型[28]。3種不同情況如下:①當q≤ks時,降雨全部滲入土壤,地表不發生徑流;②當fp>q>ks時,雨水全部入滲,fp隨著入滲深度的增加呈遞減趨勢,但這個階段不會發生徑流;③當q≥fp時,雨水不能全部入滲到土壤中,所以會產生徑流。
本文在模擬整個降雨過程中,默認入滲率大于降水強度,將降水強度設置為流量邊界,從而模擬入滲問題。
1.3.1接觸模擬
采用ABAQUS軟件中的相互作用模塊,在該模塊中設置擋土墻與土體的接觸面,包含接觸面的法向作用和切向作用。其中,接觸面的法向模型采用硬接觸,即兩物體只有在相互壓緊的狀態下才能傳遞法向壓力,若兩物體之間有間隙則不能傳遞法向壓力;接觸面的切向作用考慮為無摩擦,即默認古城墻的墻背光滑。
1.3.2滲流模擬
當接觸界面兩側都存在孔壓自由度時,則計算的接觸壓力為有效壓力,即不考慮孔隙水壓力的影響。如果接觸界面只有一側包含孔壓自由度,則不會出現流體流入或穿過接觸界面的現象,在這種情況下,接觸壓力是總壓力,即考慮孔隙水壓力的影響。
根據《水工擋土墻設計規范》,抗傾覆穩定安全系數K的計算公式為
(3)
式中∑Mv、∑MH分別為對擋土墻基底前趾的抗傾覆力矩和傾覆力矩。
蘭溪市某古城墻是蘭江防洪堤的一部分,其斷面為梯形,高6 m,頂寬0.6 m,底寬1.2 m,墻后土體為砂性土,地基為不透水巖基,橫斷面如圖1所示。墻后土體的材料屬性:干密度為1 400 kg/m3,黏聚力為5 kPa,有效內摩角為25 °,彈性模量為100 MPa,泊松比為0.32,飽和滲透系數為1.354 m/h,初始孔隙比為1.0。古城墻體的材料屬性:密度為2 400 kg/m3,彈性模量為2 000 MPa,黏聚力為260 kPa,有效內摩角為51°,泊松比為0.1,飽和滲透系數為0.02 m/d。由于地基為不透水巖基,因此不考慮地基對古城墻的影響。邊界條件設置:模型右側邊界施加水平方向位移約束,邊坡底部為全約束邊界和不透水邊界。接觸設置:古城墻與土體采用面面接觸,接觸面法向模擬采用硬接觸,并設置摩擦系數為0.3。古城墻體和土體采用本構模型和屈服準則,分別為理想彈塑性本構模型和Mohr-Coulomb強度準則;選用CPE6M六結點修正二次平面應變三角形單元,計算網格布置如圖2所示。

圖1 古城墻橫斷面(單位:m)

圖2 有限元模型的計算網格
圖3中黑色實心圓表示由土壓力計測得的離地面4.5 m、離古城墻墻頂3 m處土壓力。由圖3可知,ABAQUS模擬值在深度小于1 m時土壓力為零,在深度大于1 m時土壓力隨著深度的增加而增大,呈三角形分布,符合土力學理論。土壓力計的實測值與ABAQUS模擬值接近,而庫倫理論值偏大,這主要因為理論計算沒有考慮黏聚力對土的影響,且忽略了墻體與填土之間相互作用以及位移協調作用。因此,基于ABAQUS軟件對降雨入滲條件下古城墻墻后土體孔隙水壓力和墻體穩定性進行模擬是可行的。

圖3 水平應力隨深度的變化
在邊坡除險加固中,用倒掛井作為防滲處理設施,可減少大量的土石方開挖,且防滲效果好。本文對古城墻采用倒掛井聯合支護的方式,截斷堤身的滲流通道,并探究古城墻孔隙水壓力變化過程。有防滲設施的古城墻橫斷面如圖4所示,倒掛井直徑為1.6 m,深度約為6 m,倒掛井井壁材料為拱頂及拱座倒掛筋,井內回填C25W6F50砼,古城墻與倒掛井通過每間隔6 m設置一道混凝土拉梁進行連接。倒掛井的材料屬性:密度為2 400 kg/m3,彈性模量為20 GPa,泊松比為0.25,飽和滲透系數為0.01 m/d。數值模擬中需計算模型所有節點的應力、位移、孔隙水壓力、流速以及接觸應力。

圖4 有倒掛井的古城墻橫斷面(單位:m)
為了探究有倒掛井的古城墻孔隙水壓力變化過程,采用瞬態方式進行降雨入滲全過程的滲流場與應力場耦合分析,邊坡頂部為流量邊界,雨水入滲強度為1 354 mm/h,歷時24 h,倒掛井與土體、土體與古城墻的連接選用面面接觸。有倒掛井的古城墻孔隙水壓力遭遇降雨歷時的等值線如圖5所示。為了能準確定位浸潤線,將孔隙水壓力為負值的區域設置為灰色,即藍色等值線以上的區域,則孔隙水壓力為0的等值線就是浸潤線。與此同時,也能更好地區分飽和區與非飽和區。由圖5可知,隨著降雨歷時的增加,孔隙水壓力不斷增大,浸潤線逐漸升高,古城墻后土體的飽和區域不斷擴大,非飽和區域不斷縮小,直至第5小時,土體達到飽和狀態,雨水不能入滲到土體中,以致發生地表徑流。

圖5 有倒掛井的古城墻墻背孔隙水壓力等值線
反濾層常用于壩基、壩體、重力式碼頭和板樁碼頭中,可以防止土顆粒隨滲流析出,并且能保證滲流順暢,以達到排水減壓的效果。本節采用反濾層加排水管以降低孔隙水壓力,并分析在有、無排水設施情況下,古城墻孔隙水壓力和流速的變化過程。具有排水設施的古城墻如圖6所示,反濾層從內向外由直徑0.5~1.5 mm的砂粒、直徑5~25 mm的碎石和直徑25~100 mm的碎石組成,滲透系數從內向外分別為0.05 m/d、0.3 m/d和0.8 m/d,排水管設置在距離底部1 m處,直徑為50 mm。

圖6 有排水設施的古城墻橫斷面(單位:m)
為了探究有排水設施的古城墻孔隙水壓力變化過程,采用瞬態方式進行降雨入滲全過程的滲流場與應力場耦合模擬,邊坡頂部為流量邊界,雨水入滲強度為1 354 mm/h,歷時24 h,計算中土體與古城墻的連接選用面面接觸。由圖7可見,在墻后水位到達排水管處之前,孔隙水壓力增長趨勢較快,當水位超過排水管時,孔隙水壓力的增長速率和浸潤線升高的趨勢明顯變緩,墻后水流不斷被排出。

圖7 有排水設施的古城墻墻背孔隙水壓力等值線
為了更好地分析排水設施對古城墻后滲流場的影響,在相同計算條件下,對無排水設施的古城墻進行了降雨滲流分析,并將2種情況的計算結果作了比較。對比圖7(d)和圖8可知,由于沒有排水設施,墻后的地下水無法排出,在經歷相同的降雨歷時后,無排水設施的古城墻土體的孔隙水壓力遠大于具有排水設施情形。

圖8 持續降雨24 h后無排水設施古城墻孔隙水壓力
由圖9可見,降雨滲流可通過排水管排出,從而形成滲流通道,而無排水設施的古城墻因土體已處于飽和狀態以致于后續的降雨滲流溢出。

圖9 降雨24 h后古城墻流速矢量比較
在有排水設施、無排水設施和倒掛井3種情形下,降雨24 h后的古城墻接觸應力與深度的關系如圖10所示。從圖10中可以看出,不同情形接觸應力隨深度變化的規律基本一致,隨深度的增加而增大,近似呈三角形分布,但三者在數值上存在差異,尤其與有排水設施差異較大。有排水設施的古城墻后的地下水位比其他2種情形都低,其對古城墻的應力也都小于另外2種情形。另外,倒掛井是采用“堵”水的方式,在計算結果上與無排水設施情形相近。
由表1可知,由于排水設施的存在,墻后孔隙水壓力對古城墻的壓力大幅度減小,相比其他2種情形,有排水設施古城墻受到的總壓力和傾覆力矩最小,也是最穩定和最安全的。有倒掛井情形,墻背所受的壓力和傾覆力矩都有所減小,但是效果并不明顯,因此,設倒掛井并不能有效提高古城墻的穩定性。

表1 降雨24 h后古城墻抗傾穩定性
由圖11可見,2 h時有、無排水設施B點水平位移增長的量級是相同的,是線性增長,但有倒掛井的增長緩慢,主要是因為倒掛井的消減土壓力作用。無排水設施和有倒掛井情形古城墻B點(圖1)水平位移先隨時間增加而增大,而后趨于平緩。這是由于在降雨5 h后,墻后地下水位已上升至墻頂,墻后土體處于飽和狀態,其孔隙水壓力和土壓力將一直保持不變。而具有排水設施的墻后地下水在降雨2 h后位于排水管處,因此可看出有排水設施城墻B點位移起初呈快速增長的趨勢,但在2 h后變得緩慢。比較3種情形的B點水平位移,由于無排水設施和有倒掛井情形無排水通道,地下水無法及時排出,墻體所受壓力較大,因此B點水平位移均比有排水設施情形大。

圖11 不同情形古城墻B點水平位移隨時間的變化
古城墻無排水設施時,歷經24 h降雨后,墻后孔隙水壓力增大,致使墻背所受壓力增大,降低了古城墻的穩定性;有排水設施時,墻后孔隙水壓力大幅度減小,增強了古城墻的穩定性;有防滲設施時,墻背所受壓力并不能有效降低,安全系數未明顯增大。