液化場地橋梁足尺樁抗震簡化分析方法①

張效禹， 唐　亮,3， 凌賢長， 李　惠

(1. 哈爾濱工業大學土木工程學院，黑龍江 哈爾濱 150090；

2. 結構工程災變與控制教育部重點實驗室 哈爾濱工業大學，黑龍江 哈爾濱 150090；

3.地質災害防治與地質環境保護國家重點實驗室 成都理工大學，四川 成都 610059)

?

液化場地橋梁足尺樁抗震簡化分析方法①

張效禹1,2， 唐亮1,2,3， 凌賢長1,2， 李惠1，2

(1. 哈爾濱工業大學土木工程學院，黑龍江 哈爾濱 150090；

2. 結構工程災變與控制教育部重點實驗室 哈爾濱工業大學，黑龍江 哈爾濱 150090；

3.地質災害防治與地質環境保護國家重點實驗室 成都理工大學，四川 成都 610059)

摘要：采用國際VELACS項目中離心機試驗標定的內華達砂的動力計算參數，建立液化場地足尺樁-土動力相互作用分析的三維有限元模型；獲得不同幅值的正弦波作用下樁-土動力相互作用的p-y曲線，修正并發展一種可用于液化場地樁-土動力相互作用分析的宏單元模型，并基于非線性文克爾地基梁模型建立橋梁足尺樁抗震分析的數值模型與簡化方法，通過有限元分析結果驗證該簡化方法的正確性。

關鍵詞：宏單元模型； 簡化方法； 足尺樁-土動力相互作用； 非線性有限元分析； 液化場地

0引言

強震作用下場地樁-土相互作用的分析是進行液化場地橋梁樁基結構抗震分析的基本依據。雖然模型試驗(離心機試驗和振動臺試驗)能夠有效地進行樁-土動力相互作用分析，但存在耗時長、費錢且考慮的影響因素有限等不足。目前數值方法已成為橋梁樁基抗震設計的重要手段之一，它可以避免試驗方法的諸多缺點，并且較好地模擬地震下樁的動力特性。然而數值方法也具有建模途徑較復雜、計算參數不易選取且計算機時較長等缺點，在實際工程應用中受到很大限制[1-2]。因此發展一種高效且可靠的液化場地橋梁樁基抗震分析方法，即動力特性的分析方法，具有重要意義。

非線性文克爾地基梁法作為一種簡化方法，可以很好地模擬地震下土體的強非線性、輻射阻尼及樁-土界面處分離和滑動等特性，已在工程實踐中得到廣泛應用。然而將非線性文克爾地基梁法應用到液化場地樁基抗震分析中的研究尚不成熟，現有研究工作更多是通過一維自由場地分析獲得砂層孔壓比衰減修正，得到彈簧剛度，如Kagawa等[3]、Liyanapathirana等[4]及Sarkar等[5]基于文克爾地基梁模型建立的液化場地橋梁樁基抗震分析的簡化方法。

地震作用下樁-土動力相互作用是一個極其復雜的過程，對樁-土動力相互作用實施精細化模擬的關鍵在于如何準確刻畫其動力特性。通常將土體按距離樁的遠近分為近場土體、遠場土體和自由場土體。地震作用下近場土體表現出明顯的非線性特性，且在強地震作用下樁和土體之間還可能出現裂縫及樁土相對滑動的現象。遠場土體距離樁較遠，受到樁的作用較小，土體表現出近似彈性的行為，對遠場土體模擬的關鍵在于描述其輻射阻尼效應。自由場土體幾乎不受樁的影響，但其濾波作用會使樁受到不同的激勵作用。目前國內外學者提出了多種宏單元模型：如Matlock模型[6]、Naggar模型[7]及Boulanger模型[8]等，試圖準確表述模擬地震作用下樁-土相互作用的動力特性。其中Boulanger提出的宏單元模型由彈性單元、塑性單元和裂縫單元串聯組成，能夠更全面且合理地模擬樁-土相互作用的物理過程。目前液化場地樁基的抗震研究工作多集中于小比例模型試驗及相應的理論分析。本文擬通過修正并發展Boulanger模型，提出一種可進行液化場地橋梁足尺樁-土動力相互作用分析的簡化方法。

