地震作用下飽和土地基的彈塑性動力固結(jié)有限元分析

靳聰聰，顧曉偉

地震作用常導(dǎo)致地基液化，進而引起基礎(chǔ)的不均勻沉降和上部結(jié)構(gòu)的破壞，尤其是對于沿江河和近海岸填筑起來的陸地，在地震荷載下更容易引起人工回填砂或其下部砂土層超孔隙水壓力迅速聚集產(chǎn)生液化［1］。如1995年日本Ms7．2級阪神地震導(dǎo)致沿大阪灣北邊海岸線的大部分地區(qū)發(fā)生天然土層和人工填土的大面積液化現(xiàn)象；神戶海岸邊從東到西約20 km長的地區(qū)發(fā)生液化，海岸線向海移動2 m～ 3 m 和地表下沉近 1 m［2－3］。

國內(nèi)外諸多學(xué)者從現(xiàn)場調(diào)查、室內(nèi)試驗到數(shù)值方法分析對地震荷載作用下的地基動力分析進行了大量的研究［4－6］。其中數(shù)值方法可以考慮多種復(fù)雜加載條件，對試驗和理論方法難以解決的問題進行分析。動力分析方法從最初的總應(yīng)力法到現(xiàn)在考慮地震過程中土體孔隙水壓力擴散和消散的擬有效應(yīng)力方法和動力固結(jié)方法。其中動力固結(jié)方法在理論上比其它兩種方法更為嚴密，該方法能考慮水土耦合作用，即孔隙水壓力的產(chǎn)生、消散和擴散與土變形是完全耦合，因此可合理的描述動力作用下的土體應(yīng)力應(yīng)變的真實過程，并在工程中得到廣泛應(yīng)用［7－10］。

采用PZC模型對近海地基土層動力特性進行模擬，并用編制的可視化土體三軸試驗?zāi)M程序SM2D對PZC模型參數(shù)進行標定，采用國際上著名的飽和土動力固結(jié)分析程序 SWANDYNE II［11－13］的三維程序DYNE3WAC對該地基進行動力響應(yīng)分析。

1 動力固結(jié)分析方法

1．1 動力固結(jié)控制方程

地震作用下的飽和地基響應(yīng)問題，需充分考慮土體和液體各自性態(tài)和兩者間的耦合作用，飽和土體的動力固結(jié)理論的基本控制方程形式如下：整體平衡方程為：

式中：σij為總應(yīng)力；ui為土體位移；wi為孔隙流體相對于土體位移；bi為單位質(zhì)量上的力；ρf為流體密度；ρ為土體飽和密度。流體平衡方程為：

式中：p為孔隙水壓力；kf為動力滲透系數(shù)；n為孔隙率。流體連續(xù)方程可表達為

式中：εii為體應(yīng)變；α和Q為土體和流體的壓縮性相關(guān)的系數(shù)，有如下關(guān)系式：

式中：Ks為土顆粒體積模量；Kb為土骨架體積模量。根據(jù)孔隙率通過土體模量 Ks和流體模量 Kf得到，

式中：Kf為流體體積模量。土體本構(gòu)方程可表達為

式中：λ，μ為土體 Lame常數(shù)；σ′ij為有效應(yīng)力；δij為Kronecker張量運算符號。

對于包括地震工程在內(nèi)的大部分中低振動頻率的工程問題，可忽略流體相對加速度，根據(jù)式（1）～式（6），經(jīng)推導(dǎo)可得到僅含基本未知量 u和p的飽和土體動力固結(jié)控制方程（u－p形式）：

1．2 有限元離散

動力固結(jié)控制方程需進行離散后構(gòu)造有限元格式。首先對動力固結(jié)方程進行空間離散，然后再在時間域上采用適當?shù)姆椒ㄟM行求解。取土體位移和流體孔壓的形函數(shù)，然后對動力固結(jié)控制方程進行空間離散。令ui＝Nuˉui，p＝Npˉp，其中ˉui，ˉp分別為節(jié)點上的土骨架位移和孔隙水壓，Nu，Np分別為固相和液相的插值函數(shù)。結(jié)合格林公式對其進行分部積分，并引入邊界條件，最終可得到矩陣形式的動力固結(jié)方程：

