基于PM4Sand塑性模型的土體液化側移計算分析

楊彥鑫，陳文龍，馬建林，林大富

(西南交通大學 土木工程學院，四川 成都 610031)

在地震作用下，土體液化側移表現為液化土體產生大變形導致地表開裂及地表土層沿液化土層發生滑動而產生側移這2種形式[1]。土體液化會加劇地震造成的破壞，如2011新西蘭基督城地震、2008年汶川地震等[2-3]。

研究土體液化側移的方法大致分為簡化分析法、經驗公式法、模糊數學法、模型試驗法和數值模擬法。文獻[4]基于極限平衡理論，提出以屈服加速度為基礎的震后永久性位移計算公式。簡化分析方法主要是根據Newmark滑塊法以及改進后的Newmark滑塊法計算土體液化側移。利用數理統計方法，大量震后調查數據被用于研究土體液化側移。其中，文獻[5-6]采用多元線性回歸法(MLR，Multiple Linear Regression)提出的經驗公式應用最為廣泛。經驗計算公式包含了飽和砂土液化的主要影響因素，在一定程度上能夠有效地預測或計算土體液化側移，但不能體現土體內部發生液化側移的機理，因此經驗公式法的應用存在一定的局限性。鑒于飽和砂土液化過程的復雜性和不確定性，文獻[7-8]將神經網絡模型應用于土體液化側移計算中，提出了相應的經驗計算公式。文獻[9]假定液化土層和非液化土層之間無摩擦力，采用彈性梁分析方法計算地震作用下的地表位移。在數值計算方面，文獻[10]為描述土體液化特性建立了一個多屈服面塑性本構模型，并結合室內試驗和離心機試驗驗證了該本構模型的適用性。為了定量描述飽和砂土液化后的大變形，文獻[11]建立了彈塑性循環本構模型，并通過試驗證明了其有效性。文獻[12-13]通過開源有限元軟件OpenSees程序模擬飽和砂土液化，并討論了影響液化的主要因素。

數值計算方法可綜合考慮地震輸入特性、地形地質條件等，可作為一種預測砂土液化側移的有效手段。本文基于有限差分程序FLAC 2D，引入邊界面塑性PM4Sand砂土模型，對已發生液化側移場地進行計算，并與實測值對比，驗證該模型計算液化場地側移值的適用性。

1 PM4Sand模型

Dafalias等[14]于2004年提出DM04模型，該模型是依據土的臨界狀態理論、遵循應力控制原則的砂土塑性模型。DM04模型能夠準確地模擬土體單向受荷時的力學行為。當土體受循環荷載作用采用DM04模型模擬時，會出現土體達到某一階段后，其強度不再隨著循環荷載作用繼續軟化的現象。這與實際情況不符，因此采用DM04模型模擬土體在動力荷載作用下的受力行為存在缺陷[15]。2010年Boulanger[16]首次提出PM4Sand V1模型并應用于FLAC中，該模型是在DM04模型的基礎上進行了改進，以彌補DM04模型模擬液化引起的地表大位移的不足。PM4Sand V1模型做了如下改進：①定義塑性體積應變控制的細觀組構的形成與破壞模式，添加組構演化和累積組構形成變量，建立了塑性模量與組構的關系。②考慮了組構演化對砂土剪脹性的影響，因此更加有效地控制體縮及體脹時的砂土剪脹性變化規律，使其與Bolton剪脹關系式一致。③提出考慮應力比和組構歷史的彈性模量關系式，重新定義了初始反向應力歷史(荷載效應)記錄，采用相對狀態參數指數(Relative State Parameter Index)表達臨界狀態，解決了屈服面和剪脹面的形狀函數以適應平面應變問題。

PM4Sand模型中有6項主要參數，21項次要參數。其中3項主要參數需要輸入，其余參數依據模型手冊要求設置為默認值。將實測標準貫入錘擊數修正到有效覆蓋應力為100 kPa，落錘能量比為60%的修正標準貫入錘擊數，記為(N1)60。根據(N1)60可確定相對密實度Dr、剪切模量系數G0的計算式為

式中：(N1)60需要根據土體中細粒含量對其進行修正[17-18]；Cd為考慮粒徑效應(最大孔隙比與最小孔隙比差)的參數，一般取建議值46。

體縮率hp0在計算時用于調整塑性模量與彈性模量比，當相對密度為38%，55%，78%時，hp0可分別取0.32，0.40，0.50[18]。hp0需根據動力觸探值試算獲得，對應為單元體在15次循環荷載作用下對應的目標循環應力比。

2 算例

