基于改進(jìn)SW滑楔模型剛性抗滑樁極限滑坡推力智能優(yōu)化計(jì)算方法

陳昌富，戴宇佳，梁冠亭，朱劍鋒

(1.湖南大學(xué)巖土工程研究所，湖南長沙 410082;2.寧波大學(xué)建筑工程與環(huán)境學(xué)院，浙江寧波 315211)

工程中遇到的滑坡大多都是沿著一個或多個軟弱層面而產(chǎn)生滑動(圖1)。對于此類滑坡，如果滑動面相對較淺，通常是采用剛性抗滑樁加固。一般剛性樁的位移量較小，如果滑體在樁位處的剩余下滑力大于樁后土體所能提供的最大被動土壓力，則隨著滑體沿軟弱層面不斷滑移，樁后土體會逐漸形成一個三維被動滑動土楔(圖1)，此時，作用于樁上的極限滑坡推力就變?yōu)橛杀粍尤S滑動土楔所產(chǎn)生的推力(即被動土壓力)。因此，對于沿一軟弱層面的滑坡，在滑體處于極限平衡狀態(tài)下，作用于剛性抗滑樁上的極限滑坡推力，應(yīng)是由滑體所產(chǎn)生的剩余下滑力和由樁后三維滑動土楔所產(chǎn)生的被動土壓力兩者中較小者。剩余下滑力的計(jì)算目前已比較成熟，并已形成規(guī)范［1～2］，而對于樁后三維滑動土楔所產(chǎn)生的被動土壓力的計(jì)算問題，目前尚未完善，因此本文嘗試借用最近國內(nèi)外關(guān)注較多的應(yīng)變楔(SW)模型［3］來解決。

圖1 邊坡分析簡圖Fig．1 Diagram for the slope analysis

應(yīng)變楔(SW)模型最早由Norris等［3］提出，并用其分析了樁頂受水平集中荷載作用樁，得到了在三維土楔體作用下樁身受力分布;李彬［4］考慮了應(yīng)力集中現(xiàn)象，假定作用在楔形體面上的應(yīng)力變化Δσh在水平切面上呈三角形分布，并運(yùn)用于群樁受力分析;李忠誠等［5］將應(yīng)變楔模型應(yīng)用于被動樁，建立了被動拱-主動楔模型，實(shí)現(xiàn)了三維樁土相互作用計(jì)算;Youngho Kim等［6］采用三維非線性有限元分析，驗(yàn)證了楔體模型的合理性，并進(jìn)行了多組水平受荷樁試驗(yàn)，并在應(yīng)變楔模型基礎(chǔ)上考慮楔體兩側(cè)的土體受力，得到了樁上受力沿樁長呈雙曲線分布;Ashour等［7］將應(yīng)變楔模型與彈性地基梁法結(jié)合，運(yùn)用到邊坡穩(wěn)定分析中，通過迭代求出彈性樁上的推力分布。

現(xiàn)有SW模型中，均假定滑動楔體底部滑面為平面，這與實(shí)際情況往往不符。并且，由于邊坡坡角越大，最大側(cè)向變形越大［8］，SW模型未考慮斜坡的影響，計(jì)算結(jié)果不夠準(zhǔn)確。為此，本文結(jié)合傳遞系數(shù)法，以SW模型為基礎(chǔ)，假定滑楔底面為對數(shù)螺旋線面，得到改進(jìn)SW滑楔模型。同時，基于傾斜薄層單元法，引入智能優(yōu)化算法——粒子群優(yōu)化(PSO)算法，計(jì)算極限滑坡推力Pult。結(jié)合算例分別探討樁位、樁徑以及土體參數(shù)的變異系數(shù)對極限滑坡推力Pult的影響。

1 受力分析模型

基于SW模型，為便于分析，改進(jìn)SW滑楔模型作了如下假定:

(1)剛性樁加固的滑坡，樁身不發(fā)生位移，邊坡滑體沿軟弱層面產(chǎn)生足夠的滑移，使樁后土體達(dá)到被動極限平衡狀態(tài)，形成被動滑楔體(圖1)。

