基于混合粒子群算法和能量法主動土壓力計算

陳昌富，唐仁華，梁冠亭

（湖南大學 巖土工程研究所，長沙 410082）

1 引 言

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

2 基于能量法和對數螺線滑裂面的主動土壓力計算

能量法計算土壓力需用到以下兩個結論[8]：

（1）流動法則：土體處于塑性流動或剪切滑動狀態時，滑動面上任一點的應變速度矢量v 與該點處的滑動線成φ 角。

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

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

圖1 對數螺線滑動面計算圖 Fig.1 Calculation diagram of log-spiral sliding surface

假定滑動面為通過墻底的對數螺線滑動面，其動徑方程為

式中：rθ為0r 繞中心旋轉θ 角后的動徑長度(m)；0θ為0r 與水平面的夾角(°)；0r 為對數螺線通過邊坡的旋入點到螺線中心的長度(m)。

當θ =hθ 時，

式中：hr 為對數螺線通過邊坡的旋出點至螺線中心的長度(m)；hθ 為hr 與水平面的夾角(°)。

過O 點作AB 的平行線與CA 的延長線交與點E；作AC 的平行線與BA 的延長線交于點F，則

由式（3）可以推出

此外，對數螺線參數還滿足式（5）中的幾何關系

式中：xo為對數螺線中心橫坐標(m)；yo為對數螺線中心縱坐標(m)；L 為滑裂面起點到坡頂的距離；xAB為L 在x 軸負方向的投影。

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

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

圖2 滑體發生微小轉動計算圖 Fig.2 Calculation diagram of small rotation on slide soil

圖3 滑體發生微小轉動計算圖 Fig.3 Calculation diagram of small rotation on slide soil

超載q 所作的功為

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

沿對數螺線微分長度 rθdθ / cosφ 上的消耗功為( rθdθ / cosφ ) cv cosφ = crθvd θ，內部總消耗功為

圖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 面上單元土體速度矢量合成圖，可見其平均速度即為中點Z 處的速度。土壓力合力P 與墻面法線夾角為δ (外摩擦角)，位于法線下方，與Z 點速度方向夾角為π-(hθ -zθ )-β -δ ，則擋土墻對AC 面土體的反力所做功為

由于擋土墻對土體的反力與土體的應變速率方向相反，故其所做之功為負值。由外力作的功等于土體內部消耗的能量，消去角速度可得：

式（17）即為主動土壓力計算公式。顯然，這是一個二元函數的極值問題，對于主動土壓力，應取其最大值。

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

式中：

不同深度處的主動土壓力強度pH為

主動土壓力的作用點距墻底的高度為

式中：Pa為主動土壓力，kaγ，kac，kaq分別為單獨考慮重度，黏聚力，超載時的土壓力系數。

3 基于混合粒子群算法主動土壓力優化計算

基于能量法主動土壓力計算可歸結為如下優化模型：

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

式（25）是一個含有復雜的目標函數和約束條件的非線性規劃問題，本文采用適用性強且具有良好的全局搜索能力的粒子群算法[9](PSO)求解。

在PSO 算法中，每個粒子都具有3 個參數：位置矢量X，代表目標函數中自變量的值；速度矢量V，代表粒子飛翔的力向和距離，用來引導粒子不斷朝著更好的解空間改變當前位置；適應度fitness，由被優化的目標函數決定，目標函數值越小，則適應度越好。在每一次迭代中，粒子通過跟蹤2 個最優解來更新自己，一是粒子本身所找到的最優解，另一個是整個種群目前所找到的最優解。

混合粒子群算法是將自然選擇機制與粒子群算法結合，在每次迭代過程中將整個粒子群按適應值排序，用群體中最好的一半的粒子的速度和位置替換最差的一半的位置和速度，同時保留原來每個個體所記憶的歷史最優值。

基于以上原理和模型，主動土壓力混合粒子群優化算法計算步驟如下：

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

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

（3）將各微粒的位置和適應值存儲在各微粒的個體最好位置pbest 中，將所有pbest 中適應值最優個體的位置和適應值存儲于群體最好位置gbest 中；

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

