董志勇，龐鑫杰，周建芬

(1.浙江工業(yè)大學(xué)土木工程學(xué)院，浙江 杭州 310023；2.浙江水利水電學(xué)院水利與環(huán)境工程學(xué)院，浙江 杭州 310018)

蘭溪古城墻沿蘭江修建，始建于北宋，經(jīng)歷了多次修復(fù)和重建，于2017年被列為浙江省省級文物保護(hù)對象。該城墻擔(dān)負(fù)著老城區(qū)的防洪安全重任，同時也是當(dāng)?shù)刂匾娜宋木坝^。在梅汛期和臺汛期，由于連綿陰雨或大暴雨的作用，古城墻存在安全穩(wěn)定性問題。如2017年洪水期間，古城墻出現(xiàn)多處冒水、管涌等現(xiàn)象，且部分坍塌，洪水灌入蘭溪老城，導(dǎo)致城區(qū)被淹，約10萬人受災(zāi)。

目前，計算流體力學(xué)和計算固體力學(xué)的快速發(fā)展實(shí)現(xiàn)了滲流場與應(yīng)力場耦合。土體滲流與應(yīng)力耦合理論最早是根據(jù)多孔介質(zhì)的力學(xué)特征描述的[1-2]。太沙基(Terzaghi)和比奧(Biot)是多孔介質(zhì)中耦合問題研究的理論奠基者。20世紀(jì)初，太沙基提出了有效應(yīng)力原理、一維固結(jié)理論，并進(jìn)一步擴(kuò)展到三維，建立了太沙基-倫杜立克方程，為土力學(xué)中滲流場與應(yīng)力場的耦合奠定了理論基礎(chǔ)。比奧推導(dǎo)出三維流固耦合方程，建立了土體骨架變形與孔隙水壓力之間的關(guān)系，并在求解過程中直接體現(xiàn)了滲流場與應(yīng)力場的相互作用[3]。但是比奧固結(jié)理論設(shè)計參數(shù)較多，難以準(zhǔn)確測得，且得出方程的解析解也較困難。因此，在一段時間內(nèi)比奧固結(jié)方程沒有在工程中得到廣泛應(yīng)用。隨著現(xiàn)代計算技術(shù)的發(fā)展，尤其是有限元法的發(fā)展，比奧固結(jié)理論才得以應(yīng)用于工程實(shí)踐中。國外，Noorishad等[4-5]提出了裂隙滲流與應(yīng)力的耦合分析模型；Asaolollahi等[6]在凍土地區(qū)隧道涌水問題中，考慮了滲流場、應(yīng)力場以及溫度場三場之間的相互作用，并探究了在三場耦合作用下隧道的穩(wěn)定性問題。國內(nèi)，平揚(yáng)等[7]以比奧固結(jié)理論為基礎(chǔ)，研究在開挖過程中滲流場與應(yīng)力場的變化規(guī)律及其對基坑穩(wěn)定的影響。王媛[8]根據(jù)比奧固結(jié)理論、伽遼金有限元法，聯(lián)立了滲流與應(yīng)力方程，建立了以位移和孔隙水壓力為未知量的方程，可以同時計算出滲流場與應(yīng)力場，并以重力壩為例驗(yàn)證了該方法和計算程序的可行性。柴軍瑞等[9]從均質(zhì)土壩的滲透特性出發(fā)，分析了均質(zhì)土壩滲流場與應(yīng)力場相互作用的力學(xué)機(jī)制，提出了均質(zhì)土壩滲流場與應(yīng)力場間接耦合分析的連續(xù)介質(zhì)數(shù)學(xué)模型，并討論了該數(shù)學(xué)模型的有限元數(shù)值解法。