(2)被動滑楔體以yOz面對稱(圖2)，其兩側(cè)滑面為從樁背沿y軸向外擴(kuò)展的平面，底部破裂面為繞x軸轉(zhuǎn)動的對數(shù)螺線曲面(圖3)。

圖2 計(jì)算模型示意圖Fig．2 Schematic diagram of the calculating model

圖3 對數(shù)螺旋線滑動面Fig．3 Slip surface composed of the log-spiral surface

(3)將被動滑楔體斜分為薄層單元，單元間沿x軸方向切向?qū)娱g力因?qū)ΨQ相互抵消，只考慮豎向作用的層間力。

(4)作用于樁上的推力Pi沿深度呈線性分布。

1.1 破裂面參數(shù)的確定

如圖3所示，以樁頂處O點(diǎn)為坐標(biāo)原點(diǎn)，建立直角坐標(biāo)系，將土體分成n個薄層單元。

對數(shù)螺旋線滑動面方程為:

式中:ri——第i個薄層滑面對數(shù)螺線的半徑;

r0——滑面對數(shù)螺線的初始半徑;

a——系數(shù)，可取 a=tan2(45°-2φ/3)［9］;

φ——土的內(nèi)摩擦角;

θi——第i個薄層單元的旋轉(zhuǎn)角度。

設(shè)對數(shù)螺線在第i個薄層的斜率為ki，滑動面近似為直線［10］，對應(yīng)的滑動面傾角為 αi，初始半徑 r0與水平面的夾角為ζ，則:

1.2 樁上滑坡推力計(jì)算

樁上滑坡推力計(jì)算模型如圖2。取其中第i個薄層單元進(jìn)行受力分析(圖4)。斜薄層單元在其斜平面上呈梯形狀，其上底寬等于樁的寬度，下底寬設(shè)為li。若設(shè)斜薄層單元水平投影長度為bi，則根據(jù)幾何關(guān)系可知:

圖4 薄層單元受力分析Fig．4 Thin-layer element analysis

式中:D——樁徑(m);

φm——土的傘角。

由于底部第n個薄層單元在縱剖面上近似為三角形，于是有:

薄層單元沿y軸、z軸方向上力的平衡方程為:

式中:Pi——樁上滑坡推力，與水平面夾角為δ;

Ei，Ei+1——豎向作用的層間力;

Nsi，Tsi——滑楔側(cè)面法向力與切向力;

Nbi，Tbi——滑楔底面法向力與側(cè)向力。Nbi表現(xiàn)為拉力，取為0。

滑動面上土體處于極限平衡狀態(tài)，切向力與法向力服從Mohr-Coulomb定律:

式中:ci，φi——滑動面上土體抗剪強(qiáng)度指標(biāo);

Asi，Abi——滑楔側(cè)面和底面的表面積。

聯(lián)立式(8)～(11)，可得:

假定樁側(cè)剪應(yīng)力強(qiáng)度完全發(fā)揮，τi=(τult)i，綜合考慮粘性土中應(yīng)變軟化、應(yīng)力歷史及不排水強(qiáng)度比對樁側(cè)抗剪強(qiáng)度的影響，根據(jù)文獻(xiàn)［11］知:

由于破裂面為任意假定對數(shù)螺旋線面，取不同的破裂面可求得不同的滑坡推力合力P值，可通過優(yōu)化方法求出最大滑坡推力Pmax，即為極限滑坡推力Pult，對應(yīng)的破裂面即為最危險破裂面。

2 PSO搜索破裂面

基于本文提出的改進(jìn)SW滑楔模型，樁上極限滑坡推力Pult計(jì)算可歸結(jié)為如下優(yōu)化模型:

本文采用具有良好全局搜索能力的粒子群算法［12］(PSO)搜索破裂面。PSO模擬鳥群的捕食行為，每個優(yōu)化問題的解都是搜索空間中的一個“粒子”。PSO初始化為一群粒子(隨機(jī)解)，通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個最優(yōu)解來更新自己，一是粒子本身所找到的最優(yōu)解，另一個是整個種群目前所找到的最優(yōu)解。

