考慮群樁效應的液化場地樁-土動力p-y曲線研究

惠舒清， 張效禹， 唐亮， 崔杰， 滿孝峰， 凌賢長,3

(1.哈爾濱工業大學 土木工程學院，黑龍江 哈爾濱 150090； 2.廣州大學 土木工程學院，廣東 廣州 510006； 3.青島理工大學 土木工程學院，山東 青島 266033)

歷次震害調查表明[1]，強震下場地液化是導致樁基橋梁結構震害的重要原因。因此，國際上對強震下液化場地樁基橋梁抗震問題十分關注，并陸續開展了卓有成效的研究[2-3]。工程實踐與研究表明[4-5]，可液化場地條件下，樁基橋梁結構地震反應分析的關鍵在于精確模擬樁-土動力相互作用?，F有樁-土動力相互作用分析方法中，p-y曲線法具有無需選取土體本構模型、計算參數容易獲取、計算時長很短等特點，在工程設計中得到廣泛使用。

現階段，已有多位學者提出了用于樁-土相互作用分析的p-y曲線，如Reese等[6]提出了砂土的p-y曲線、Matlock等[7]提出了水下軟粘土和硬粘土的p-y曲線、Lombardi等[8]提出了土體液化后的p-y曲線。然而，現有關于砂土p-y曲線的研究多是針對單樁展開的，其對群樁基礎是否適用尚需進一步探討。為此，本文針對液化場地群樁-土-結構耦合體系地震反應離心機振動臺試驗，建立了可靠的三維有限元數值模型；隨后，在反復的數值模擬與分析基礎上，構建了考慮樁間距效應的液化場地群樁-土動力相互作用p-y曲線，并通過離心機試驗結果對其正確性進行了檢驗。

1 液化場地橋梁群樁-土地震相互作用離心機試驗數值模擬

1.1 離心機試驗簡介

本文針對Wilson[4]開展的液化場地橋梁群樁-土地震相互作用離心機試驗進行數值模擬工作，離心加速度為30g，文中給出的參數均為按照離心加速度確定的原型尺寸，試驗布置見圖1。試驗中，地基土采用內華達砂(Nevada sand)，總厚度20.5 m，分為上、下兩層。上層為相對密度約35%～40%的飽和松砂，厚度為9.1 m。下層為相對密度約80%的飽和密砂，厚度為11.4 m。水位線位于地表處。試驗采用真空砂雨法制備地基模型土。試驗中，基礎為2×2群樁基礎，其中，樁采用鋁制管樁，樁長14.5 m，樁徑0.67 m，壁厚19 mm，樁間距2.67 m(約為4倍樁徑)，彈性模量70 GPa，屈服應力290 MPa，屈服彎矩5.3×103kN·m。承臺尺寸為長×寬×高：4.6 m×4.6 m×2.3 m；柱墩高10.9 m。上部結構(即模型頂部質量塊)配重2 000 kN。試驗以1995年日本Kobe地震中的地震記錄作為基底激勵，如圖2所示。

圖1 離心機試驗Fig.1 Centrifuge test setup

圖2 模型輸入地震動Fig.2 Base input motion

1.2 數值模型

本文數值模擬借助有限元計算平臺OpenSees完成。針對上述離心機試驗，考慮到試驗體的對稱性，對稱建立了液化場地群樁-土-上部結構動力相互作用三維有限元數值模型(圖3)，模型長×寬×高：51 m×10.5 m×20.5 m。

圖3 有限元模型Fig.3 Finite element model

模型中，飽和砂土選用土-水完全耦合的六面體單元模擬。該單元是基于Biot動力方程，將飽和砂土模擬為液固兩相介質。單元本質上仍為20節點六面體等參數單元(即為20-8節點單元)，其中，角上的節點包括1個孔壓自由度和3個位移自由度，棱上節點只含3個位移自由度。

砂土本構模型采用Elgamal等[9]提出的多屈服面彈塑性本構模型。該本構模型能夠準確地模擬砂土的液化動力特性。在多屈服面彈塑性本構模型理論框架下，模型采用數量可自定的圓錐形屈服面。這些屈服面通過定義不同的切向剪切模量，表示不同的土體應力-應變關系。此外，通過定義加載-卸載流動法則，該模型可以有效地模擬砂土的剪脹效應。該本構模型已通過離心機試驗結果進行了大量驗證，具體詳見文獻[9]。砂土本構模型參數采用Lu等[10]給出的本構模型建議值，見表1。

