李 倩 ，凌天清 ，韓林峰 ，張瑞剛

（1.重慶交通大學(xué)土木工程學(xué)院，重慶 400074；2.重慶交通大學(xué)重慶市山區(qū)道路建設(shè)與維護技術(shù)重點實驗室，重慶 400074；3.水利部長江水利委員會河湖保護與建設(shè)運行安全中心，湖北 武漢 430010）

加筋土擋墻因節(jié)省成本、易于施工、設(shè)計靈活、承受大變形且不易破壞的優(yōu)點而廣泛適用于土木工程中.加筋土擋墻必須具有足夠的穩(wěn)定性，以抵抗各種可能的外部荷載作用.目前，地震荷載作用下加筋土擋墻的穩(wěn)定性分析多采用擬靜力法，該方法假設(shè)墻后填土的破裂面為平面，將擋土墻所受到地震加速度等效為水平和豎直方向的單向加速度[1].在應(yīng)用擬靜力法進行穩(wěn)定性計算時，由于加筋土擋墻的水平成層特點，一般結(jié)合水平條分法（HSM）進行分析，該方法由Lo 等[2]基于極限平衡理論首次提出.Shahgholi 等[3]提出了一種簡化的HSM，將整個加筋區(qū)分成若干個水平土條，只考慮每個土條的豎向平衡和整體水平向平衡（不考慮力矩平衡），將N層水平土條的平衡方程個數(shù)減少為 2N+1 個.Nouri 等[4]提出了一種不同的HSM，稱為 5N-1 公式，采用該方法計算了加筋土擋墻維持內(nèi)部穩(wěn)定所需的筋材臨界長度.Shekarian 等[5]提出了一種用于加固回填土剛性擋土墻抗震分析的HSM 方法.

雖然擬靜力法應(yīng)用廣泛，但該方法沒有考慮地震波特性，為此，提出了擬動力計算方法.擬動力計算方法多采用正弦波模擬地震加速度，但只考慮垂直向上傳播的正弦波，沒有考慮向上傳播的波與自由地面反射波之間的相互作用，由于該方法不考慮地面對加速度的放大影響，需在分析中假設(shè)一個合適的放大系數(shù).Choudhury 等[6]采用擬動力法確定了地震作用下加筋土擋墻抵抗直接滑移破壞和傾覆破壞所需的加筋長度.Bellezza[7]將均質(zhì)無黏性土體假設(shè)為Kelvin-Voigt 固體，通過求解一維波動方程得到填土體內(nèi)的加速度，并計算了墻體所受土壓力.徐鵬等[8]分析了黏土地基對加筋土擋墻穩(wěn)定性的影響，提出了黏土地基加筋土擋墻的設(shè)計建議.

采用HSM 進行穩(wěn)定性分析的關(guān)鍵是確定加筋土擋墻的破裂面位置和形態(tài)，Sil 等[9]采用圓弧狀破裂面分析了均勻黏性土質(zhì)邊坡在地震作用下的動力性能.Prater[10]假定滑動面為圓形和對數(shù)螺旋形，研究了均質(zhì)土坡的屈服加速度.在實際情況下，加筋土擋墻破裂面的位置往往是未知的，破裂面也不是光滑的，學(xué)者們研究了加筋土結(jié)構(gòu)破裂面的確定方法.蔣薇等[11]基于極限平衡理論采用優(yōu)化算法求得加筋土邊坡臨界滑動面.林永亮等[12]以極限分析理論為基礎(chǔ)，假定破裂面為對數(shù)螺旋形，研究了加筋土邊坡的臨界高度.鄧東平等[13]以簡化Janbu 法為基礎(chǔ)采用滑動面搜索法計算了邊坡的整體穩(wěn)定性.Hu 等[14]利用變尺度混沌優(yōu)化算法確定了非圓弧滑移面.Malkawi 等[15]采用蒙特卡羅優(yōu)化方法計算了邊坡臨界滑動面.Sarma 等[16]采用一種附加應(yīng)力容許準(zhǔn)則的極限平衡方法確定了邊坡的臨界滑動面.萬文等[17]將加速混合遺傳算法應(yīng)用于邊坡最危險滑動面求解之中.Ausilio 等[18]采用極限分析原理結(jié)合對數(shù)螺旋破裂面研究了加筋土擋墻的地震穩(wěn)定性.Nimbalkar等[19]以水平條分法為基礎(chǔ)結(jié)合回歸分析對加筋土擋墻進行了擬動力穩(wěn)定性分析，Basha 等[20]采用對數(shù)螺旋破裂面對加筋土擋墻進行了擬動力計算.