根據(jù)以上模型和原理，滑坡推力混合粒子群算法計(jì)算步驟如下:

(1)設(shè)定粒子群數(shù)目，初始化每個粒子的位置矢量 Xi=［xi，1，xi，2］和速度矢量 Vi=［vi，1，vi，2］;

(2)根據(jù)每個粒子當(dāng)前的Xi計(jì)算其適應(yīng)度。若gi(X)＞0成立，則fitness(i)=f(Xi);否則，將一個遠(yuǎn)大于極限滑坡推力的值賦給適應(yīng)度函數(shù)作為該粒子的適應(yīng)值以使其受到懲罰，例如fitness(i)=10000;

(3)對每個微粒，將它的適應(yīng)值與其經(jīng)歷過的最好位置pbest的適應(yīng)值作比較，如果較好，則更新pbest，將它的適應(yīng)值和全局所經(jīng)歷最好位置gbest的適應(yīng)值作比較，如果較好，則更新gbest;

(4)根據(jù)式(18)和(19)更新各粒子的速度和位置，取每個微粒的最好適應(yīng)值對應(yīng)的位置作為當(dāng)前的最好位置，更新gbest;

(5)設(shè)定最大迭代次數(shù)，達(dá)到最大迭代次數(shù)后則輸出結(jié)果，即滑坡推力最大時的ζ、θn值以及極限滑坡推力Pult大小。

3 算例與參數(shù)分析

采用文獻(xiàn)［7］中給出的兩個分析算例(圖5)，邊坡的水平投影距離為L，抗滑樁間距S=3m，樁位至坡腳水平投影距離為Lx。邊坡巖土體參數(shù)指標(biāo)如圖5，其中滑面抗剪強(qiáng)度參數(shù)與土體的相同。

圖5 邊坡示意圖(引自文獻(xiàn)［5］)Fig．5 Schematic diagram of slope(cited in［5］)

3.1 樁位的影響

取樁徑D=1.2m，分別采用傳遞系數(shù)法、SW法和本文方法計(jì)算不同樁位處的極限滑坡推力Pult，計(jì)算結(jié)果如圖6所示。

由圖6可知，當(dāng)滑體處于極限平衡狀態(tài)時，Pult隨樁位(Lx/L)變化呈現(xiàn)中間大兩頭小的趨勢。在樁位靠近坡頂處，由于滑體產(chǎn)生的剩余下滑力相對較小，樁后不能發(fā)展出如圖1所示穿出滑體的被動三維滑動土楔，此時的Pult就是滑體所產(chǎn)生的剩余滑坡推力(本文采用傳遞系數(shù)法計(jì)算)，而如果采用SW法計(jì)算就會高估作用于樁上的極限滑坡推力(最大誤差達(dá)76%)。當(dāng)樁位處在坡的中下部，滑體產(chǎn)生的剩余下滑力相對較大，樁后土體能夠發(fā)展出如圖1所示被動土楔，此時Pult就等于被動土楔所產(chǎn)生的被動土壓力，若仍采用傳遞系數(shù)法計(jì)算作用于樁上極限滑坡推力，則計(jì)算結(jié)果偏大，且最大誤差處在坡中，相對誤差為16%。

圖6 極限滑坡推力Pult隨樁位Lx/L的變化Fig．6 Variation in the ultimate landslide thrust Pultwith Lx/L

圖6進(jìn)一步說明，當(dāng)樁間距較大，樁位設(shè)置在坡的中部時，作用于樁上的實(shí)際極限滑坡推力小于按規(guī)范［13］推薦的傳遞系數(shù)法所得計(jì)算結(jié)果，因此，根據(jù)規(guī)范法計(jì)算的滑坡推力來設(shè)計(jì)抗滑樁，偏于保守。而且，已有研究結(jié)果表明［14］:與有限元數(shù)值法和Morgenstern-Price法相比，采用規(guī)范的傳遞系數(shù)法計(jì)算滑坡推力，所設(shè)計(jì)的抗滑樁也偏于保守。