表1 砂土的參數Table 1 Model parameters for sand

選用線彈性的梁-柱單元模擬。樁和柱墩材性計算參數與試驗中樁和柱墩的物理特性參數保持一致。樁的尺寸效應，通過在樁單元與土單元直接布置剛度極大的桿單元實現(即剛性連接)，見圖4，其中，桿單元的抗彎剛度設置為樁的10 000倍[13]。樁-土之間接觸面的模擬，借鑒Rahmani和Pak提出的方法[3, 11]，通過在樁-土之間設置零厚度單元實現，如圖4所示。為了提高計算的穩定性與收斂性，模擬承臺仍然采用20-8節點六面體實體單元模擬。不過，計算中約束了其節點上孔壓自由度。上部結構通過在柱墩頂端節點施加集中質量實現模擬。

圖4 樁-土相互作用模擬Fig.4 Modeling of soil-pile interaction

模型中，將基底和側面設置為不透水邊界，保證在該面內孔隙水流速為零[13]。地表處定義為自由面，保證水位線處孔壓為零。側面邊界選用OpenSees中的Shear beam邊界。該邊界可以實現位于同一高度上與振動方向垂直的2個側面邊界產生的水平向位移保持一致，模擬試驗中層狀剪切土箱的土體剪切[14](圖1)。

1.3 數值模型可靠性驗證

圖5～圖7對比了離心機試驗和數值計算的自由場土體的超孔壓和加速度時程、上部結構加速度時程、承臺加速度和樁的彎矩時程。

圖5 超孔壓試驗值與計算值對比Fig.5 Comparison of the experimental and computed results of excess pore pressure

可以看出，試驗和計算獲得的孔壓吻合很好。計算得到的2.9 m埋深處土體在2～4 s內加速度幅值略大于試驗值，但總體上土體加速度峰值的試驗值與計算值相差不大。震動開始時，上部結構加速度的試驗值與計算值吻合度較低，但幅值相差不大；相比而言，震動后期上部結構加速度的試驗值與計算值吻合度較高。在2～4 s和6～8 s時段內，承臺加速度試驗值與計算值的幅值整體變化趨勢稍有差異，但均在合理范圍內。樁1彎矩時程的計算值與試驗值吻合較好，說明該模型可以很好地反應液化場地樁的動力響應規律。

圖6 加速度試驗值與計算值對比Fig.6 Comparison of the experimental and computed results of acceleration

圖7 彎矩試驗值與計算值對比Fig.7 Comparison of the experimental and computed results of bending moment

如上所述，雖然各響應的試驗值與計算值在某個時間段可能存在一定的差異性，但是總體上計算和試驗得到的曲線吻合很高，能夠較真實地重現離心機試驗中群樁和上部結構的動力響應基本規律與典型特征。

2 群樁-土動力相互作用p-y曲線

用于描述樁-土相互作用的p-y曲線中，土反力用p表示，相對位移用y表示。目前，靜力條件下砂土p-y曲線的研究已經比較成熟，特別是API(american petroleum institute)規范中建議的砂土p-y曲線已廣泛應用。但是，動力作用下液化砂土的p-y曲線的研究還不夠成熟，且現有p-y曲線的研究均未考慮群樁效應。群樁效應是指群樁基礎承受荷載后，由于承臺、樁、土的相互作用使其樁側阻力、樁端阻力、沉降等性狀發生變化而與單樁明顯不同，承載力往往不等于各單樁承載力之和這一現象。鑒于此，本節在上述群樁數值模型的基礎上，進行液化場地群樁-土地震相互作用p-y曲線的研究。

2.1 p乘因子法的局限性

API規范中建議的砂土p-y曲線已被廣泛應用于工程設計。目前，關于液化砂土p-y曲線的研究，以Boulanger等[12]提出的“p-乘因子法”(p-multiplier)應用最為廣泛。該方法以API規范中未考慮土體液化效應的砂土p-y曲線為基礎，通過對土反力p進行折減來考慮土體的液化效應。

為了研究“p-乘因子法”在液化場地群樁-土動力相互作用分析中的適用性，根據唐亮[13]的研究結果，需要對比“p-乘因子法”提出的p-y曲線與有限元計算得到的動力p-y曲線骨干線的區別。按照唐亮[13]建議的方法，通過擬合不同幅值下動力p-y曲線滯回圈的頂點，可以獲得其動力p-y曲線的骨干線，這里分別探討土體液化前后“p-乘子法”的有效性。

