加筋土擋墻地震穩(wěn)定性破裂面隨機(jī)搜索法

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

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

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

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

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

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

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

1.1 水平土條的劃分方法

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

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

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

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

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

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

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

滑動(dòng)楔體水平方向的整體平衡方程為

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

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

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

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

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

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

式中： γ 為土的重度.

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

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

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

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

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

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

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

1.3 多線段破裂面隨機(jī)搜索

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

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

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

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

1.4 多線段破裂面計(jì)算特點(diǎn)

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

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

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

2 多線段破裂面的合理性驗(yàn)證與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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