3.2 樁徑D的影響

將樁設(shè)置于Lx/L=0.5處，樁徑D從0.9 m逐漸增大至1.6 m，計(jì)算水平極限推力Pult，結(jié)果如圖7所示。

圖7 極限滑坡推力Pult隨樁徑的變化Fig．7 Variation in the ultimate landslide thrust Pult with pile diameter

由圖7可知，Pult與樁徑D呈線性關(guān)系，隨D的增大而增大。與SW法相比，Pult隨樁徑變化的分布一致，樁徑增大對Pult影響較小，Pult增長幅度小于SW法計(jì)算結(jié)果的增長幅度。

4 土體參數(shù)變異性對Pult的影響

以圖5邊坡(a)為例，將樁設(shè)置于Lx/L=0.5處，土體的重度γ1，粘聚力c1，以及內(nèi)摩擦角φ1的變異系數(shù)分別為 δγ1，δc1，δφ1，分析上述 3 個參數(shù)的變異性對極限滑坡推力Pult的變異系數(shù)δPult的影響。

假定各參數(shù)服從正態(tài)分布，選取土體參數(shù)不同的變異系數(shù)分別計(jì)算Pult的變異系數(shù)δPult(圖8)。

圖8 不同變異系數(shù)對Pult的影響Fig．8 Influences of variable coefficients to Pult

由圖8 可知，Pult的變異系數(shù) δPult隨土體參數(shù) γ1、c1和 φ1的變異系數(shù) δγ1、δc1和 δφ1的增大而增大。當(dāng) δγ1=0.025 ～0.1、δc1=0.1 ～0.4、δφ1=0.1 ～0.25 時，極限滑坡推力Pult的變異系數(shù)δPult在0.06～0.2之間變化。由此法得到的極限滑坡推力Pult變異系數(shù)δPult，可用于抗滑樁嵌固深度的可靠性分析，從而可簡化抗滑樁可靠性設(shè)計(jì)計(jì)算。

5 結(jié)論

(1)本文借鑒應(yīng)變楔(SW)模型，假定破裂面為對數(shù)螺旋線面，并引入傾斜薄層單元法，推導(dǎo)了極限滑坡推力計(jì)算公式，并利用粒子群優(yōu)化(PSO)算法搜索計(jì)算極限滑坡推力Pult。

(2)算例計(jì)算結(jié)果表明:當(dāng)樁間距較大，樁位設(shè)置在坡的中部時，作用于樁上的極限滑坡推力Pult小于按規(guī)范推薦的傳遞系數(shù)法所得計(jì)算結(jié)果，而且Pult與樁徑基本上呈線性關(guān)系。

(3)假定各參數(shù)服從正態(tài)分布，根據(jù)對土體參數(shù)不同變異系數(shù)進(jìn)行抽樣計(jì)算，可得到極限滑坡推力Pult的變異系數(shù) δPult，由此法得到的變異系數(shù) δPult結(jié)果可直接用于剛性抗滑樁可靠性分析計(jì)算。

(4)考慮土體應(yīng)力-應(yīng)變關(guān)系，可建立非極限滑坡推力與Pult、土體位移之間的關(guān)系式，于是可求得非極限滑坡推力，由于本文篇幅有限，作者將另文闡述。

［1］ 中華人民共和國住房和城鄉(xiāng)建設(shè)部.GB50330-2013建筑邊坡工程技術(shù)規(guī)范［S］.北京:中國建筑工業(yè)出版社，2014.［Ministry of Housing and Urban-Rural Development of the People＇s Republic of China.GB50330-2013 Technical code for building slope engineering［S］. Beijing:China Architecture ＆Building Press，2014.(in Chinese)］

［2］ 中華人民共和國住房和城鄉(xiāng)建設(shè)部.GB50007-2011建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范［S］.北京:中國建筑工業(yè)出版社，2012.［Ministry of Housing and Urban-Rural Development of the People＇s Republic of China.GB50007- 2011 Code fordesign ofbuilding foundation［S］. Beijing:China Architecture ＆Building Press，2012.(in Chinese)］