2.1 背景簡介

金銀島(Treasure Island)位于美國舊金山地區，建造于20世紀20年代，連接舊金山和奧克蘭市，并毗鄰Yerba Buena島，是一座在砂層和舊金山海灣泥之上用水力疏浚和抓斗挖泥式土建造的人工島嶼，并在砂石填料層上建造堆石圍堰作為支擋結構。

圖1 金銀島典型剖面(單位：m)

金銀島的地質條件與舊金山灣區地質環境類似，從上至下分為砂土填料、天然灘砂、全新世灣泥以及更新世灣泥4層，表層為松散砂質填土。其中填土的靜力觸探值約為1～5 MPa，松散砂質填土容易發生液化。在1989年的Loma Prieta地震中，金銀島距離地震破裂帶北部最近為80 km。在地震作用下觀測到金銀島多處發生液化、砂涌和地表沉降現象。其中最大的側移出現在島的東側，側移值為30.5 cm。金銀島地下勘察資料完善，震后液化側移有記錄數據，因此本文選取該場地作為研究對象[19]。為將PM4Sand應用于側移計算,選取如圖1所示的典型剖面作為分析對象。該剖面位于金銀島北街和第三街的交界處，根據現場對基準點的測量，可以確定該處的水平側移為25.4 cm。該地層截面主要由3部分組成，自上而下分別為砂土填料、全新世灣泥和更新世灣泥。

2.2 土體參數

根據文獻[20-21]獲得金銀島從地表至基巖的地質勘察資料，具體信息見表1。將表1的剖面作反卷積計算以獲得分析所用的地震輸入。圖2為地表至基巖的剪切波波速分布?？芍?，地表下0.9 m處的砂性填土的剪切波波速平均值為165 m /s，容易液化。

表1 金銀島土體基本參數

圖2 剪切波波速分布

2.3 地震輸入參數

將金銀島記錄的實測地震波作為出露地震波，利用Shake 2000軟件進行反卷積計算，得到加載至FLAC模型底部的地震輸入，分別記為TRI000,TRI090地震波，見圖3。

圖3 地震波

3 建模與計算結果分析

利用FLAC 2D有限差分程序，建立土體剖面模型。先將所有材料參數設置為摩爾-庫倫模型(參數取值見表2)，設置大應變計算模式，設置較大的黏聚力生成初始應力。設置地下水位，各土層的孔隙率和滲透率，水的密度和體積模量，進行滲流計算，生成孔隙水壓力，使系統達到平衡。系統平衡以后，將黏聚力和內摩擦角參數設置為真實值，可液化材料(水位以下的吹填土)設置為PM4Sand模型(參數取值見表3)，其余土體仍然采用摩爾-庫倫模型，并對模型中的非液化土體設置模量衰減曲線和阻尼比曲線。

表2 土體參數

表3 PM4Sand模型參數

注: (N1)60,cs為經細粒含量修正后，標準大氣壓下的純凈砂動力觸探值。

模型中土體的動剪切模量G和阻尼比曲線見圖4。利用FLAC中的Hardin模型表達模量衰減曲線和阻尼比曲線,即將模型用Hardin參數表示，將Shake 2000中得到的輸入地震波加載至模型底部。由圖4可知，在Hardin模型中，隨著土體應變增大，動剪切模量減小，阻尼比增大。

圖4 土體的動剪切模量曲線和阻尼比曲線

設置靜止邊界和自由場邊界減少模型邊界上的地震波反射，取中心頻率f為15 Hz，設置瑞利阻尼0.02吸收地震波的高頻部分。根據文獻[22]模型中砂土的動力觸探值(N1)60為10.0，砂土的細粒含量為10%，根據文獻[23]，將(N1)60修正為純砂動力觸探值(N1)60,cs=11.0，計算得到PM4Sand參數。

為了判斷模型是否在地震動力作用下達到平衡，需計算系統最大不平衡力(見圖5)，即所有節點的外力與內力之差的最大值?？芍?，作用在系統上的最大不平衡力在40 s后趨于穩定，表面系統達到平衡狀態。

圖5 系統不平衡力

圖6 孔隙水壓力時程曲線

有效應力、孔隙水壓力、有效孔壓比均可表征系統是否發生液化。這里利用砂土中的孔隙水壓力的變化表征砂土是否發生液化。砂土層孔隙水壓力隨時間變化曲線見圖6。可知，初始時刻可液化層在飽和狀態下的孔隙水壓約為60 kPa。監測砂性填土底部孔隙水壓力變化，在2條地震波分別作用下，孔隙水壓發展并達到穩定，在TRI000地震波作用14.8 s時場地發生液化，孔隙水壓力趨于穩定；在TRI090地震波作用13.0 s 時場地發生液化，之后孔隙水壓力趨于穩定?？紫端畨毫Φ淖兓砻髟诘卣鹱饔孟聢龅鼐l生了液化。