在上述有限元模型基礎上，選取幅值0.10g、0.15g、0.20g、0.25g、0.30g、0.35g、0.40g和0.50g的正弦波對模型進行激勵。圖8(a)給出了孔壓比為0.3時的動力p-y曲線骨干線繪制過程，即基于計算得到的動力p-y曲線滯回圈，確定其滯回圈頂點，擬合各個頂點得到其骨干線。圖8(b)對比了埋深5.7 m處有限元計算得到的動力p-y曲線和“p-乘因子法”提出的p-y曲線的區別，發現兩者整體趨勢基本相同。但是“p-乘因子法”的土反力稍大于數值計算結果。

圖8 動力p-y曲線骨干線Fig.8 Backbone curves of dynamic p-y

圖9為土體孔壓比為0.6和0.9時的對比情況，發現其規律與孔壓比為0.3時的相同，即兩者變化趨勢相近。但是，“p-乘因子法”計算得到的土反力偏大。因此，“p-乘因子法”可以描述液化前群樁-土相互作用p-y曲線的基本形式，但仍需對其土反力大小做進一步修正。需要說明的是，其他孔壓比及不同埋深處的動力p-y曲線具有相同的規律，這里不再一一列出。圖10對比了液化后土體動力p-y曲線與“p-乘因子法”提出的p-y曲線，發現兩者相差較大。液化后的土體動力p-y曲線骨干線呈現上“凹”形式，這是不排水條件下，循環荷載作用時飽和砂土發生的剪脹現象導致，然而“p-乘因子法”提出的p-y曲線并未考慮這一特征。綜上所述，“p-乘因子法”并不能準確描述液化場地群樁-土地震相互作用p-y曲線特征。鑒于此，本文將基于精細化的有限元計算模型，改進“p-乘因子法”現有局限性，提出更精確的液化場地群樁-土地震相互作用p-y曲線。

圖9 p-乘因子法與動力p-y曲線對比Fig.9 Comparison between p-multiplier and dynamic p-y curve

圖10 p-乘因子法與液化后土體動力p-y曲線對比Fig.10 Comparison between p-multiplier and dynamic p-y curve for the liquefied soil

2.2 改進的p-y曲線

1)土體液化前p-y曲線的修正。

上述研究表明，“p-乘因子法”能夠描述液化前群樁-土相互作用動力p-y曲線骨干線的基本形式，但其土反力值偏大。究其原因在于，該法并不是直接針對群樁基礎提出的，同時未考慮樁間距效應的影響。為此，基于驗證的群樁-土地震相互作用分析的三維有限元模型，以樁間距和樁徑的比值(s/d)為控制指標，對“p-乘因子法”構建的p-y曲線進行修正。

圖11為不同樁間距下(s/d分別等于3～8)計算得到動力p-y曲線與“p-乘因子法”中的p-y曲線的對比。同時，針對“p-乘因子法”中土反力偏大的問題，分別選取不同的樁間距修正系數m對其進行修正。

圖11 不同樁間距下改進的p-y曲線Fig.11 Improved p-y curves under the different pile spaces

由圖可知，修正后的p-y曲線與計算得到的動力p-y曲線滯回圈頂點位置吻合較好。此外，隨著樁間距逐漸增大，修正系數m逐漸減小。當s/d=8時，不需要采取任何修正，“p-乘因子法”也能夠很好地刻畫計算得到的動力p-y曲線滯回圈頂點位置。換言之，當樁間距大于8倍樁徑后，“p-乘因子法”能夠直接用于液化場地群樁-土動力相互作用分析中。對于樁間距小于8倍樁徑的情況，圖12給出了樁間距與修正系數的關系，可見兩者之間存在較好的線性相關性。因此，擬合得到修正系數m的計算公式：

圖12 修正系數m與樁間距的關系Fig.12 Relationship between correction coefficient m and pile space

(1)

式中：m為樁間距修正系數(s/d≤7)，s/d≥8時，取m=1；s為樁間距，m；d為樁徑，m。

據此，基于“p-乘因子法”，結合式(1)，得到了考慮樁間距效應的土體液化前p-y曲線表達式：

p=mmppAPI=m(1-0.9ru)A×

(2)

式中：m為樁間距修正系數(m≤7)，m≥8時，取m=1；mp為孔壓修正系數；ru為孔壓比；pAPI為API公式中的土反力p。