本文以一個隨機變量為基礎(chǔ)，采用隨機搜索法結(jié)合水平條分法，對加筋土擋墻破裂面進行搜索計算，生成了一個多線段形式的破裂面，得到了加筋土擋墻臨界破裂面的位置.考慮墻后填土為黏彈性介質(zhì)，采用黏彈性地震加速度計算地震慣性力，對加筋土擋墻進行了地震穩(wěn)定性擬動力分析.

1 加筋土擋墻破裂面的確定方法

1.1 水平土條的劃分方法

水平條分法計算模型如圖1 所示，圖中：點A、C分別為坡面頂點和坡腳點；點B為墻后填土表面的破裂點；α為擋墻傾角；H為擋墻高度；D為水平土條厚度；n為筋材總層數(shù)；Vi、Vi+1分別為作用于第i個水平土條上、下兩側(cè)豎直條間力；Hi、Hi+1分別為作用于第i個水平土條上、下兩側(cè)水平條間力；qhi、qvi分別為第i個水平土條水平、豎直向地震慣性力；zi和Wi分別為第i個水平土條的位置深度和質(zhì)量；Si、Ni分別為第i個水平土條破裂面上的切向與法向力；βi為第i個水平土條破裂面與水平方向的夾角；Ti為第i個水平土條中筋材的拉力；li=D/sin βi為第i個水平土條的破裂面長度.

圖1 水平條分法計算模型Fig.1 Calculation model of horizontal slice method

根據(jù)極限平衡原理，第i個水平土條的水平向與豎直向平衡方程為

考慮平衡條件，水平土條在破裂面處存在式（3）表示關(guān)系.

式中： τf為失效剪應(yīng)力；τr為維持平衡所需剪應(yīng)力；Fs為安全系數(shù)，假定對所有水平土條安全系數(shù)取相同值.

根據(jù)Mohr-Coulomb 破壞準(zhǔn)則，第i個水平土條破裂面上Si和Ni滿足式（4）所示關(guān)系.

式中： φ 為土的內(nèi)摩擦角；c為土的黏聚力.

滑動楔體水平方向的整體平衡方程為

將式（1）～（4）代入式（5），令Fs=1，筋材總拉力為

將墻后填土視為黏彈性介質(zhì)，采用Kelvin-Voigt模型，根據(jù)橫波與縱波在平面內(nèi)傳播的運動方程可得Kelvin-Voigt 黏彈性介質(zhì)運動方程為[21-23]

式中：uh為土體水平位移；uv為土體豎向位移；t為時間；z為土體深度；ρ 為土體密度；G為土體剪切模量；λ 為拉梅常數(shù)；ηs為橫波傳播黏性系數(shù)；ηp為縱波傳播黏性系數(shù).

Bellezza[21]利用阻尼比公式和歸一化的橫波頻率公式，并假設(shè)擋墻模型下方的剛性層受到諧波激勵，地表的剪切應(yīng)力和法向應(yīng)力為0，擋墻底部與剛性基底的位移一致，得出地震波在Kelvin-Voigt 黏彈性土體中的水平加速度為

式中：ωs為橫波頻率；kh為水平地震加速度；

同理，地震波在Kelvin-Voigt 黏彈性土體中的豎向加速度為

式中： γ 為土的重度.

