深厚軟土區樁柱式橋墩臨界荷載簡化計算方法

劉曉明 ，李曦 ，劉齊建 ，劉濤 ，殷悅 ，曾磊

（1.湖南大學 土木工程學院，湖南 長沙，410082；2.中交公路規劃設計院有限公司，北京，100088；3.廣州市高速公路有限公司，廣東 廣州，510032）

橋梁樁基屈曲問題是橋梁設計中不可避免的關鍵問題之一，對于深厚軟土地區尤其如此.20世紀20年代以來，學者們對樁基穩定性問題進行了大量的理論、試驗和數值分析研究，得到了豐富的成果.趙明華等［1-3］基于能量法推導出樁柱式橋墩在不同邊界條件下的計算長度公式，開展了室內模型箱試驗驗證其準確性，并提出了有限元分析樁基屈曲問題的方法.彭錫鼎［4］基于伽遼金法考慮了樁側土抗力，推導出兩端固支樁基的計算長度公式.鄒新軍等［5］對樁基屈曲問題引入新型無單元Galerkin 數值分析法，建立了無單元網格計算程序.雷勇等［6］依據能量法建立了穿越單層溶洞基樁總勢能方程，并引入尖點突變理論，導出了分岔集方程，提出了穿越單層溶洞基樁屈曲臨界荷載計算方法.但是，要將這些科研成果應用于工程，需要開發復雜的程序或軟件，因此設計中往往傾向于簡化方法.對于樁柱式橋墩，基于歐拉公式導出的簡化計算公式［7］簡單明了，是目前設計中最常用的計算方法［因其被寫入《建筑樁基技術規范》（JGJ 94—2008）［8］，下文簡稱該方法為規范法］，但是分析表明，它在理論上存在缺陷，計算結果并不準確.

經過多年發展，基于特征值屈曲分析的有限元方法已經在結構穩定性分析中得到廣泛的應用［9］，其結果也被證明與實際相符.本文擬在深厚軟土區樁柱式橋墩臨界荷載分析中引入特征值屈曲分析的有限元方法對樁柱式橋墩進行分析，在對比說明規范法的適用范圍后，基于機器學習方法整理出簡化計算公式供設計同行參考.

1 有限元法的特征值屈曲分析

1.1 計算原理

特征值屈曲分析的有限元方法是將樁身劃分為若干單元，考慮每個單元所受的節點荷載、樁側土抗力、自重、摩阻力，建立相應于該單元的樁土體系總勢能方程，利用勢能駐值原理，獲得各單元的剛度矩陣［k］e，將單元剛度矩陣組成與樁頂荷載P有關的總剛度方程后得到：

式中：［K］為總剛度矩陣，由單元剛度矩陣［k］e組裝而成；｛δ｝為結構在荷載作用下的位移；｛P｝為結構上的荷載.

式（1）中，總剛度矩陣［K］的逆矩陣等于其伴隨矩陣除以系數的行列式，當總剛度系數行列式等于零時，即：

此時，位移｛δ｝將趨于無窮大，說明結構喪失了穩定性，將發生失穩破壞，故式（2）通常稱為樁身屈曲穩定的有限元特征方程.求解該n階方程可得到n個特征值，對應n階的屈曲模態，在實際工程中，當結構的荷載達到特征值中的最小值時即發生破壞，故其中最小特征值即樁身的臨界荷載Pcr.

1.2 現場試驗簡介

采用作者早年完成的位于洞庭湖深厚軟土區的茅草街大橋超長嵌巖灌注樁屈曲試驗對有限元模型進行標定.現場地質條件如表1 所示［10］，樁周土體水平抗力系數m≈5 500 kN·m-4，樁徑1 m，計算寬度b1=1.8 m，混凝土等級C30，彈性模量E=31 GPa，計算可得樁土變形系數α=0.381 9 m-1.限于篇幅，試驗過程不再贅述，可參閱相應文獻.

表1 樁周土層力學指標Tab.1 Soil mechanical indicators around the pile

試驗加載至17 280 kN 并維持加載約5 min 時，樁身沉降突然增加，地下同時發出一聲很悶的巨響，之后試樁偏位，地面出現裂縫，樁身沉降無法穩定，可判定樁身折斷，其屈曲的臨界荷載可取16 720 kN.

1.3 有限元計算

使用有限元軟件ABAQUS 對該現場試驗進行建模分析.模型采用三維可變形的C3D8R單元建立，樁基本構模型采用彈性模型，彈性模量E設為3.1 GPa；土層的本構模型采用摩爾庫侖模型，具體參數取值如表1 所示.樁土之間設置接觸面，法向行為設置為硬接觸、切向行為設置為罰函數，摩擦系數設為0.5.土體四周及底面設置位移邊界條件以約束土體位移，樁底邊界條件為嵌固，樁端自由.同時為提高計算精度，將樁與樁周土體進行進一步的加密劃分，模型建立如圖1所示.