式中：B為應(yīng)變矩陣；D為彈性矩陣；t為應(yīng)力邊界條件上的荷載；q為邊界處的法向流量。［M］為土體的質(zhì)量矩陣；［K］為土體的剛度矩陣；［Q］為耦合矩陣，范穎飽和土體中土體和流體的耦合作用；［H］為流體的滲透矩陣；［S］為流體的壓縮矩陣；Fu和Fp為固相的荷載和液相的荷載的右端項。

SWANDYNE－II采用無條件穩(wěn)定的廣義Newmark法對飽和土體的動力固結(jié)控制方程進行時域離散。對位移u采用GN22二階積分方案，t＝tn＋1時刻的加速度ˉu¨n＋1、速度ˉu·n＋1和位移ˉun＋1分別按下式計算：

式中：β1和β2是積分參數(shù)。

對孔壓p采用一階積分方案，t＝tn＋1時刻的孔壓變化率ˉ·p和孔壓ˉp分別按下式計算：

n＋1n＋1

式中：ˉβ1為積分參數(shù)。

如此逐步計算，可以得到整個時程中每個時步的位移和孔壓。

2 PZC彈塑性本構(gòu)模型

動力固結(jié)分析模擬能力主要取決于程序中所采用的動力本構(gòu)模型。采用PZC彈塑性本構(gòu)模型描述土體動力特性［14－15］。該模型直接確定塑性流動方向、加載方向和塑性模量，能夠較好描述土體主要靜動力特性，尤其是砂土非線性，剪脹剪縮及循環(huán)累計殘余變形，并能考慮卸載時塑性應(yīng)變，解決普通彈塑性模型在屈服面內(nèi)無塑性應(yīng)變、無法計算循環(huán)荷載引起的塑性應(yīng)變積累的問題［16］。

PZC模型共有靜動力參數(shù)13個，可以由排水或者不排水條件下的常規(guī)三軸單調(diào)、循環(huán)加卸載試驗得到，采用編制的可視化土體三軸試驗?zāi)M程序SM2D，對PZC模型參數(shù)進行標定。

PZC模型中彈性剪切模量如下所示：

式中：G0為模型參數(shù)；e為孔隙比；pa為大氣壓，其值為 101．325 kPa。

彈性體積模量可采用和式（13）類似的形式：

式中：K0是模型參數(shù)。

根據(jù)剪脹比dg和應(yīng)力比η的試驗結(jié)果發(fā)現(xiàn)剪脹方程近似表示為：

式中：φc為砂土殘余內(nèi)摩擦角；θ為羅德角。剪脹方程可用來計算塑性流動方向：

卸載時的塑性流動方向：

依據(jù)剪脹比定義塑性流動方向可較合理反映不同條件下砂土剪脹剪縮特性。

在PZC模型中，模型采用非相適應(yīng)流動法則，其加載方向與塑性流動方向的是分別定義的，加載方向nf具有和塑性流動方向 ng相似的形式：

式中：αf和 Mf為模型參數(shù)。

加載時的塑性模量HL表達式如下所示：

式中：β0與β1為模型參數(shù)，考慮剪切硬化的影響，ζ為塑性等效剪應(yīng)變累計值，表達式如下：

隨ζ不斷增加，Hs逐漸趨近為零。Hv、Hs兩者趨近于零時，土體接近于破壞狀態(tài)，因此可以反映不同密實度砂土在單向加載時的主要力學(xué)特性。

土體卸載完成后再加載時，土體為超固結(jié)土，需記錄土體的應(yīng)力歷史，故引入Hdm描述邊界面塑性映射規(guī)則，其表達式如下：

式中：γdm為模型參數(shù)；ζ表達式如下：

