基于混合粒子群算法和能量法主動(dòng)土壓力計(jì)算

陳昌富，唐仁華，梁冠亭

（湖南大學(xué) 巖土工程研究所，長(zhǎng)沙 410082）

1 引 言

自Terzaghi[1]使用對(duì)數(shù)螺線滑動(dòng)面模型求解擋土墻土壓力后，此模型經(jīng)歷了幾十年試驗(yàn)和實(shí)踐的考證，被普遍認(rèn)為是最接近實(shí)際的破壞機(jī)構(gòu)[2]。Chen[3]較早地系統(tǒng)研究了極限分析法在土壓力中的應(yīng)用，隨后許多學(xué)者運(yùn)用塑性極限上限方法來(lái)計(jì)算土壓力，取得了較多成果。范文等[4]和楊建民[5]利用多三角形破壞機(jī)構(gòu)，分別基于統(tǒng)一強(qiáng)度理論和非線性強(qiáng)度準(zhǔn)則求得了土壓力的上限解。楊小禮[6]運(yùn)用非線性強(qiáng)度準(zhǔn)則分析了墻后填土面水平時(shí)動(dòng)、靜主動(dòng)土壓力上限解。Chen[3]分別采用6 種破壞機(jī)構(gòu)結(jié)合最速下降法求得了土壓力系數(shù)，其中2 條直線中間夾對(duì)數(shù)螺線破壞機(jī)構(gòu)的計(jì)算結(jié)果與廣泛接受的Sokoloskii[7]解誤差最小，Chen 認(rèn)為，是對(duì)數(shù)螺線夾層起到了關(guān)鍵作用。但上述機(jī)構(gòu)假定對(duì)數(shù)螺線的中心固定于坡頂不變，而實(shí)際情況并非如此。本文采用接近實(shí)際情況的單純對(duì)數(shù)螺線滑動(dòng)面破壞機(jī)構(gòu)，不約束對(duì)數(shù)螺線中心的位置，推導(dǎo)出土壓力及其系數(shù)計(jì)算公式，然后運(yùn)用基于自然選擇的混合粒子群算法對(duì)最危險(xiǎn)滑動(dòng)面進(jìn)行全局搜索并計(jì)算主動(dòng)土壓力，最后將本文方法計(jì)算結(jié)果與已有解答和實(shí)測(cè)值進(jìn)行對(duì)比分析。

2 基于能量法和對(duì)數(shù)螺線滑裂面的主動(dòng)土壓力計(jì)算

能量法計(jì)算土壓力需用到以下兩個(gè)結(jié)論[8]：

（1）流動(dòng)法則：土體處于塑性流動(dòng)或剪切滑動(dòng)狀態(tài)時(shí)，滑動(dòng)面上任一點(diǎn)的應(yīng)變速度矢量v 與該點(diǎn)處的滑動(dòng)線成φ 角。

（2）土體在塑性流動(dòng)中的能量消散為cvcosφ。c、φ 分別為土體的黏聚力和內(nèi)摩擦角。

如圖1 所示的擋土墻，墻高為H，墻背與水平面成β 角，墻后填土與水平面成α 角。填土上作用均布荷載q。以坡頂處A 點(diǎn)為坐標(biāo)原點(diǎn)，向左和向上分別定為x，y 坐標(biāo)軸的正方向。

圖1 對(duì)數(shù)螺線滑動(dòng)面計(jì)算圖 Fig.1 Calculation diagram of log-spiral sliding surface

假定滑動(dòng)面為通過(guò)墻底的對(duì)數(shù)螺線滑動(dòng)面，其動(dòng)徑方程為

式中：rθ為0r 繞中心旋轉(zhuǎn)θ 角后的動(dòng)徑長(zhǎng)度(m)；0θ為0r 與水平面的夾角(°)；0r 為對(duì)數(shù)螺線通過(guò)邊坡的旋入點(diǎn)到螺線中心的長(zhǎng)度(m)。

當(dāng)θ =hθ 時(shí)，

式中：hr 為對(duì)數(shù)螺線通過(guò)邊坡的旋出點(diǎn)至螺線中心的長(zhǎng)度(m)；hθ 為hr 與水平面的夾角(°)。