1液化場地橋梁足尺樁-土動力相互作用數值模擬與分析

采用唐亮等[9-12]提出的液化場地樁-土動力相互作用分析的基本數值建模途徑與計算方法，采用國際VELACS(Verification of liquefaction analysis using centrifuge studies)項目標定的相對密度40%的Nevada(內華達)砂計算參數，建立液化場地橋梁足尺樁基地震反應分析三維有限元模型(圖1)。內華達砂的計算參數見表1[11-12]。

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

計算參數取值參考剪切模量G0(kPa,Pr=80kPa)7.85×104參考體積模量B0(kPa,Pr=80kPa)221850峰值剪應變gmax(Pr=80kPa)0.1密度/(kg·m-3)1973內摩擦角?/(°)31.4參考平均有效圍壓Pr/kPa80收縮系數c20.3膨脹系數d10.4膨脹系數d22液化系數y1/kPa10液化系數y20.01液化系數y33滲透系數/(m·s-1)6.6×10-5壓力相關參數np0.5相位轉換角?pt/(°)26.5

模型中，樁徑0.6 m，樁長12 m，彈性模量E=3.0×107kPa，樁頂施加240 kg重的質量點以模擬上部結構在地震中的慣性效應。樁采用線彈性梁-柱單元模擬，土層為10 m厚的飽和內華達砂，水位線位于地表處。模型長30.6 m×寬15 m×高10 m。樁與土采用剛性連接桿模擬。砂土選用的本構模型、邊界條件和計算收斂準則等其他建模的技術細節詳見文獻[9-10]。

2動力p-y曲線

2.1不同加載幅值下的動力p-y曲線

采用上述有限元模型，模型基底輸入不同幅值的1 Hz正弦波，獲得液化場地樁-土動力相互作用的p-y曲線，其中p為樁和土之間的相互作用力(簡稱“土反力”)，y為樁和土之間的相對位移。圖2為不同幅值(0.05g、0.10g、0.15g和0.25g)的1 Hz正弦波輸入下埋深4 m處樁與土相互作用的動力p-y曲線。圖3為不同幅值的正弦波下埋深6 m處砂層超孔壓時程。

圖2　加載幅值對埋深4 m處土體動力p-y曲線的影響Fig.2　Effect of loading amplitudes on dynamic p-y curves at the depth of 4 m

圖3　加載幅值對埋深6 m處砂層超孔壓時程影響Fig.3　Effect of loading amplitudes on excess pore pressure　　　 time-histories of sand at the depth of 6 m

由圖2、圖3可知，加載幅值較小時(0.05g)，土體的超孔壓較小，p-y曲線滯回圈幾乎呈線性；隨著加載幅值的增加，土體發生液化(孔壓比ru=1)，滯回圈面積逐漸越大；當幅值由0.15g增至0.25g時滯回圈面積并未發生明顯增大，這應該主要是由于土反力p接近極限土反力pult所致。此外，隨著加載幅值的增加，土體超孔壓增長率逐漸變大。

2.2動力p-y曲線的骨干曲線

借鑒唐亮[1]提出的方法，引入構造土的動應力-動應變關系曲線骨干線的基本思想，將不同加載幅值下的動力p-y曲線頂點擬合得到其骨干曲線，據此發展液化場地樁-土動力相互作用分析的簡化方法。由上述可知，土體孔壓比對動力p-y曲線影響顯著，土體單位樁長的極限承載力pult在土體受到震動過程中隨孔壓比的變化而改變。因此需要構建不同孔壓比下動力p-y曲線的骨干曲線。將孔壓比按每0.2一級劃分六段：0～0.2、0.2～0.4、0.4～0.6、0.6～0.8、0.8～1.0和1.0。

動力p-y曲線的骨干曲線構建過程為：

(1) 孔壓比ru=1.0