［3］ ASHOUR M，NORRIS G，Pilling P.Lateral loading of a pile in layered soil using the strain wedge model［J］.Journal of Geotechnical And Geoenvironmental Engineering，ASCE，1998，124(4):303-305.

［4］ 李彬.改進(jìn)SW法及其在樁基水平受荷分析中的應(yīng)用［D］.杭州:浙江大學(xué)，2002.［LI B.Analysis Of Laterally Loaded Piles With Modified SW Model［D］.Hangzhou:Zhejiang University，2002.(in Chinese)］

［5］ 李忠誠，楊敏.被動樁土壓力計(jì)算的被動拱-主動楔模型［J］.巖石力學(xué)與工程學(xué)報，2006，25(增刊2):4241-4247.［LI Z C，YANG M.Passive ArchingActive Wedge ModelOfSoilPressure Calculation In Passive Piles［J］.Rock Mechanics and Engineering，2006，25(Sup2):4241-4247.(in Chinese)］

［6］ YOUNGHO Kim，SANGSEOM Jeong，M ASCE，et al.Wedge Failure Analysis of Soil Resistance on Laterally Loaded Piles in Clay［J］. Journalof Geotechnical and Geoenvironmental Engineering，2011，137(7):678-694.

［7］ MOHAMED Ashour，HAMED Ardalan.Analysis of pile stabilized slopes based on soil-pile interaction［J］.Computers and Geotechnics，2012，39:85-97.

［8］ 魏煥衛(wèi)，孫建平，陳啟輝，等.基坑邊坡變形的理論計(jì)算方法［J］.水文地質(zhì)工程地質(zhì)，2006，33(2):75-79.［WEI H W，SUN J P，CHEN Q H，et al.Theory Calculation Method of Slope Displacement［J］.Hydrogeology＆ Engineering Geology，2006，33(2):75-79.(in Chinese)］

［9］ 趙樹德，趙樹惠.用對數(shù)螺旋線滑動面計(jì)算擋土墻主動土壓力［J］.西安建筑科技大學(xué)學(xué)報(自然科學(xué)版)，2002，34(4):343-345.［ZHAO S D，ZHAO S H.Calculation of active earth pressure of retaining wall with logistic spiral sliding surface［J］.Journal ofXi’an University ofArchitecture ＆Technology(Natural Science Edition)，2002，34(4):343-345.(in Chinese)］

［10］ 陳昌富，曾玉瑩，肖淑君，等.基于薄層單元法主動土壓力計(jì)算的復(fù)合遺傳算法［J］.巖土力學(xué)，2006，27(3):398-403.［CHEN C F，ZENG Y Y，XIAO S J，et al.Complex genetic method of active earth pressure based on thin-layer element method［J］.Rock and Soil Mechanics，2006，27(3):398-403.(in Chinese)］

［11］ API(1993)Recommendedpracticeforplanning，designing and constructing fixed offshore platforms—working stress design［S］.API RP2A，20th ed.Washington DC:American Petro-leum Institute.

［12］ 龔純，王正林.精通Matlab最優(yōu)化計(jì)算［M］.北京:電子工業(yè)出版社，2009.［GONG C，WANG Z L.Proficient in the optimization based on matlab［M］.Beijing:Publishing House of Electronics Industry，2009.(in Chinese)］

［13］ 交通部第二勘察設(shè)計(jì)院.公路設(shè)計(jì)手冊(路基)［M］.北京:中國鐵道出版社，2001.［Second Survey and Design Institute of the Ministry of Communications.Highway Design Manual(Subgrades)［M］.Beijing:China Railway Press，2001.(in Chinese)］

［14］ 程建軍，廖小平，王浩，等.滑坡推力計(jì)算方法的對比研究與應(yīng)用［J］.水文地質(zhì)工程地質(zhì)，2008，35(1):44-48.［CHENG J J，LIAO X P，WANG H，etal. Comparative research and application of landslide thrust calculation methods ［J］.Hydrogeology ＆ Engineering Geology，2008，35(1):44-48.(in Chinese)］