過(guò)O 點(diǎn)作AB 的平行線與CA 的延長(zhǎng)線交與點(diǎn)E；作AC 的平行線與BA 的延長(zhǎng)線交于點(diǎn)F，則

由式（3）可以推出

此外，對(duì)數(shù)螺線參數(shù)還滿足式（5）中的幾何關(guān)系

式中：xo為對(duì)數(shù)螺線中心橫坐標(biāo)(m)；yo為對(duì)數(shù)螺線中心縱坐標(biāo)(m)；L 為滑裂面起點(diǎn)到坡頂?shù)木嚯x；xAB為L(zhǎng) 在x 軸負(fù)方向的投影。

作用在滑動(dòng)土體上的外力有均布荷載載q，滑動(dòng)土體的重力W 已及擋土墻對(duì)滑動(dòng)土體的反作用力p。假設(shè)滑動(dòng)土體繞對(duì)數(shù)螺線中心產(chǎn)生微小的角速度ω，由假定的對(duì)數(shù)螺線方程可知，螺線上每點(diǎn)的切線方向與動(dòng)徑成(π/2 φ+ )角，結(jié)合流動(dòng)法則，則滑動(dòng)面上每點(diǎn)的應(yīng)變速度方向均與動(dòng)徑垂直，即

如圖2、3 所示，由于超載和重力的方向都是豎直向下，故其與速度方向夾角等于動(dòng)徑與水平線的夾角。即超載或重力做的功等于力與角速度ω 的乘積，再乘以超載或者單元土體形心到對(duì)數(shù)螺線中心O 點(diǎn)的水平距離。

圖2 滑體發(fā)生微小轉(zhuǎn)動(dòng)計(jì)算圖 Fig.2 Calculation diagram of small rotation on slide soil

圖3 滑體發(fā)生微小轉(zhuǎn)動(dòng)計(jì)算圖 Fig.3 Calculation diagram of small rotation on slide soil

超載q 所作的功為

如圖4 所示，滑動(dòng)土體重力作的功可由疊加法求出。以W1表示OBC 對(duì)假定轉(zhuǎn)動(dòng)中心O 作的功，W2表示OAB 對(duì)假定轉(zhuǎn)動(dòng)中心O 作的功，W3表示OAC 對(duì)假定轉(zhuǎn)動(dòng)中心O 作的功。OAB 形心到O 點(diǎn)的水平距離可由OAM 與ABM 組合求出，OAC 形心到O 點(diǎn)的水平距離可由OCN 與ACN 組合求出。則W1，W2，W3分別為

沿對(duì)數(shù)螺線微分長(zhǎng)度 rθdθ / cosφ 上的消耗功為( rθdθ / cosφ ) cv cosφ = crθvd θ，內(nèi)部總消耗功為

圖4 AC 面上單元土體速度矢量圖 Fig.4 Soil element’s velocity vector on surface AC

如圖4（a），由余弦定理知

由正弦定理知

式 中：A=sin-1[ H sin(β + θh)/(2 ozsin β)]；G = sin(β +θh)/(2sin β)；J = ( r0/ H )2e2(θh-θ0)tanφ；U=cos(β +θ )/ H sin β。 h

圖4（b）右為AC 面上單元土體速度矢量合成圖，可見(jiàn)其平均速度即為中點(diǎn)Z 處的速度。土壓力合力P 與墻面法線夾角為δ (外摩擦角)，位于法線下方，與Z 點(diǎn)速度方向夾角為π-(hθ -zθ )-β -δ ，則擋土墻對(duì)AC 面土體的反力所做功為

由于擋土墻對(duì)土體的反力與土體的應(yīng)變速率方向相反，故其所做之功為負(fù)值。由外力作的功等于土體內(nèi)部消耗的能量，消去角速度可得：

式（17）即為主動(dòng)土壓力計(jì)算公式。顯然，這是一個(gè)二元函數(shù)的極值問(wèn)題，對(duì)于主動(dòng)土壓力，應(yīng)取其最大值。

特別地，當(dāng)α =0°、β =90°，即墻背豎直、墻后填土面水平時(shí)，整理公式后可得：

式中：

不同深度處的主動(dòng)土壓力強(qiáng)度pH為

主動(dòng)土壓力的作用點(diǎn)距墻底的高度為