根據上述研究可知，孔壓比ru=1.0時動力p-y曲線滯回圈表現出明顯的滯后性和變形累積等特性，滯回圈頂點的位置很難確定。因此選取p-y曲線滯回圈每個加載周期的最大土反力pmax和最大樁土相對位移ymax點作為其頂點，以0.10g 1Hz正弦波輸入下埋深4m處土體p-y曲線為例(圖4)，選取點(p1，y1)、點(p2，y2)作為動力p-y曲線滯回圈頂點。將美國石油工程協會(AmericanPetraleumInstitute，API)規范[13]中推薦的標準砂土p-y曲線進行折減，擬合埋深4m處的骨干曲線[圖5(a)]。折減系數Cz=0.16時，得到的曲線與滯回圈頂點的排列形式吻合較好，即埋深4m處土體的動力p-y曲線的骨干曲線為：

(1)

圖4　0.10g 1 Hz正弦波輸入下埋深4 m處　　　土體動力p-y曲線時程Fig.4　Dynamic p-y time history at the depth of 4 m under　　　 the input of sine wave with amplitude of 0.1g and　　　 frequency of 1 Hz

圖5　埋深4 m和8 m處土體動力p-y曲線的　　　骨干曲線(ru=1.0)Fig.5　Backbone curve of the dynamic p-y curve at　　　 the depth of 4 m and 8 m (ru=1.0)

隨后選取埋深8 m處土體的p-y曲線，驗證折減系數Cz=1.6是否適用于其他埋深處的骨干曲線[圖5(b)]。由圖5(b)可以看出，雖然土反力p沒有達到極限土反力pult，但是折減得到的骨干曲線初始部分與滯回圈頂點的排列形式吻合較好，因此可假定折減系數Cz=0.16同樣適用于埋深8 m處土體。限于篇幅，不再對其他埋深處逐個進行檢驗(詳見參考文獻[2])。因此孔壓比ru=1.0時的動力p-y曲線骨干曲線可以通過API規范中推薦的砂土p-y曲線乘以折減系數0.16得到，即

(2)

式中：pu為砂土的極限承載力(kN/m)，取式(3)和式(4)中較小者；A為荷載類型系數，可根據API規范確定[13]；H為計算點深度(m)。按楔體失效理論計算砂土極限承載力pu為：

(3)

按流動失效理論計算砂土極限抗力pu為：

(4)

式中:γ為土體有效容重(N/m3)；D為樁徑(m)；C1、C2、C3為計算參數，根據API規范確定[13]。

(2) 孔壓比ru<1.0

上述分析可知，當孔壓比ru<1.0時，土體動力p-y曲線滯回圈面積很小，加載幅值增至0.25g時，土反力p與極限土反力pult仍然相距甚遠，無法擬合獲得骨干曲線。然而根據Liu等[14]提出的Cu因子法可知，強度衰減乘因子Cu隨孔壓比改變呈線性變化(Cu=1-ru)。為構建出統一的動力p-y曲線的骨干曲線，以Cu因子法為基礎，假定孔壓比ru<1.0時，折減系數隨孔壓比改變呈線性變化，即折減系數Cz=1-0.84ru，孔壓比ru<1.0時的動力p-y曲線的骨干曲線式為：

(5)

以埋深4 m處孔壓比ru=0.2和ru=0.6時的動力p-y曲線滯回圈頂點，驗證式(5)的有效性(圖6)。可見，雖然土反力p和樁土相對位移y均很小，但基本上可以描述出動力p-y曲線的骨干曲線初始剛度。因此式(5)可適用于孔壓比ru<1.0時的動力p-y曲線的骨干曲線。

圖6　埋深4 m處土體動力p-y曲線的　　　骨干曲線驗證Fig.6　Verification of backbone curve of the dynamic　　　 p-y curve at the depth of 4 m

綜上，土體動力p-y曲線的骨干曲線表達式為：

(6)

3修正的宏單元模型

Boulanger等[8，15]提出的宏單元模型用于刻畫樁-土動力相互作用的物理過程時，需要確定四個計算參數：模型的極限承載力pult、荷載p=0.5pult時的位移y50、最大摩擦力與pult之比Cd和阻尼系數c。通過對這些參數進行修正，將宏單元模型應用于液化場地的樁-土動力相互作用分析中。

3.1基本元件

Boulanger的宏單元模型由彈性單元(p-ye)、塑性單元(p-yp)和裂縫單元(p-yg)串聯組成(圖7)。

圖7　宏單元模型Fig.7　Macro-element model