國內(nèi)外學(xué)者的廣泛關(guān)注和研究，推動了滲流-應(yīng)力耦合理論的應(yīng)用和發(fā)展。與此同時，邊坡滲流穩(wěn)定分析也得到了迅速發(fā)展。毛昶熙等[10]用土體單元所受的滲透力替代土條周邊的孔隙水壓力，從而開辟了新的分析土坡穩(wěn)定性的途徑。這種方法考慮了滲流場與應(yīng)力場的相互作用，從而提高了穩(wěn)定性分析的精度，并在幾個實(shí)際工程中得到了理想的結(jié)果。隨后，黃俊等[11]提出用土條所受的滲透力代替孔隙水壓力的邊坡穩(wěn)定計算方法，在實(shí)際工程中得到了較好的結(jié)果。汪自力等[12]在有限元計算基礎(chǔ)上，也用土體單元所受的滲透力代替其周邊的孔隙水壓力，采用數(shù)學(xué)規(guī)劃中的單純形法自動尋找最小安全系數(shù)；李吉慶等[13]以極限平衡理論為基礎(chǔ)，對比分析了畢肖普條分法和有限元單元法，指出了條分法在計算滲流作用下邊坡安全系數(shù)時存在的缺點(diǎn)，并論證了以土條受到的滲透力來代替土條邊界水壓力的有限元單元法的優(yōu)越性，并通過模型土壩及幾個實(shí)例土壩加以驗(yàn)證；徐衛(wèi)亞等[14]綜述了國內(nèi)外降雨型堆積體滑坡滲流穩(wěn)定性研究進(jìn)展，并對今后研究作了展望；劉永濤等[15]結(jié)合EEMD和RVM算法的優(yōu)點(diǎn)，構(gòu)建了土石壩滲流量時間序列的預(yù)測模型，為土石壩滲流分析提供了新的方法。在研究非飽和土方面，張信貴等[16]根據(jù)飽和與非飽和土體抗剪強(qiáng)度及滲流特性，推導(dǎo)了飽和與非飽和土體在應(yīng)力場和滲流場耦合作用下的有限元方程。韋立德等[17]建立一個考慮飽和與非飽和滲流場與應(yīng)力場耦合的三維強(qiáng)度折減有限元模型。張玉軍等[18-20]通過推導(dǎo)非飽和土中滲流與應(yīng)力耦合有關(guān)的方程獲得解析解。田東方等[21]闡述了坡面徑流-非飽和滲流場與應(yīng)力場耦合的重要意義和作用，提出了滲流場與應(yīng)力場間接耦合的計算方法。

隨著有限元商業(yè)軟件的發(fā)展，許多學(xué)者利用ABAQUS、FLAC3D、GEO-SLOPE、PLAXIS、北京理正等研究邊坡穩(wěn)定性問題。徐海奔[22]利用ABAQUS有限元軟件對飽和與非飽和滲流規(guī)律進(jìn)行了初步探討。唐曉松等[23-25]通過PLAXIS、ADINA等有限元軟件對滲流作用下的邊坡穩(wěn)定性進(jìn)行了分析，并利用GEO-SLOPE程序進(jìn)行驗(yàn)算；Desai等[26-27]利用有限元軟件模擬庫水位驟降下邊坡的穩(wěn)定性。綜上，雖然國內(nèi)外學(xué)者對滲流場與應(yīng)力場耦合分析研究較多，但在實(shí)際加固設(shè)計中考慮較少，尤其是耦合作用對擋土墻穩(wěn)定性的影響。本文在降雨入滲條件下，基于ABAQUS軟件，對蘭溪市某防洪古城墻邊坡措施的滲流場和應(yīng)力場進(jìn)行數(shù)值模擬，比較分析了不同加固措施對古城墻穩(wěn)定性的影響。

1 數(shù)值模擬方法

邊坡內(nèi)滲流場與應(yīng)力場的相互作用、相互影響較為復(fù)雜。一方面，滲流場會因外界條件的改變而變化，致使土體孔隙水壓力分布發(fā)生變化，從而導(dǎo)致土體和作用在土體上的有效應(yīng)力與滲透力發(fā)生變化，即改變了應(yīng)力場；另一方面，應(yīng)力場的變化會導(dǎo)致土體發(fā)生變形，從而改變土體的孔隙率和滲透系數(shù)，即改變了滲流場。

1.1 流固耦合分析的基本原理

由于巖土體的孔隙結(jié)構(gòu)比較復(fù)雜，通常把流體和固體看作一個理想的連續(xù)體，用這個連續(xù)體代替復(fù)雜的孔隙結(jié)構(gòu)，建立滲流與應(yīng)力耦合的控制方程。ABAQUS軟件控制方程就是基于連續(xù)介質(zhì)的概念將多孔介質(zhì)分為固相和液相。固相應(yīng)力平衡可用虛功原理表示：