式中：t 為迭代的步數；pi為個體最優適應值對應的[xi,1， xi,2]值，即第i 個粒子最優適應值對應的的[0θ ，hθ ]值；pg為群體最優適應值對應的[0θ ，hθ ]值；w 為慣性權重因子，文中取w=0.5；c1、c2為正的學習因子，文中取c1= c2=2；r1、r2為[0，1]范圍內的均勻隨機數。

（5）對每個微粒，將其適應值與其經歷過的最好適應值比較，如果較好，則將其對應的位置作為當前的最好位置，然后比較當前所有pbest 和gbest的值，更新gbest；

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

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

4 結果對比分析

當墻背豎直時，對于無超載的砂性土，將c =0，q=0，β =90°代入式(20)，運用上述優化方法計算土壓力系數kaγ隨填土面傾角α 的變化結果，并與文獻[3]中結果進行對比。

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

當填土面水平時(α =0°)，土壓力系數kaγ隨墻背傾角β 的變化與文獻[3]結果對比如表2 所示。

表1 墻背豎直時與文獻[3]結果對比 Table 1 Comparison between results obtained by reference [3] and with a vertical wall

表2 填土面水平時與文獻[3]結果對比 Table 2 Comparison between results obtained reference [3] and with horizontal backfills

從表2 中選擇部分數據繪成更直觀的圖5。

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

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

表3 墻背向內傾斜30°時土壓力系數對比 Table 3 Comparison of active earth pressure coefficients while the wall angle is 30°

從表3 中可以看出，當墻背傾角較大時，本文比文獻[3]更接近文獻[10]的計算結果，筆者認為這是由不同的破壞機構導致的。由于文獻[3]沒有對對數螺線中心進行搜索，所以在一些特定情況下不能搜索到土壓力系數的最優值（最大值），從而小于實際值。顯然，此時的庫倫理論精度更小，因為墻背傾角較大時，真實滑動面與平面滑動面的假定存在較大的差異。

5 實例驗證

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

計算結果見表4。同時給出了建筑地基基礎規范[12]推薦方法的計算結果與實測結果。

表4 本文計算結果 Table 4 Result of the article

表5 本文計算結果與實測值對比 Table 5 Comparison of results obtained by the article and measured date

顯然，本文對黏性土的計算結果比規范推薦方法計算結果更接近實測值。

6 結 論

（1）本文基于能量法和混合離子群優化算法，假定滑動面為對數螺線滑裂面，得到了一種新的主動土壓力及其系數計算方法。

（2）對于黏性土，通過與一重力式擋墻試驗結果對比分析發現，本文結果與實測值誤差僅為5.4%，表明本文方法計算結果可靠。

（3）對于無超載砂性土的土壓力系數，本文方法計算結果與經典的極限分析上限解在墻面傾角較小時基本一致，但當墻面傾角較大時，經典解明顯偏小，與本文解最大差別達54.71%，而本文解與基于最優性原理的極限平衡解較接近。

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

[2] 趙樹德， 趙樹惠. 用對數螺線滑動面計算擋土墻主動土壓力[J]. 西安建筑科技大學學報， 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] 范文， 沈珠江， 愈茂宏. 基于統一強度理論的土壓力極限上限分析[J]. 巖土工程學報， 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] 楊建民. 基于非線性破壞準則的主動土壓力上限計 算[J]. 巖土力學， 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] 顧慰慈. 擋土墻土壓力計算手冊[M]. 北京： 中國建材出版社， 2005.

[9] 龔純， 王正林. 精通Matlab 最優化計算[M]. 北京： 電子工業出版社， 2009.

[10] 朱大勇， 錢七虎. 極限平衡法計算土壓力系數的新途徑[J]. 土木工程學報， 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] 四川省建筑科學研究院. 重力式擋土墻的試驗研究[R]. 成都： 四川省建筑科學研究院， 1982.

[12] 中華人民共和國國家標準. 建筑地基基礎設計規范[S].北京： 中國建筑工業出版社， 2002.