彈性單元由彈性元件和阻尼元件并聯而成。彈性元件的力p與位移關系為

(7)

式中：Ke為彈性元件的切線模量；ye為彈性元件的位移。

塑性單元表示土體的塑性變形。塑性單元和彈性單元串聯，二者具有相等的力p。塑性元件中p-yp曲線見圖8 (a)，力p按式(8)確定：

(8)

塑性單元的切線模量Kp定義為：

(9)

Matlock(1970)建議，軟質黏土取C=10、n=5和Cr=0.35。API規范建議砂土取C=0.5、n=2和Cr=0.2。本文選取API規范給出的建議值。

裂縫單元由一個非線性的閉合元件(pc-yg)和一個非線性的摩擦元件(pd-yg)并聯組成。裂縫單元的力p=pc+pd。閉合元件控制著樁-土界面處裂縫的張開和閉合，作用機理與Matlock等[6]提出的黏土p-y屬性的裂縫一致。采用摩擦元件描述樁土相對滑動過程中樁受到的側向摩擦力。閉合元件pc-yg曲線見圖8 (b)，摩擦元件pd-yg曲線見圖8 (c)。閉合元件的力pc為：

(10)

摩擦元件的力pd為：

(11)

圖8　宏單元模型荷載-位移曲線Fig.8　Load-displacement curves of macro-element model

(12)

(13)

裂縫單元的切線模量按式(14)確定:

(14)

3.2組合模式

將上述三種單元串聯得到刻畫樁-土動力相互作用的宏單元模型[圖8 (d)]。模型的位移見式(15)，模型切線模量見式(16)。

(15)

(16)

模型中各元件相互協調發揮作用：當土體在彈性范圍內變形，只有彈性元件起作用。隨荷載增加，土體進入塑性，彈性元件和塑性元件共同作用。荷載卸載時，樁-土界面處出現裂縫，裂縫單元、彈性單元和塑性單元共同工作。當樁向裂隙移動時，只有摩擦元件發揮作用。

3.3模型參數

通過前述建立的液化場地樁-土相互作用的動力p-y曲線骨干曲線表達式[見式(6)]確定pult和y50，骨干曲線中土反力p的極值即為模型的極限承載力pult，當p=0.5pult時對應的位移即為y50。參照Liyanapathirana[16]的方法，取阻尼系數c=ρvS(vS為剪切波速) 。

最大摩擦力與pult的比值Cd參照Brandenberg等[15]給出的建議值取Cd=1.0。

4簡化分析方法與數值模型

4.1數值模型

與三維有限元模型保持一致，基于文克爾地基梁模型，建立樁的簡化分析數值模型(圖9)。樁采用線彈性梁-柱單元模擬，分24個單元。樁頂節點模擬為集中質量點以模擬上部結構。樁與土的相互作用采用上述動力宏單元模型代替。將動力p-y模型采用2節點(從屬彈簧節點和固定彈簧節點)構成的零長度單元表示，將其從屬彈簧節點與樁節點的水平方向自由度捆綁在一起。在固定的彈簧節點處輸入土體自由場的位移時程作為邊界激勵條件，完成樁-土動力相互作用的模擬與分析。

圖9　簡化分析的數值模型示意圖Fig.9　Numerical model for simplified method

4.2外部輸入激勵

簡化分析方法為兩部分：(1)進行自由場地的地震反應分析，獲得土體的位移和孔壓時程，作為簡化模型的外部激勵；(2)輸入外部激勵，進行液化場地樁-土動力相互作用的分析。選取El Centro波(圖10)進行自由場土體的動力反應分析，得到其位移時程和孔壓比時程(圖11)。采用簡化數值模型，輸入自由場土體的位移時程作為外部激勵，并通過考慮孔壓比改變pult和y50，進行樁-土動力相互作用分析。

圖10　模型基底輸入運動Fig.10　Input at model base

4.3簡化方法的正確性檢驗

針對簡化數值模型，按照上述參數取值的方法確定各個參數的值，在每一個動力p-y模型的固定彈簧節點輸入與其相同埋深處自由場地土體的位移時程，以El Centro波輸入下樁的動力特性進行簡化分析方法的正確性檢驗。