(1)

式中：σ為有效應(yīng)力；δε為虛變形速率；δv為虛速度場；T為單位面積的表面力；f為體積力(不包括流體質(zhì)量)；s為固相的飽和度；n為固相材料的孔隙率；ρw為流體密度；g為重力加速度。

應(yīng)用有限元法的位移模式和拉格朗日公式可以將虛功方程離散化，得到固相的有限元網(wǎng)格，并使液相通過網(wǎng)格。流體必須滿足連續(xù)性方程，使得某時段流入或流出的流量等于流體體積的變化率，滲流連續(xù)性方程為

(2)

式中：ρe為流體的參考密度；vf為滲流速度；n為S面外法線方向。

式(2)應(yīng)用了反向歐拉法得到近似積分，將孔隙水壓力看作一個變量，并對其進(jìn)行有限元離散。另外，孔隙內(nèi)流體遵循2個定律，分別在不同條件下采用：一是Darcy定律，用于流速較低的情況；二是Forchheimer定律，用于流速較高的情況。

1.2 降雨入滲邊界條件處理

降雨入滲隨時間和空間變化，通常用土體允許入滲容量fp、土體飽和滲透系數(shù)ks和降水強(qiáng)度q這3個因素之間的關(guān)系描述降雨入滲模型[28]。3種不同情況如下：①當(dāng)q≤ks時，降雨全部滲入土壤，地表不發(fā)生徑流；②當(dāng)fp>q>ks時，雨水全部入滲，fp隨著入滲深度的增加呈遞減趨勢，但這個階段不會發(fā)生徑流；③當(dāng)q≥fp時，雨水不能全部入滲到土壤中，所以會產(chǎn)生徑流。

本文在模擬整個降雨過程中，默認(rèn)入滲率大于降水強(qiáng)度，將降水強(qiáng)度設(shè)置為流量邊界，從而模擬入滲問題。

1.3 墻體與土體接觸和滲流計算

1.3.1接觸模擬

采用ABAQUS軟件中的相互作用模塊，在該模塊中設(shè)置擋土墻與土體的接觸面，包含接觸面的法向作用和切向作用。其中，接觸面的法向模型采用硬接觸，即兩物體只有在相互壓緊的狀態(tài)下才能傳遞法向壓力，若兩物體之間有間隙則不能傳遞法向壓力；接觸面的切向作用考慮為無摩擦，即默認(rèn)古城墻的墻背光滑。

1.3.2滲流模擬

當(dāng)接觸界面兩側(cè)都存在孔壓自由度時，則計算的接觸壓力為有效壓力，即不考慮孔隙水壓力的影響。如果接觸界面只有一側(cè)包含孔壓自由度，則不會出現(xiàn)流體流入或穿過接觸界面的現(xiàn)象，在這種情況下，接觸壓力是總壓力，即考慮孔隙水壓力的影響。

1.4 擋土墻穩(wěn)定性計算

根據(jù)《水工擋土墻設(shè)計規(guī)范》，抗傾覆穩(wěn)定安全系數(shù)K的計算公式為

(3)

式中∑Mv、∑MH分別為對擋土墻基底前趾的抗傾覆力矩和傾覆力矩。

2 古城墻滲流分析

2.1 古城墻概況

蘭溪市某古城墻是蘭江防洪堤的一部分，其斷面為梯形，高6 m，頂寬0.6 m，底寬1.2 m，墻后土體為砂性土，地基為不透水巖基，橫斷面如圖1所示。墻后土體的材料屬性：干密度為1 400 kg/m3，黏聚力為5 kPa，有效內(nèi)摩角為25 °，彈性模量為100 MPa，泊松比為0.32，飽和滲透系數(shù)為1.354 m/h，初始孔隙比為1.0。古城墻體的材料屬性：密度為2 400 kg/m3，彈性模量為2 000 MPa，黏聚力為260 kPa，有效內(nèi)摩角為51°，泊松比為0.1，飽和滲透系數(shù)為0.02 m/d。由于地基為不透水巖基，因此不考慮地基對古城墻的影響。邊界條件設(shè)置：模型右側(cè)邊界施加水平方向位移約束，邊坡底部為全約束邊界和不透水邊界。接觸設(shè)置：古城墻與土體采用面面接觸，接觸面法向模擬采用硬接觸，并設(shè)置摩擦系數(shù)為0.3。古城墻體和土體采用本構(gòu)模型和屈服準(zhǔn)則，分別為理想彈塑性本構(gòu)模型和Mohr-Coulomb強(qiáng)度準(zhǔn)則；選用CPE6M六結(jié)點(diǎn)修正二次平面應(yīng)變?nèi)切螁卧嬎憔W(wǎng)格布置如圖2所示。

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