式中：Pa為主動(dòng)土壓力，kaγ，kac，kaq分別為單獨(dú)考慮重度，黏聚力，超載時(shí)的土壓力系數(shù)。

3 基于混合粒子群算法主動(dòng)土壓力優(yōu)化計(jì)算

基于能量法主動(dòng)土壓力計(jì)算可歸結(jié)為如下優(yōu)化模型：

式中：X 為優(yōu)化設(shè)計(jì)變量；f (X)為目標(biāo)函數(shù)，計(jì)算土壓力時(shí)f (X)=-Pa，計(jì)算土壓力系數(shù)時(shí)f (X)=-kaγ，其中Pa與kaγ分別按式（17）與式（20）計(jì)算；gj(X)( j=1， 2， 3，4)為約束條件，其中xAB按式（5）計(jì)算。

式（25）是一個(gè)含有復(fù)雜的目標(biāo)函數(shù)和約束條件的非線性規(guī)劃問(wèn)題，本文采用適用性強(qiáng)且具有良好的全局搜索能力的粒子群算法[9](PSO)求解。

在PSO 算法中，每個(gè)粒子都具有3 個(gè)參數(shù)：位置矢量X，代表目標(biāo)函數(shù)中自變量的值；速度矢量V，代表粒子飛翔的力向和距離，用來(lái)引導(dǎo)粒子不斷朝著更好的解空間改變當(dāng)前位置；適應(yīng)度f(wàn)itness，由被優(yōu)化的目標(biāo)函數(shù)決定，目標(biāo)函數(shù)值越小，則適應(yīng)度越好。在每一次迭代中，粒子通過(guò)跟蹤2 個(gè)最優(yōu)解來(lái)更新自己，一是粒子本身所找到的最優(yōu)解，另一個(gè)是整個(gè)種群目前所找到的最優(yōu)解。

混合粒子群算法是將自然選擇機(jī)制與粒子群算法結(jié)合，在每次迭代過(guò)程中將整個(gè)粒子群按適應(yīng)值排序，用群體中最好的一半的粒子的速度和位置替換最差的一半的位置和速度，同時(shí)保留原來(lái)每個(gè)個(gè)體所記憶的歷史最優(yōu)值。

基于以上原理和模型，主動(dòng)土壓力混合粒子群優(yōu)化算法計(jì)算步驟如下：

（1）設(shè)定粒子群數(shù)目為40，即i∈[1，40]。對(duì)每個(gè)粒子的當(dāng)前位置矢量 Xi=[xi,1，xi,2]和速度矢量Vi= [vi,1，vi,2]在變量范圍內(nèi)[0， π-β ]隨機(jī)賦予初始值；

（2）根據(jù)每個(gè)粒子當(dāng)前的iX 計(jì)算其適應(yīng)度。若gi(X)＞0 成立，則fitness(i)=f (iX )；否則，將一個(gè)遠(yuǎn)大于正常土壓值的大數(shù)賦給適應(yīng)度函數(shù)作為此粒子的適應(yīng)值以使其受到懲罰，例如fitness(i)=10 000；

（3）將各微粒的位置和適應(yīng)值存儲(chǔ)在各微粒的個(gè)體最好位置pbest 中，將所有pbest 中適應(yīng)值最優(yōu)個(gè)體的位置和適應(yīng)值存儲(chǔ)于群體最好位置gbest 中；

（4）按式（26）和（27）更新每個(gè)微粒的速度和位置：

式中：t 為迭代的步數(shù)；pi為個(gè)體最優(yōu)適應(yīng)值對(duì)應(yīng)的[xi,1， xi,2]值，即第i 個(gè)粒子最優(yōu)適應(yīng)值對(duì)應(yīng)的的[0θ ，hθ ]值；pg為群體最優(yōu)適應(yīng)值對(duì)應(yīng)的[0θ ，hθ ]值；w 為慣性權(quán)重因子，文中取w=0.5；c1、c2為正的學(xué)習(xí)因子，文中取c1= c2=2；r1、r2為[0，1]范圍內(nèi)的均勻隨機(jī)數(shù)。

（5）對(duì)每個(gè)微粒，將其適應(yīng)值與其經(jīng)歷過(guò)的最好適應(yīng)值比較，如果較好，則將其對(duì)應(yīng)的位置作為當(dāng)前的最好位置，然后比較當(dāng)前所有pbest 和gbest的值，更新gbest；