2) 土體液化后p-y曲線的修正。

研究發現，液化后土體動力p-y曲線骨干線呈現出一種“上凹”的形式。因此，為得到液化后土體p-y曲線的表達式，可選用指數函數形式對計算得到的動力p-y曲線滯回圈頂點進行擬合，即液化后土體p-y曲線的統一表達式為：

p=AyB

(3)

式中：p為土反力；y為樁-土相對位移；A、B為曲線形狀控制參數。

以埋深5.7 m處為例，如圖13所示，通過液化后土體動力p-y曲線滯回圈頂點擬合得到其p-y曲線的表達式：

圖13 埋深5.7 m處液化后土體動力p-y曲線Fig.13 Dynamic p-y curve for the liquefied soil at 5.7 m depth

p=1.081y1.01

(4)

埋深5.7 m處的p-y曲線的曲線形狀控制系數A和B的值分別為1.08和1.01。其他埋深處的曲線形狀控制參數如表2所示。

為了得到便于推廣使用的液化后土體p-y曲線計算公式，基于表2中曲線形狀控制系數A、B和土體埋深之間的關系，擬合得到系數A和B的計算公式：

表2 不同埋深處p-y曲線形狀系數A和B

A=-0.01z2+0.17z+0.36

(5)

B=0.02z2-0.35z+2.56

(6)

式中：z為土體埋深，m；A為曲線形狀控制參數；B為曲線形狀控制參數。

通過以上分析，最終建立一組可描述液化后樁-土動力相互作用特性的p-y曲線，該曲線主要由深度決定，可依據式(4)～(6)聯合計算獲得。

2.3 改進的p-y曲線正確性檢驗

Boulanger等[14]基于文克爾地基梁假設和離心機試驗結果，提出了能夠考慮樁-土分離與滑動的力學模型，較為全面地刻畫地震作用下樁-土地震相互作用的物理過程。但是，該模型采用的是API規范推薦的p-y曲線，未考慮群樁效應的影響。本文將提出的改進后的p-y曲線，融入到Boulanger等[14]提出的文克爾地基梁模型中，進行改進后的p-y曲線的正確性檢驗。

針對上述離心機試驗，基于文克爾地基梁法，建立如圖14所示的液化場地群樁-土地震相互作用分析模型。模型中，樁和柱墩采用線性梁-柱單元模擬，樁和柱墩的力學參數和幾何參數與試驗體保持一致?？紤]承臺具有較大的抗彎剛度，不易破壞。因此，選用剛性連接模擬承臺，其抗彎剛度與試驗體保持一致。上部橋梁結構通過在柱墩頂端節點施加集中質量進行模擬。樁-土相互作用通過Boulanger等[14]提出的文克爾地基梁模型模擬，其p-y曲線采用本文提出的改進后的p-y曲線。自由場土體的孔壓和位移時程借助軟件Cyclic 1D進行計算。

圖14 數值模型Fig.14 Numerical modeling

圖15和圖16分別驗證了上部結構加速度和樁1峰值彎矩計算結果的正確性，圖中明顯看到計算值和試驗值吻合較好。因此，采用本文提出的改進后的p-y曲線可以較為精確的描述地震荷載作用下群樁-土-上部結構體系的動力響應。同時，相對有限元法，簡化分析方法避免了土體本構模型的選取問題，且計算時間短、計算參數也容易確定。

圖15 上部結構加速度時程試驗值與計算值對比Fig.15 Computed and experimental acceleration time histories of superstructure

圖16 樁1峰值彎矩試驗值與計算值對比Fig.16 Computed and experimental peak bending moment of pile 1

3 結論

1) 通過有限元數值模型，計算獲得以不同幅值的正弦波為模型激勵的群樁-土動力相互作用p-y曲線，發現“p-乘因子法”高估了土體液化前的樁上土反力，忽略了液化后砂土的剪脹特性。

2) 通過對比不同樁間距下，“p-乘因子法”得到的p-y曲線和動力p-y曲線的滯回圈骨干線，分析得到樁間距小于8倍樁徑時，樁間距是影響“p-乘因子法”準確性的重要因素。

3) 基于傳統的“p-乘因子法”中p-y曲線，以樁間距和樁徑的比值(s/d)為控制指標，提出考慮樁間距效應的液化場地群樁-土動力相互作用p-y曲線，提高了p-y曲線法對于群樁分析的準確性。