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

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

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

經(jīng)過多年發(fā)展，基于特征值屈曲分析的有限元方法已經(jīng)在結(jié)構(gòu)穩(wěn)定性分析中得到廣泛的應(yīng)用［9］，其結(jié)果也被證明與實(shí)際相符.本文擬在深厚軟土區(qū)樁柱式橋墩臨界荷載分析中引入特征值屈曲分析的有限元方法對樁柱式橋墩進(jìn)行分析，在對比說明規(guī)范法的適用范圍后，基于機(jī)器學(xué)習(xí)方法整理出簡化計算公式供設(shè)計同行參考.

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

1.1 計算原理

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

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

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

此時，位移｛δ｝將趨于無窮大，說明結(jié)構(gòu)喪失了穩(wěn)定性，將發(fā)生失穩(wěn)破壞，故式（2）通常稱為樁身屈曲穩(wěn)定的有限元特征方程.求解該n階方程可得到n個特征值，對應(yīng)n階的屈曲模態(tài)，在實(shí)際工程中，當(dāng)結(jié)構(gòu)的荷載達(dá)到特征值中的最小值時即發(fā)生破壞，故其中最小特征值即樁身的臨界荷載Pcr.

1.2 現(xiàn)場試驗(yàn)簡介

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

表1 樁周土層力學(xué)指標(biāo)Tab.1 Soil mechanical indicators around the pile

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

1.3 有限元計算

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

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

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

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

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

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

計算表明，采用有限元方法計算結(jié)果與現(xiàn)場試驗(yàn)實(shí)測值接近，說明有限元方法計算樁基屈曲臨界荷載是可靠的.

1.4 按規(guī)范簡化公式計算

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

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

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

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

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

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

2.1 數(shù)值實(shí)驗(yàn)方案

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

圖4 模型參數(shù)示意圖Fig.4 Model parameter schematic diagram

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

表2 參數(shù)取值Tab.2 Parameter values

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

表3 m與E對應(yīng)取值Tab.3 Corresponding values between m and E

2.2 計算結(jié)果分析

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

圖5 規(guī)范法Pcrs-有限元法Pcra計算結(jié)果Fig.5 Results between code method Pcrs and FEM Pcra

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

3 臨界荷載簡化計算公式

3.1 隨機(jī)森林算法簡介

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

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

3.2 隨機(jī)森林模型預(yù)測結(jié)果

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

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

結(jié)果顯示，訓(xùn)練集與測試集的計算結(jié)果均處于1∶1 線附近，測試集的擬合優(yōu)度也達(dá)到了0.99 以上，深度學(xué)習(xí)的效果非常好，證明輸入?yún)?shù)與輸出參數(shù)之間具有很強(qiáng)的相關(guān)性，可使用公式的形式對其進(jìn)行擬合.

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

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

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

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

圖8 擬合曲面Fig.8 Fitted surface

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

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

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

4 結(jié)論

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

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

2）現(xiàn)有常用簡化計算方法在樁的長細(xì)比較大、時計算結(jié)果誤差較大，適用性差.

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