卸載時的塑性模量根據(jù)當前應(yīng)力是否到達臨界狀態(tài)區(qū)分為以下兩種情況：

式中：Hu0和γu為模型參數(shù)。

由以上各項定義得到PZC模型中的塑性應(yīng)變增量與應(yīng)力增量之間的關(guān)系：

進而可以推導(dǎo)出PZC模型的彈塑性矩陣Dep：

通過PZC彈塑性本構(gòu)模型公式分析可知，其彈性部分參數(shù)為 K0、G0、p0塑性勢面的參數(shù)為αg和Mg，屈服面的參數(shù)為αf和 Mf，塑性模量的參數(shù)為H0、β0、β1、Hu0、γu、γdm。

3 算例分析

3．1 三維有限元建模

在近海岸填筑起來的陸域，在地震荷載下容易導(dǎo)致人工回填砂的超孔隙水壓力快速累積并產(chǎn)生液化。以某近海岸填筑區(qū)域地基為例，對該場地一部分進行自由場三維動力分析。該處地層有約5 m厚的人工回填砂和10 m厚的全新世砂土，取長寬均為30 m的一塊場地進行液化自由場分析。根據(jù)該場地地層建立三維有限元網(wǎng)格，采用DYNE3WAC對該場地進行靜動力分析，并對不同深度的5個典型土質(zhì)點進行時程分析，該處陸域有限元模型如圖1所示。

采用DYNE3WAC進行有限元分析，土體單元為八節(jié)點六面體單元，用水土兩相飽和孔隙介質(zhì)模擬。

圖1 自由場地基有限元網(wǎng)格示意圖

該模型的邊界條件為：底邊作為固定邊界；遠置的側(cè)向邊界設(shè)定為固定；其余兩個側(cè)面的法向位移為零。

孔壓邊界條件為：表面結(jié)點的孔隙壓力為零；底邊和側(cè)邊均為不透水邊界。

3．2 參數(shù)標定

地基中兩種土層均采用PZC彈塑性本構(gòu)模型進行計算模擬。人工填土根據(jù)動三軸試驗結(jié)果進行參數(shù)調(diào)整，采用MATLAB／GUI編制的土體三軸試驗?zāi)M輔助程序SM2D對試驗進行模擬，并與試驗結(jié)果比較，進而標定本構(gòu)模型參數(shù)，人工填土的有效應(yīng)力路徑和應(yīng)力應(yīng)變關(guān)系模擬曲線如圖2所示。

圖2 人工填土PZC模型參數(shù)擬合用戶界面程序

同理，對全新世砂土使用SM2D程序模擬結(jié)果，并與試驗結(jié)果比較，確定本構(gòu)模型參數(shù)，砂土有效應(yīng)力路徑和應(yīng)力應(yīng)變關(guān)系模擬曲線如圖3所示。

圖3 砂土有效應(yīng)力路徑與應(yīng)力應(yīng)變關(guān)系

根據(jù)SM2D程序標定PZC模型的靜動力參數(shù)如 表2所示。

表2 模型PZC 彈塑性本構(gòu)模型參數(shù)

3．3 動力輸入

在計算中，對模型施加如圖4所示的雙向地震加速度荷載，主震峰值加速度為0．244 g和0．051 g，持時為40 s。為了研究模型孔隙水壓力在地震荷載作用下產(chǎn)生和消散的過程，計算時間共計100 s，計算時間步長為0．01 s。

圖4 地震水平及豎向加速度時程曲線

3．4 超靜孔壓時程分析

三維地基中的P1到P5的5個不同深度處孔隙水壓力在地震動作用下的擴散時程曲線和消散時程曲線如圖5所示。

圖5 地基不同深度超靜孔壓力形成的時程曲線