圖11　自由場土體位移時程與孔壓比時程曲線Fig.11　Displacement time histories and pore pressure　　　 ratio time histories of free field soil

樁的位移是體現樁動力特性的重要指標，它可以通過微分得到樁的加速度。因此樁的位移基本上能體現出樁的動力反應規律。圖12為簡化方法和有限元法計算得到的不同位置樁的位移時程對比情況。由圖可見，兩種方法得到的位移時程曲線整體波動趨勢基本相同，振動前3 s結果吻合較好，3 s以后簡化方法計算得到的結果稍大于有限元法的結果。總之，簡化方法基本上能夠較好地模擬出更為精細化的有限元法，且其計算時間遠小于有限元法的計算機時，所需確定的計算參數也較少，相比而言更適合工程實踐的應用。

圖12　簡化方法和有限元法計算結果Fig.12　Results of simplified method and finite element method

5結論

基于液化場地樁-土動力相互作用分析的數值建模途徑與計算方法，采用離心機試驗標定的內華達砂的動力計算參數，建立液化場地橋梁足尺樁-土動力相互作用三維有限元模型，得到：

(1) 分析不同加載幅值對土體動力p-y曲線影響。結果表明，加載幅值較小時，土體不發生液化，p-y曲線滯回圈面積很小；加載幅值越大，其滯回圈面積越大；當幅值大到一定程度，滯回圈面積不再隨幅值的增加發生明顯變化，顯示土體達到了極限狀態。

(2) 通過對API規范中推薦的砂土p-y曲線的計算公式進行修正，得到不同孔壓比下液化砂土動力p-y曲線的骨干曲線表達式。

(3) 通過新的液化砂土動力p-y曲線的骨干曲線表達式，改進Boulanger提出的宏單元模型，并將其應用于液化場地的樁-土動力相互作用分析中。

(4) 基于改進的Boulanger模型，提出液化場地樁-土動力相互作用簡化分析方法，并驗證簡化方法的正確性。

參考文獻(References)

[1]唐亮.液化場地樁-土動力相互作用p-y曲線模型研究[D].哈爾濱: 哈爾濱工業大學, 2010.

TANG Liang.p-yModel of Dynamic Pile-soil Interaction in Liquefying Ground[D].Harbin:Harbin Institute of Technology,2010.(in Chinese)

[2]張效禹.液化場地橋梁足尺樁-土動力相互作用簡化分析方法[D].哈爾濱:哈爾濱工業大學,2011.

ZHANG Xiao-yu.Simplified Method of Dynamic Full-scale Pile-soil Interaction in Liquefying Ground[D].Harbin:Harbin Lnstitute of Technology,2011. (in Chinese)

[3]Kagawa T,Kraft L M.Seismicp-yResponses of Flexible Piles[J].Journal of the Geotechnical Engineering Division,1980,106(GT84):899-918.

[4]Liyanapathirana D S,Poulos H G.Pseudostatic Approach for Seismic Analysis of Piles in Liquefying Soil[J].Journal of Geotechnical and Geoenvironmental Engineering,2005,131(12):1480-1487.

[5]Sarkar R,Maheshwari B K.Influence of Soil Nonlinearity and Liquefaction on Dynamic Response of Pile Groups[C]//The 14thWorld Conference on Earthquake Engineering.Beijing,China:[s.n.],2008.

[6]Matlock H,Foo S H C.Simulation of Lateral Pile Behavior Under Earthquake Motion[C]//Speciality Conference on Earthquake Engineering and Soil Dynamics.1978,2: 600-619.

[7]Naggar M H,Novak M.Nonlinear Analysis for Dynamic Lateral Pile Response[J].Soil Dynamics and Earthquake Engineering,1996,15(4):233-244.

[8]Boulanger R W,Curras J,Kutter B L,et al.Seismic Soil-pile-structure Interaction Experiments and Analyses[J].Journal of Geotechnical and Geoenvironmental Engineering,1999,125(9):750-759.

[9]唐亮,凌賢長,徐鵬舉,等.液化場地樁-土地震相互作用振動臺試驗數值模擬[J].土木工程學報,2012,45(增刊1):302-306, 311.