圖1 樁-土模型示意圖Fig.1 Pile-soil model schematic diagram

模型的材料屬性、網格、接觸屬性、邊界條件建立后模擬現場試驗的各級加載，創建靜力學功能下的geo分析步對地應力進行平衡，之后創建靜力通用分析步，將荷載分級加于樁頂，得到各級荷載下樁身軸力，繪制深度-軸力曲線，并將現場試驗實測數據以散點形式在圖上對應位置標出，如圖2 所示.圖2曲線顯示，現場試驗實測數據點基本落于深度-軸力曲線上，說明了模型參數的正確性以及模型計算的可行性與準確性.

圖2 深度-軸力曲線Fig.2 Depth-axial force curves

然后分析該地層條件下樁的臨界荷載，創建線性攝動功能下的buckle 分析步計算樁的臨界荷載，計算得其臨界荷載Pcr=17 006 kN，樁發生屈曲破壞時，樁、土模態如圖3所示.

圖3 樁、土變形示意圖Fig.3 Pile and soil deformation diagram

計算表明，采用有限元方法計算結果與現場試驗實測值接近，說明有限元方法計算樁基屈曲臨界荷載是可靠的.

1.4 按規范簡化公式計算

常用的簡化計算方法［7］（下文稱：常用簡化公式）過程如下：首先按式（3）確定計算長度lp的取值.

式中：l0為柱高，m；α為樁土變形系數，m-1.

然后將lp代入式（4），即可得到規范法下的臨界荷載值.

式中：E為樁彈性模量，Pa；I為樁橫截面極慣性矩，m4.

對于茅草街大橋試驗基樁，根據式（3）可得lp=7.33 m，再根據式（4）可得臨界荷載值Pcr=223 610 kN，這與實測值（Pcr=16 720 kN）相差很大，表明常用簡化公式計算結果并不準確，這是因為該公式假定樁基在4/α深度處已經嵌固，故將計算長度定為柱高與嵌固深度之和再乘邊界條件系數，樁嵌固深度以下的部分不受影響.事實上，不同種類土的力學性質差異很大，這一假定顯然是不嚴謹的，因此本文將對該方法的適用范圍進行研究，然后提出新的簡化計算方法.

2 常用簡化公式適用范圍分析

2.1 數值實驗方案

一般來說，深厚軟土區的地層分布特點為：上層的軟土層（包括淤泥、淤泥質土、流塑性黏土、軟塑黏性土等），中層的硬土層（包括砂層、粉砂層、硬塑黏土層、卵石層等），下層的巖石層.圖4 所示地層可代表典型深厚軟土區的地層.本文基于該圖所示地層研究典型軟土區橋梁樁基臨界荷載變化規律.

圖4 模型參數示意圖Fig.4 Model parameter schematic diagram

考慮對臨界荷載有影響的土層水平抗力系數m值、樁徑d、樁長l、柱高l0這4 種因素進行數值實驗.其中，樁長l由軟土層中樁長l1與硬土層中樁長l2兩部分組成（l=l1+l2）；m由軟土層中m1與硬土層中m2兩部分組成.改變各參數的取值水平，分別采用常用簡化公式法和基于特征值分析的有限元法計算各參數組合下樁柱式橋墩臨界荷載值，模型參數示意圖如圖4 所示.模型中各參數取值范圍如表2 所示，計算時將各個因素取不同水平組合進行計算.

表2 參數取值Tab.2 Parameter values

由于有限元軟件中沒有直接定義土的水平抗力系數的參數，故需要將土層水平抗力系數m值與彈性模量E值進行換算，其方法為將有限元軟件各彈性模量下樁的水平載荷試驗產生的樁頂水平位移與規范法各m值下計算的樁頂水平位移進行對比，將計算結果非常接近的彈性模量E與m值對應起來，形成若干組數據，運用最小二乘法進行回歸分析，以得到彈性模量E與m值的對應關系［11］，本文模型將樁周土層簡化為單層土，這樣能將土層參數與某一確定m值對應起來，本文中m值與彈性模量E的對應關系如表3所示.

表3 m與E對應取值Tab.3 Corresponding values between m and E

2.2 計算結果分析

將相同參數組合下，兩種不同方法計算所得臨界荷載繪制成散點圖（規范法Pcrs-有限元法Pcra），如圖5所示（長細比=l1/d+l2/d）.

圖5 規范法Pcrs-有限元法Pcra計算結果Fig.5 Results between code method Pcrs and FEM Pcra