作用在第i個水平土條上的水平向地震慣性力和豎直向地震慣性力分別如式（11）和式（12）所示.

將式（11）、（12）代入式（6），可求出筋材總拉力，筋材總拉力可用一類似于土壓力的無量綱參數(shù)K表示，即

1.2 多線段破裂面形狀及生成

假定加筋土擋墻滑動破裂面由多個傾斜線段組成，這些線段在一個平面內(nèi)以不同的長度和角度相互連接，如圖2 所示.多線段破裂面的傾斜線段數(shù)與筋材層數(shù)相同，若加筋層數(shù)有n層，即有n個水平土條，則滑動破裂面中的傾斜線段有n條，隨著層數(shù)的增加，多線段滑動破裂面逐漸逼近曲線形態(tài).

圖2 多線段破裂面形式Fig.2 Form of polyline fracture surface

多線段破裂面采用產(chǎn)生隨機角的方式生成，如圖3 所示.Lc為填土表面破裂點與坡面頂部之間的距離，∠CBE為生成隨機角所需的最大角度，其值可由已知幾何條件求得.圖3 中擋墻加筋層數(shù)有5 層，則程序采用產(chǎn)生的隨機數(shù)生成5 條虛線，每條虛線與相應(yīng)的水平土條界面的交點為1、2、3、4 以及坡腳點C，將所有的交點相連可得到多線段破裂面形狀.利用產(chǎn)生的隨機角計算得到 θ1～ θ4的值，得到每條傾斜線段的傾斜角度，可計算每層水平土條中的筋材拉力.

圖3 多線段破裂面生成方法Fig.3 Generation method of polylinefracture surface

1.3 多線段破裂面隨機搜索

加筋土擋墻多線段破裂面搜索程序采用Fortran語言編制.程序首先采用random_seed 函數(shù)產(chǎn)生隨機數(shù)所需的種子，然后在隨機角生成循環(huán)步的初始時刻，根據(jù)種子產(chǎn)生隨機數(shù)，采用random_number函數(shù)將各隨機數(shù)與最大角度相乘，便可得到各隨機角的角度值，進而可得到每條傾斜線段的傾斜角度.根據(jù)傾斜線段利用水平條分法便可計算各水平土條中的筋材拉力.

如圖4 所示，多線段破裂面搜索程序有兩層循環(huán)計算，外層循環(huán)為墻后水平填土表面破裂點位置循環(huán)，內(nèi)層循環(huán)為隨機角度循環(huán).為計算準(zhǔn)確，在初始計算時需合理選擇初始破裂點，在選定初始破裂點之后程序便開始產(chǎn)生隨機角并計算每層筋材的拉力，進一步計算所有筋材的總拉力，為保證計算的精確性，需適當(dāng)設(shè)置隨機角循環(huán)次數(shù).

圖4 破裂面搜索計算流程Fig.4 Flow chart of fracture surface search calculation

在一次內(nèi)循環(huán)計算結(jié)束時記錄筋材總拉力最大值，最大值所對應(yīng)的破裂面形狀即為此次循環(huán)所得破裂面形狀.一次內(nèi)循環(huán)計算結(jié)束后，程序便選擇新的破裂點進行下一次外循環(huán)計算，并記錄此次外循環(huán)中筋材總拉力的最大值，如此往復(fù)，每選擇一次破裂點便會計算得到一次筋材總拉力的最大值.程序最后對計算得到的筋材總拉力進行比較，最大值所對應(yīng)的破裂面形狀即為所求的加筋土擋墻臨界破裂面.

1.4 多線段破裂面計算特點

在以往的研究中，加筋土擋墻地震穩(wěn)定性計算多是預(yù)先假定滑動面形狀或采用數(shù)學(xué)優(yōu)化算法對破裂面位置進行求解.相比于預(yù)先確定破裂面形狀的計算方法，本文提出的滑動破裂面隨機搜索計算方法更具靈活性，且無需進行優(yōu)化計算.