TANG Liang,LING Xian-zhang,XU Peng-Ju,et al.Numerical Simulation of Shaking Table Test for Seismic Soil-pile Interaction in Liquefying Ground[J].China Civil Engineering Journal,2012,45(Supp1):302-306,311.(in Chinese)

[10]唐亮,凌賢長,艾哈邁德·艾格瑪.液化側向流動場地樁基動力反應振動臺試驗三維有限元數值模擬方法[J].土木工程學報,2013,46(增刊1):180-184.

TANG Liang,LING Xian-zhang,Elgamal A.Three-dimensional Finite Element Analysis of Shake-table Test for Dynamic Pile Behavior in Liquefaction-induced Lateral Spreading Ground[J].China Civil Engineering Journal,2013,46(Supp1):180-184.(in Chinese)

[11]Mazzoni S,Mckenna F,Scott M H,et al.Opensees Command Language Manual[R].Berkley:Berkley University of California,2007.

[12]Elgamal A,Lu J,Forcellini D.Mitigation of Liquefaction-induced Lateral Deformation in a Sloping Stratum:three-dimensional Numerical Simulation[J].Journal of Geotechnical and Geoenvironmental Engineering,2009,135(11):1672-1682.

[13]American Petraleum Institute.Recommended Practice for Planning,Designing and Constructing Fixed Offshore Platforms[M].API Recommended Practice 2A-WSD (RP2A-WSD),2000.

[14]Liu L,Dobry R.Effect of Liquefaction on Lateral Response of Piles by Centrifuge Model Tests[J].National Center for Earthquake Engineering Research Nceer,1995, 9(1):7-11.

[15]Brandenberg S J,Zhao M,Boulanger R W,et al.p-yPlasticity Model for Nonlinear Dynamic Analysis of Piles in Liquefiable Soil[J].Journal of Geotechnical and Geoenvironmental Engineering,2012,139(8):1262-1274.

[16]Liyanapathirana D S,Poulos H G.Sesmic Lateral Response of Piles in Liqufying Soil[J].Journal of Geotechnical and Geo-Environmental Engineering,2005,1(131):1466-1479.

Simplified Seismic Analysis Method for Bridge

Full-scale Piles in Liquefied Ground

ZHANG Xiao-yu1, 2, TANG Liang1, 2, 3, LING Xian-zhang1,2， LI Hui1,2

(1.SchoolofCivilEngineering,HarbinInstituteofTechnology,Harbin150090,Heilongjiang，China;

2.KeyLaboratoryofStructuresDynamicBehaviorandControl,MinistryofEducation,HarbinInstituteofTechnology,

Harbin150090,Heilongjiang，China3.StateKeyLaboratoryofGeohazardPreventionand

GeoenvironmentProtection,ChengduUniversityofTechnology,Chengdu610059,Sichuan，China)

Abstract：Most research in the seismic analysis of pile foundations in liquefied ground is based on small-scale model testing and theoretical analysis. This paper presents the development of a three-dimensional finite element model for dynamic pile-soil interaction using the dynamic calculation parameters of Nevada sand calibrated by conducting the centrifuge test outlined in the Verification of Liquefaction Analysis by Centrifuge Studies (VELACS) project.On the basis of this numerical model, the dynamicp-ycurves under various shaking amplitudes are studied, and a macro-element model for analysis of dynamic pile-soil interaction in liquefied ground is developed. In addition, a simplified method for seismic analysis of the piles on the basis of the Winkler model is presented and verified by the results of numerical simulation.

Key words：macro-element model; simplified method; dynamic full-scale pile-soil interaction; nonlinear finite element analysis; liquefied ground

DOI:10.3969/j.issn.1000-0844.2015.04.0944

中圖分類號:TU43

文獻標志碼：A

文章編號:1000-0844(2015)04-0944-08

作者簡介：張效禹(1991-)，男(漢族)，博士研究士，研究方向：巖土工程防災減災。E-mail:zxy_hit@163.com。

基金項目：國家自然基金項目(51578195，51378161)；黑龍江省應用技術研究與開發項目(GZ13A009)；國家重點基礎研究發展973計劃項目(2012CB026104)；地質災害防治與地質環境保護國家重點實驗室項目(SKLGP2013K011)

收稿日期：①2014-08-20