計算結果顯示，長細比為7.5～20時，兩種方法的計算結果接近；長細比為20～45 時，兩種方法的計算結果已經較為接近；長細比為60～45 時，兩種方法的計算結果偏差很大.這說明隨著樁長細比的增大，規范法計算結果與有限元計算的結果偏差越來越大，故需要對規范法的計算公式進行修正.

3 臨界荷載簡化計算公式

3.1 隨機森林算法簡介

近年來，人工智能發展迅速，機器學習逐漸應用到越來越多的專業領域，其中隨機森林算法（Random Forest）對于數據的適應能力強，實現方法較為簡便，被廣泛地應用于數據處理當中.隨機森林是一種可用于回歸分析的非參數集成學習算法［12］，起源于對于決策樹的研究［13-14］，通過隨機選擇訓練樣本與特征變量的方法生成多個CART 決策樹［15］，之后通過投票計分的方式將這些決策樹的訓練結果結合起來，以獲得最終的回歸模型，訓練好的模型將可以通過輸入自變量參數，較為準確地預測出因變量參數.

評價隨機森林模型的回歸效果，可采用平均絕對值誤差（MAE）、平均均方根誤差（RMSE）和擬合優度（R2）這3個評價指標［13］，本文采取擬合優度（R2）指標進行評價.評估方法采用“留出法”，將用有限元方法計算出來的數據劃分為兩個數據集合，把其中一個集合作為訓練集S，另一個作為測試集T，在訓練集S上訓練好模型之后，再在測試集T上評估其預測精度，若兩個集合的預測精度均在0.99以上，則說明模型訓練效果較好，模型預測精度較高.

3.2 隨機森林模型預測結果

為得到深厚軟土區樁柱式橋墩計算長度的簡化公式，以第2 節中的數值模型的計算結果為基礎，將結果換算為計算長度，與輸入參數構成數據集.隨機選取數據集66.7%的數據作為訓練集S，33.3%的數據作為測試集T［16］，以長細比、柱高、樁土變形系數α值、土層水平抗力系數m值作為輸入特征變量，采用隨機森林回歸算法對計算長度進行預測.將隨機森林模型得到的預測值與有限元法的計算值繪制散點圖，如圖6所示.

圖6 隨機森林模型預測值與有限元計算值比較Fig.6 Comparison of predicted values from random forest model and FEM calculations

結果顯示，訓練集與測試集的計算結果均處于1∶1 線附近，測試集的擬合優度也達到了0.99 以上，深度學習的效果非常好，證明輸入參數與輸出參數之間具有很強的相關性，可使用公式的形式對其進行擬合.

3.3 基于隨機森林模型的計算長度簡化公式

對于隨機森林模型，根據Gini 指數計算特征變量，可分析各個屬性的重要性［15］，以各因素的重要性指數為基礎可以確定擬合公式中參數的形式從而導出公式.計算本文隨機森林模型中各特征的重要性，并繪制柱狀圖如圖7所示.

圖7 各影響因素重要性Fig.7 Importance of each influencing factor

計算結果顯示，在深厚軟土區的樁柱式橋墩的穩定性問題中，樁的長細比對于計算長度的影響占主導地位，其中軟土層中樁長細比l1/d占44.79%，硬土層中樁長細比l2/d占24.76%，合計69.55%.m值的影響次之，其中軟土層m1占15.15%，硬土層m2占5.19%，合計20.34%.柱高l0對計算長度也有一定的影響，占10.11%；因此以式（3）為基本形式，納入l/d參數可得到更準確的簡化計算公式.本文采用不同形式的公式進行分析，得到樁柱式橋墩計算長度簡化公式如（5）所示，將擬合結果繪制在圖8 中，可見計算較好.

圖8 擬合曲面Fig.8 Fitted surface

根據式（5）計算各參數水平下計算長度，然后采用式（4）計算樁柱式橋墩的臨界荷載，再與有限元法計算值比較，繪制散點圖（2 520 個數據點）如圖9所示.

圖9 修正公式預測值與有限元計算值比較Fig.9 Comparison of predicted values from modified formula and FEM calculations

計算結果表明，本文提出的簡化計算方法效果良好，具有實用價值.

4 結論

本文分別采用基于特征值屈曲分析的有限元數值軟件法、規范法，對樁柱式橋墩的屈曲臨界荷載進行計算和比較，并使用隨機森林法對有限元計算結果進行分析，得到結論如下：

1）采用有限元數值軟件方法計算樁柱式橋墩的臨界荷載，其計算結果與現場試驗結果較為吻合，說明基于特征值屈曲分析的有限元法對樁柱式橋墩的臨界荷載計算是適用的.

2）現有常用簡化計算方法在樁的長細比較大、時計算結果誤差較大，適用性差.

3）基于有限元計算以及隨機森林回歸分析，擬合得到了深厚軟土區樁柱式橋墩的臨界荷載計算簡化方法，可供設計人員使用與參考.