采用隨機角度結(jié)合水平條分法的計算方法便于計算機編程，且計算方法中只有一個隨機變量，計算過程簡單，各計算因素易于掌握，相比于采用數(shù)學(xué)優(yōu)化方法求解，具有計算簡便的特點.

水平條分法不是本文所提破裂面計算方法應(yīng)用的必需條件，該破裂面計算方法也可用于豎向條分分析方法，結(jié)合Janbu 法、Sarma 法等經(jīng)典方法可進行邊坡穩(wěn)定性分析.對數(shù)螺旋破裂面的滑動軌跡與土的內(nèi)摩擦角有關(guān)，該破裂面的應(yīng)用僅限于均質(zhì)土層問題，而本文所提出的計算方法可用于非均質(zhì)土層的穩(wěn)定性分析.

2 多線段破裂面的合理性驗證與分析

2.1 與筋材總拉力和Lc/H的比較

為便于將本文方法與現(xiàn)有方法的計算結(jié)果進行對比分析，計算參數(shù)取值相同，設(shè)置 H=5 m，α=0°，φ=30°，Vs=100m/s，γ=18kN/m3，ξ=10%，kv=0，計算結(jié)果如表1 所示.

表1 筋材總拉力對比Tab.1 Comparison of total reinforcement forces kN/m

由表1 可知：地震作用下本文算法得到的筋材總拉力最大，產(chǎn)生這種差異的原因在于本文考慮了地震加速度的非線性分布.文獻[3]采用擬靜力計算方法進行計算，沒有考慮加速度放大的因素；文獻[19]的擬動力方法雖考慮了加速度相位的變化，但假設(shè)加速度變化幅值為常數(shù)；文獻[20]認為加速度沿地面線性增加，但方法不滿足邊界處的應(yīng)力條件，即加速度的放大值不受約束.

為分析計算方法的可靠性，在 kh=0.2時計算Lc/H隨 φ的變化規(guī)律，并與文獻[4,18]進行了對比，如圖5 所示.由圖可知：本文計算方法所得參數(shù)變化規(guī)律與文獻[4,18]相同，即在相同條件下，Lc/H隨 φ的增大而減小，擋墻傾角越小，Lc/H越小.本文算法計算結(jié)果偏小，所得墻后填土表面破裂點到坡面頂部的距離最小.文獻[4]采用對數(shù)螺旋破裂面結(jié)合水平條分法，文獻[18]采用對數(shù)螺旋破裂面結(jié)合能量原理對 Lc/H進行了分析.文獻[4,18]雖然采用了相同的破裂面形式，但由于采用了不同的計算方法造成計算結(jié)果略有差距.本文計算結(jié)果小于文獻[4,18]，主要是破裂面形態(tài)不同所致.

圖5 Lc/H 與φ關(guān)系的比較Fig.5 Comparison of relationship between Lc/H and φ

2.2 與對數(shù)螺旋破裂面位置的比較

對數(shù)螺旋破裂面是加筋土擋土結(jié)構(gòu)在穩(wěn)定性計算時常用到的一種破裂面形狀，如圖6 所示.對數(shù)螺旋破裂面方程可用式（14）表示.

圖6 對數(shù)螺旋破裂面形式Fig.6 Form of logarithmic spiral fracture surface

式中： θ 為極坐標(biāo)系中的旋轉(zhuǎn)角度；θ0為初始角度；θh為總角度；r為對數(shù)螺旋線的半徑；r0為初始半徑.

對數(shù)螺旋破裂面采用 θ0與 θh兩個最重要的參數(shù)進行定義.在計算時，通過確定維持平衡所需筋材拉力最大值，來搜尋對數(shù)螺旋破裂面的位置.利用MATLAB 工具箱中fminsearch 函數(shù)，采用優(yōu)化方法求解筋材拉力計算方程，使得筋材拉力達到最大值，求解過程中的約束條件為 0°＜θ0＜90°，0°＜θh＜90°.當(dāng) α=15°時，對數(shù)螺旋破裂面與多線段破裂面位置計算結(jié)果如圖7 所示.由圖可知：多線段破裂面計算方法所得加筋土擋墻的破裂面位置與對數(shù)螺旋破裂面較接近，但更靠近臨空面，隨著填土內(nèi)摩擦角的減小兩者計算結(jié)果差距略有增大，在此算例中最大相差9.6%.