（6）將整個(gè)粒子群按適應(yīng)值排序，用群體中最好的一半的粒子的速度和位置替換掉最差的一半的位置和速度，保持pbest 和gbest 不變。

（7）設(shè)定最大迭代次數(shù)為100 步，達(dá)到迭代次數(shù)后則輸出結(jié)果，即土壓力最大時(shí)的0θ 、hθ 值以及土壓力或者土壓力系數(shù)的值。

4 結(jié)果對(duì)比分析

當(dāng)墻背豎直時(shí)，對(duì)于無(wú)超載的砂性土，將c =0，q=0，β =90°代入式(20)，運(yùn)用上述優(yōu)化方法計(jì)算土壓力系數(shù)kaγ隨填土面傾角α 的變化結(jié)果，并與文獻(xiàn)[3]中結(jié)果進(jìn)行對(duì)比。

由表1 知，在不同的填土面傾角下，本文結(jié)果與文獻(xiàn)[3]結(jié)果幾乎一致(差值在2%以內(nèi))。說(shuō)明兩者的差別對(duì)填土面傾斜角度不敏感。

當(dāng)填土面水平時(shí)(α =0°)，土壓力系數(shù)kaγ隨墻背傾角β 的變化與文獻(xiàn)[3]結(jié)果對(duì)比如表2 所示。

表1 墻背豎直時(shí)與文獻(xiàn)[3]結(jié)果對(duì)比 Table 1 Comparison between results obtained by reference [3] and with a vertical wall

表2 填土面水平時(shí)與文獻(xiàn)[3]結(jié)果對(duì)比 Table 2 Comparison between results obtained reference [3] and with horizontal backfills

從表2 中選擇部分?jǐn)?shù)據(jù)繪成更直觀的圖5。

圖5 填土面水平時(shí)差值隨墻面傾角的變化(δ =φ /2) Fig.5 Difference varying with wall angle for horizontal backfills(δ =φ /2)

從圖可以看出，當(dāng)墻面與水平面的夾角為90°、80°、70°時(shí)，兩者差值接近于0。而當(dāng)墻背向內(nèi)傾斜角度大于30°時(shí)，兩者存在明顯的差別，尤其是填土摩擦角較大時(shí)，本文比文獻(xiàn)[3]的計(jì)算結(jié)果大50%左右。將兩者結(jié)果與朱大勇等[10]的極限平衡解對(duì)比如表3 所示。

表3 墻背向內(nèi)傾斜30°時(shí)土壓力系數(shù)對(duì)比 Table 3 Comparison of active earth pressure coefficients while the wall angle is 30°

從表3 中可以看出，當(dāng)墻背傾角較大時(shí)，本文比文獻(xiàn)[3]更接近文獻(xiàn)[10]的計(jì)算結(jié)果，筆者認(rèn)為這是由不同的破壞機(jī)構(gòu)導(dǎo)致的。由于文獻(xiàn)[3]沒(méi)有對(duì)對(duì)數(shù)螺線中心進(jìn)行搜索，所以在一些特定情況下不能搜索到土壓力系數(shù)的最優(yōu)值（最大值），從而小于實(shí)際值。顯然，此時(shí)的庫(kù)倫理論精度更小，因?yàn)閴Ρ硟A角較大時(shí)，真實(shí)滑動(dòng)面與平面滑動(dòng)面的假定存在較大的差異。

5 實(shí)例驗(yàn)證

四川建筑科學(xué)研究所在四川簡(jiǎn)陽(yáng)做的重力式擋土墻試驗(yàn)[11]：墻背豎直，墻高為4 m，土天然重度γ =19.34 kN/m3，不排水三軸快剪指標(biāo)φ = 16.63°、δ =8.32°、c=4.7 kPa。

計(jì)算結(jié)果見(jiàn)表4。同時(shí)給出了建筑地基基礎(chǔ)規(guī)范[12]推薦方法的計(jì)算結(jié)果與實(shí)測(cè)結(jié)果。

表4 本文計(jì)算結(jié)果 Table 4 Result of the article