由圖5可知，在地震動往復(fù)作用下，地基5個典型土質(zhì)點超靜孔壓明顯增大，且隨著深度不斷增加對應(yīng)孔壓增加。地基底部中（P1）、地基中部（P2、P3、P4）、地基上部（P5）對應(yīng)的峰值孔壓分別為 89．8 kPa、72．3 kPa、60．5 kPa、47．7 kPa、27．5 kPa。這是因為在受到地震剪應(yīng)力作用時，由于砂土的剪縮剪脹性，導(dǎo)致砂土體積收縮，引起超孔隙水壓力快速上升，其中底部砂土區(qū)域產(chǎn)生高孔壓區(qū)，并引起孔隙水向上覆土體中擴散。

地震初期超靜孔壓有明顯增長，5個典型土質(zhì)點到達峰值孔壓時刻不同，其中P5和P4峰值對應(yīng)的時間分別為 13．1 s和 13．9 s，而 P1 在 15．2 s達到峰值孔壓。P4峰值孔壓相對于P5滯后，這種滯后現(xiàn)象表明液化區(qū)域由淺層向深層發(fā)展。

在地震作用期間，由于作用時間較短，飽和土中的孔隙水不能及時排出，因此土體中超孔壓逐漸增大。當土體的超孔壓增大到一定程度后，土體將發(fā)生液化。超孔壓比（EPWPR）是土體的超孔壓與其初始豎向有效應(yīng)力的比值，一般來說，當超孔壓比大于0．95時，認為土體發(fā)生液化（見圖6）。

圖6 地基不同深度超孔壓比的時程曲線

從圖6可知，典型土質(zhì)點P5深度的土體在13 s時超靜孔壓比大于0．95，該土質(zhì)點以上區(qū)域均產(chǎn)生液化。從P5超靜孔壓比時程曲線可知，該處的超孔壓比在3 s～13 s內(nèi)迅速增大，此后基本保持高孔壓比狀態(tài)。隨著深度增加，土體的超孔壓比減小，超孔壓比上升的速度逐漸變緩了，這是因為初始圍壓較大的原因。

地震結(jié)束后的超靜孔壓的消散時程分析如圖7所示。

圖7 地基不同深度超靜孔壓比和超靜孔壓等值線圖

從圖7可以得出地震荷載作用結(jié)束后5個典型土質(zhì)點的超孔隙水壓力逐漸消散，各個土質(zhì)點的消散速率均不相同。其中P1消散最快，P5消散最慢，這是因為地震中和地震后的超靜孔壓均是自下而上擴散，因此地震后超靜孔壓的消散也是底部土質(zhì)點先發(fā)生消散，最終到第100 s，各個土質(zhì)點的超靜孔壓均消散完畢，土體恢復(fù)到震前狀態(tài)。

3．5 位移時程曲線

三維地基中的P5的豎向沉降時程曲線如圖8所示。

圖8 P5豎向沉降時程曲線

從圖8可知，在雙向地震動作用下，典型土質(zhì)點P5豎向沉降呈波動變化，地基由于產(chǎn)生塑性變形而出現(xiàn)逐漸增大的沉降，地震結(jié)束后產(chǎn)生一定的殘余變形，地基的最大沉降為0．1 m。其中P5在第13 s發(fā)生液化之后，豎向沉降波動振幅減小，且沉降變化趨勢減小，到第20 s時基本不再發(fā)生沉降。這種現(xiàn)象說明當?shù)鼗饎涌讐荷仙?dǎo)致部分地基液化，該處土體剛度明顯降低，阻尼增大，地震剪切波向上傳播受到明顯的抑制，因此沉降波動減小。

4 結(jié) 論

基于PZC彈塑性模型和動力固結(jié)理論的DYNE3WAC彈塑性動力固結(jié)三維程序能夠合理分析近海填筑陸域在地震作用下的動力響應(yīng)。利用編制的土體三軸試驗?zāi)M輔助程序SM2D對該陸域土層材料進行模擬，標定土層材料的PZC模型參數(shù)。三維程序DYNE3WAC能全面考慮水土耦合和土體的動力特性，更加合理的模擬土體的變形和孔隙水壓力發(fā)展。計算結(jié)果表明，DYNE3WAC程序能夠模擬地震動荷載下土體的沉降變形、孔隙水壓力的波動上升和消散過程。