不同地震波時地表水平加速度時程曲線見圖7?？芍斎爰铀俣葧r程被場地放大。在TRI000地震波作用下，地震波的峰值加速度由輸入的0.24g放大至0.34g，在TRI090地震波作用下，地震波的峰值加速度由輸入的0.41g放大至0.43g。

圖7 不同地震波時地表水平加速度時程曲線

峰值加速度隨著深度的變化見圖8?？芍?，當土層剪切波波速較大，接近基巖的剪切波波速時，峰值加速度將接近地震輸入的峰值加速度。而剪切波波速較小時，對地震有一定的放大作用。最大的峰值加速度出現在地表下約6.5 m處，最大的峰值加速度分別為0.512g，0.695g，對應于砂性填土。

圖8 峰值加速度隨著深度的變化

計算求得輸入地震波和地表地震波的地震反應譜，地震反應譜受到阻尼比和地震動的影響，隨著阻尼比的增大地震反應譜值變小，而地震動的振幅對反應譜值的大小影響較大。設置阻尼比為0.5%，獲得最大加速度和其對應的自振周期。圖9所示為對應2條輸入地震波的擬加速度反應譜?？芍搱龅氐淖畲蠹铀俣葘谥芷?.5～1.5 s。

圖9 不同地震波時擬加速度反應譜

圖10 不同地震波時場地最大剪應變云圖

地震結束后，場地的最大剪應變云圖見圖10。在TRI000和TRI090地震波作用下，土層最大剪應變均發生在堆石斜坡下方的砂土填料處?？芍?，場地處于塑性狀態且砂土填土區域的應變為2%～3%，砂土填土已發生液化。

在土層最左側自由面取監測點，監測點處水平位移體現該土層液化側移量。土層在2條地震波分別作用下自由面側移量時程曲線見圖11?？芍琓RI000和TRI090地震波對土層作用使自由面發生最大側移分別為446，559 mm。文獻[22]中給出金銀島在Loma Prieta地震用下液化側移值約為254 mm。由計算結果可知，預測值較液化實測值大，且地震波的方向與場地自由面的夾角會造成2條地震波所得的位移有所差異，故TRI000地震波所得位移較小，TRI090地震波所得位移較大。其中計算值分別是觀測值的1.76倍和2.2倍，預測值具有一定的安全儲備。

圖11 自由面側移時程曲線

土體的初始剪切模量會對土體變形產生影響，較高的剪切模量會導致預測值偏?。煌馏w的動剪切模量衰減曲線和阻尼比曲線難以確定，且不同的土體動剪切模量衰減曲線和阻尼比曲線會對計算結果產生一定的影響；場地土體參數難以確定，過高或者過低的參數估計會對預測結果產生影響，需要進行室內試驗對土層參數加以確定。地震輸入會對液化側移預測產生影響，該場地并沒有實測基巖地震波，本算例的地震加速度輸入是由反卷積計算所得，地震輸入會對預測結果造成一定的影響。

PM4Sand作為一種土體的塑性模型，依據可液化土動力觸探試驗值以及砂土的粒徑含量確定模型所需要的主要參數，在參數確定上簡單，同時可以有效模擬孔隙水壓力的發展與變化。利用有限差分程序獲得可液化場地的液化側移。在實際計算過程中，可同時與簡化計算方法、經驗公式法比較研究，獲得較為合理的預測值。

4 結論

應用有限差分程序，采用邊界面塑性PM4Sand模型，結合相關文獻中某場地工程地質剖面圖，得到該場地的基巖出露地震波，計算該可液化緩坡場地在地震作用下自由面的側移值。主要結論如下：

1)PM4Sand模型能夠根據可液化土的動力觸探標準貫入值，獲得可液化土在PM4Sand中的相關特征值，模擬可液化土在地震中的強度衰減、孔隙水壓力生成和發展，以判定場地是否液化并引發滑移或側移。

2)在近場地地震波作用下，計算得到的場地側移值較文獻中記錄的側移值大，非線性分析結合PM4Sand模型可以對側移值進行預測。

3)在工程地質勘察不完善或者缺少可靠地震波的情況下，FLAC分析所得的側移或側移值須與其他預測方法的預測值進行對比分析，從而獲得更可靠的側移預測值。