圖2 有限元模型的計算網(wǎng)格

2.2 土壓力計算分析

圖3中黑色實(shí)心圓表示由土壓力計測得的離地面4.5 m、離古城墻墻頂3 m處土壓力。由圖3可知，ABAQUS模擬值在深度小于1 m時土壓力為零，在深度大于1 m時土壓力隨著深度的增加而增大，呈三角形分布，符合土力學(xué)理論。土壓力計的實(shí)測值與ABAQUS模擬值接近，而庫倫理論值偏大，這主要因?yàn)槔碚撚嬎銢]有考慮黏聚力對土的影響，且忽略了墻體與填土之間相互作用以及位移協(xié)調(diào)作用。因此，基于ABAQUS軟件對降雨入滲條件下古城墻墻后土體孔隙水壓力和墻體穩(wěn)定性進(jìn)行模擬是可行的。

圖3 水平應(yīng)力隨深度的變化

2.3 有倒掛井的古城墻孔隙水壓力變化過程

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

圖4 有倒掛井的古城墻橫斷面(單位：m)

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

圖5 有倒掛井的古城墻墻背孔隙水壓力等值線

3 具有排水設(shè)施的古城墻降雨滲流分析

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

圖6 有排水設(shè)施的古城墻橫斷面(單位：m)

3.1 孔隙水壓力變化過程

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

圖7 有排水設(shè)施的古城墻墻背孔隙水壓力等值線

3.2 排水設(shè)施對滲流場的影響

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

圖8 持續(xù)降雨24 h后無排水設(shè)施古城墻孔隙水壓力

由圖9可見，降雨滲流可通過排水管排出，從而形成滲流通道，而無排水設(shè)施的古城墻因土體已處于飽和狀態(tài)以致于后續(xù)的降雨滲流溢出。

圖9 降雨24 h后古城墻流速矢量比較

4 古城墻穩(wěn)定性分析

4.1 古城墻抗傾穩(wěn)定性分析

在有排水設(shè)施、無排水設(shè)施和倒掛井3種情形下，降雨24 h后的古城墻接觸應(yīng)力與深度的關(guān)系如圖10所示。從圖10中可以看出，不同情形接觸應(yīng)力隨深度變化的規(guī)律基本一致，隨深度的增加而增大，近似呈三角形分布，但三者在數(shù)值上存在差異，尤其與有排水設(shè)施差異較大。有排水設(shè)施的古城墻后的地下水位比其他2種情形都低，其對古城墻的應(yīng)力也都小于另外2種情形。另外，倒掛井是采用“堵”水的方式，在計算結(jié)果上與無排水設(shè)施情形相近。

由表1可知，由于排水設(shè)施的存在，墻后孔隙水壓力對古城墻的壓力大幅度減小，相比其他2種情形，有排水設(shè)施古城墻受到的總壓力和傾覆力矩最小，也是最穩(wěn)定和最安全的。有倒掛井情形，墻背所受的壓力和傾覆力矩都有所減小，但是效果并不明顯，因此，設(shè)倒掛井并不能有效提高古城墻的穩(wěn)定性。

表1 降雨24 h后古城墻抗傾穩(wěn)定性

4.2 古城墻抗滑穩(wěn)定性分析

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

圖11 不同情形古城墻B點(diǎn)水平位移隨時間的變化

5 結(jié) 論

古城墻無排水設(shè)施時，歷經(jīng)24 h降雨后，墻后孔隙水壓力增大，致使墻背所受壓力增大，降低了古城墻的穩(wěn)定性；有排水設(shè)施時，墻后孔隙水壓力大幅度減小，增強(qiáng)了古城墻的穩(wěn)定性；有防滲設(shè)施時，墻背所受壓力并不能有效降低，安全系數(shù)未明顯增大。