［1］ 黃 雨，于 淼，Bhattacharya S．2011年日本東北地區(qū)太平洋近海地震地基液化災(zāi)害綜述［J］．巖土工程學(xué)報，2013，35（5）：834-840．

［2］ 陳國興，金丹丹，常向東，等．最近20年地震中場地液化現(xiàn)象的回顧與土體液化可能性的評價準則［J］．巖土力學(xué)，2013，34（10）：2737-2795．

［3］ 陳國興，顧小鋒，常向東，等．1989—2011期間8次強地震中抗液化地基處理成功案例的回顧與啟示［J］．巖土力學(xué)，2015，36（4）：1102-1118．

［4］ 胡記磊，唐小微，裘江南．基于貝葉斯網(wǎng)絡(luò)的地震液化概率預(yù)測分析［J］．巖土力學(xué)，2016，37（6）：1745-1752．

［5］ 鄧海峰，劉振紋，祁 磊，等．波浪作用下飽和砂土孔壓發(fā)展規(guī)律試驗研究［J］．水利與建筑工程學(xué)報，2017，15（3）：45-48．

［6］ 邵 琪．飽和砂土地震液化的網(wǎng)格自適應(yīng)數(shù)值方法研究［D］．大連：大連理工大學(xué)，2014．

［7］ 賈會會．地震作用下加高擴容尾礦庫的動力穩(wěn)定性研究［J］．水利與建筑工程學(xué)報，2017，15（4）：157-161．

［8］ 徐令宇，王國新，蔡 飛，等．可液化場地地震反應(yīng)完全耦合動力分析及其驗證［J］．地震工程與工程振動，2014，1（6）：136-144．

［9］ 吳永康，王翔南，董威信，等．考慮流固耦合作用的高土石壩動力分析［J］．巖土工程學(xué)報，2015，37（11）：2007-2013．

［10］ 董威信，周漢民，李全明，等．基于比奧動力固結(jié)理論的尾礦壩地震響應(yīng)分析［J］．中國安全生產(chǎn)科學(xué)技術(shù)，2016，12（12）：28-32．

［11］ 王忠濤，劉 鵬．考慮主應(yīng)力軸旋轉(zhuǎn)效應(yīng)與狀態(tài)依賴性的砂土本構(gòu)模型［C］／／中國土木工程學(xué)會土力學(xué)及巖土工程分會．中國土木工程學(xué)會第十二屆全國土力學(xué)及巖土工程學(xué)術(shù)大會論文摘要集．上海：中國土木工程學(xué)會，2015．

［12］ Chen H，Zheng J，Li Q，et al．Response of Pore Pressure and Effective Stress in Seabed to Waves Around a Permeable Submerged Breakwater［C］／／International Conference of Offshore Mechanics and Arctic Engineering，Volume 1：Offshore Technology：Offshore Geotechnics．Newfoundland：ASME．2015：V001710A013．

［13］ Li P，Song E X．A high-order time-domain transmitting boundary for cylindrical wave propagation problems in unbounded saturated poroelastic media［J］．Soil Dynamics ＆Earthquake Engineering，2013，48（5）：48-62．

［14］ 李宏恩，李 錚，徐海峰，等．Pastor－Zienkiewicz狀態(tài)相關(guān)本構(gòu)模型及其參數(shù)確定方法研究［J］．巖土力學(xué)，2016，37（6）：1623-1632．

［15］ Moscariello M，Cuomo S，Manzanal D，et al．Modelling of wetting tests for a natural pyroclastic soil［C］／／E3S Web of Conferences Paris：EDP Sciences，2016：17007．

［16］ Ye J，Jeng D，Wang R，et al．Numerical simulation of the wave-induced dynamic response of poro-elastoplastic seabed foundations and a composite breakwater［J］．Applied Mathematical Modelling，2015，39（1）：322-347．