圖7 多線段破裂面與對數(shù)螺旋破裂面位置的比較Fig.7 Position comparison of polyline fracture surfaces and logarithmic spiral fracture surfaces

本文所提多線段破裂面計算方法與對數(shù)螺旋破裂面的不同之處在于，多線段破裂面采用隨機角的方式生成，隨機角的產(chǎn)生不受破裂面方程的制約.對數(shù)螺旋破裂面在優(yōu)化過程中，主要優(yōu)化θ0與θh兩個角度參數(shù)，優(yōu)化過程中破裂面受到對數(shù)螺旋方程的限制.而采用隨機角生成的多線段破裂面，每一段破裂形狀不受破裂面方程的限制，求解過程中形式變化靈活，對筋材總拉力的求解更為精確.由此導(dǎo)致多線段破裂面位置的求解結(jié)果相對于對數(shù)螺旋破裂面更接近擋墻臨空面.

3 穩(wěn)定性影響因素分析

假定H=5m，ξ=10%，γ=18 kN/m3，α 為0°～30°，φ為25°～45°，kh為0.1～0.2，計算不同地震加速度和擋墻傾角下，φ 與K、Lc/H的關(guān)系分別如圖8、9 所示.

圖8 α不同時K 隨φ的變化規(guī)律Fig.8 Trend of K changing with φ under different α

3.1 α、kh 與K 的關(guān)系

由圖8 可知：地震作用下，K隨著 φ 的增大而減小，隨著加筋土擋墻坡度的減小而減小，隨著地震加速度的增大而增大.當(dāng) α=0°時，φ由2 5°增大到45°，kh由0.1增加到0.2 時，K增加了約19.7%和23.2%；在φ為25°且kh為0.2 情況下，擋墻傾角從0°增加到15°和從15°增加到 30°時，K分別減小了約21.3% 和37.7%；kv對K的影響不大，尤 其當(dāng)φ較大時；當(dāng)φ=45°，α=0°時，kh= 0.2，相比于kh=0.1時，K僅增加了8.3%.

3.2 α、kh與Lc/H 的關(guān)系

由圖9 可知：地震作用下，Lc/H隨著加筋土擋墻坡度的減小而減小，隨著 φ 的增大而減小，該結(jié)論也可從文中對比分析算例圖5 中得到驗證.當(dāng)不考慮α與φ的變化時，當(dāng)kh=0.1，kv對Lc/H的影響可忽略；當(dāng)kh=0.2 時，kv的變化對Lc/H的影響相比于K更明顯.

圖9 α不同時 Lc/H 隨φ的變化規(guī)律Fig.9 Trend of L c/H changing with φ under different α

4 結(jié) 論

1） 將加筋土擋墻破裂面視為多線段形式，基于水平條分法提出了破裂面生成新方法，以筋材總拉力作為目標(biāo)函數(shù)，采用角度隨機的方式對目標(biāo)函數(shù)進行搜索計算，確定加筋土擋墻破裂面的位置.

2） 所提出的方法由于考慮了地震加速度的非線性分布，計算得地震作用下筋材總拉力較大.與對數(shù)螺旋破裂面相比，多線段破裂面計算得加筋土擋墻破裂面位置更接近臨空面，當(dāng)填土內(nèi)摩擦角為25°時，兩者位置最大相差 9.6%.

3） 地震作用下，筋材總拉力隨著加筋土擋墻坡度和地震加速度的增大而增大，筋材長度隨著填土內(nèi)摩擦角的增大而減小.當(dāng)填土內(nèi)摩擦角較大時，豎向加速度的變化對筋材總拉力和筋材長度的影響較小.