表5 本文計(jì)算結(jié)果與實(shí)測(cè)值對(duì)比 Table 5 Comparison of results obtained by the article and measured date

顯然，本文對(duì)黏性土的計(jì)算結(jié)果比規(guī)范推薦方法計(jì)算結(jié)果更接近實(shí)測(cè)值。

6 結(jié) 論

（1）本文基于能量法和混合離子群優(yōu)化算法，假定滑動(dòng)面為對(duì)數(shù)螺線滑裂面，得到了一種新的主動(dòng)土壓力及其系數(shù)計(jì)算方法。

（2）對(duì)于黏性土，通過(guò)與一重力式擋墻試驗(yàn)結(jié)果對(duì)比分析發(fā)現(xiàn)，本文結(jié)果與實(shí)測(cè)值誤差僅為5.4%，表明本文方法計(jì)算結(jié)果可靠。

（3）對(duì)于無(wú)超載砂性土的土壓力系數(shù)，本文方法計(jì)算結(jié)果與經(jīng)典的極限分析上限解在墻面傾角較小時(shí)基本一致，但當(dāng)墻面傾角較大時(shí)，經(jīng)典解明顯偏小，與本文解最大差別達(dá)54.71%，而本文解與基于最優(yōu)性原理的極限平衡解較接近。

[1] TERZAGHI K. Theoretical soil mechanics[M]. New York： [s. n.]， 1943.

[2] 趙樹(shù)德， 趙樹(shù)惠. 用對(duì)數(shù)螺線滑動(dòng)面計(jì)算擋土墻主動(dòng)土壓力[J]. 西安建筑科技大學(xué)學(xué)報(bào)， 2002， 34(4)： 342－345. ZHAO Shu-de， ZHAO Shu-hui. Calculation of earth pressure of active retaining wall with logarithmic spiral sliding surface[J]. Xi'an Univ. of Arch.& Tech.， 2002， 34(4)： 342－345.

[3] CHEN W F. Limit analysis and soil plasticity[M]. New York： Elsevier Scientific Publishing Company， 1975.

[4] 范文， 沈珠江， 愈茂宏. 基于統(tǒng)一強(qiáng)度理論的土壓力極限上限分析[J]. 巖土工程學(xué)報(bào)， 2005， 27(10)： 1147－1153. FAN Wen， SHENG Zhu-jiang， YU Mao-hong. Upper-bound limit analysis of earth pressure based on unified strength theory[J]. Chinese Journal of Geotechnical Engineering， 2005， 27(10)： 1147－1153.

[5] 楊建民. 基于非線性破壞準(zhǔn)則的主動(dòng)土壓力上限計(jì) 算[J]. 巖土力學(xué)， 2009， 30(2)： 503－508. YANG Jian-min. Upper bound solution for active earth pressure based on a nonlinear failure criterion[J]. Rock and Soil Mechanics， 2009， 30(2)： 503－508.

[6] YANG Xiao-li. Upper bound limit analysis of active earth pressure with different fracture surfaces and nonlinear yield criterion[J]. Theoretical and Applied Fracture Mechanics， 2007， 47(1)： 46－56.

[7] SOKOLOSKII V V. Statics of granular media[M]. New York： Pergamon Press， 1965.

[8] 顧慰慈. 擋土墻土壓力計(jì)算手冊(cè)[M]. 北京： 中國(guó)建材出版社， 2005.

[9] 龔純， 王正林. 精通Matlab 最優(yōu)化計(jì)算[M]. 北京： 電子工業(yè)出版社， 2009.

[10] 朱大勇， 錢七虎. 極限平衡法計(jì)算土壓力系數(shù)的新途徑[J]. 土木工程學(xué)報(bào)， 2002， 33(1)： 63－68. ZHU Da-yong， QIAN Qi-hu. New approach for computation of earth pressure coefficients based on limit equilibrium method[J]. China Civil Engineering Journal， 2002， 33(1)： 63－68.

[11] 四川省建筑科學(xué)研究院. 重力式擋土墻的試驗(yàn)研究[R]. 成都： 四川省建筑科學(xué)研究院， 1982.

[12] 中華人民共和國(guó)國(guó)家標(biāo)準(zhǔn). 建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范[S].北京： 中國(guó)建筑工業(yè)出版社， 2002.