基于OpenSees的砂土本構模型對比研究

周林祿, 蘇 雷, 邱志堅, 凌賢長,3, 張安琪(. 青島理工大學土木工程學院, 山東 青島 66033; .廈門大學建筑與土木工程學院, 福建 廈門 36005;

3. 哈爾濱工業(yè)大學土木工程學院, 黑龍江 哈爾濱 150090)

0 引言

砂土液化是造成地震災害的重要原因之一,其機理是飽和砂土在循環(huán)荷載作用下有效應力逐漸減小,孔隙水壓力逐漸上升,砂土由固態(tài)轉變流動狀態(tài)。在1964年Niigata地震和Alaska地震[1-2]、1976年唐山地震[3]、1989年San Francisco地震[4]、1995年Kobe地震[5]和2018年Palu地震[6]中,由飽和砂土液化引起的工程結構破壞造成了巨大的人員傷亡和財產(chǎn)損失。因此,研究飽和砂土液化問題對地震易發(fā)區(qū)的工程結構抗震具有重要意義。針對砂土液化特性,研究者開展了各種類型的室內(nèi)和原位試驗。如Adalier等[7]通過離心模型試驗研究了超固結比(OCR)對飽和砂土液化的影響,發(fā)現(xiàn)隨著OCR的增加,引起超孔壓累積所需的加速度閾值將會增加。周燕國等[8]對含黏粒砂土進行了離心機振動臺試驗,發(fā)現(xiàn)含黏粒砂土與潔凈砂之間滲透系數(shù)的差異性導致振動過程中二者超孔壓的產(chǎn)生和消散速率顯著不同。陳國興等[9]開展了一系列循環(huán)三軸試驗,研究飽和砂礫土的兩種液化破壞機理:均等固結條件下的液化破壞和非均等固結條件下的軸向應變累積破壞。付海清等[10]基于現(xiàn)場液化試驗研究飽和砂土的動力特性,試驗結果表明相對密度對孔壓累積速率能起到明顯的控制作用。Chattaraj等[11]基于共振柱試驗和循環(huán)三軸試驗研究了砂土的動剪切模量和阻尼等特性,并建立了砂土最大動剪切模量和循環(huán)強度的關系。Dash等[12]研究了循環(huán)荷載頻率對飽和砂土抗液化性的影響,發(fā)現(xiàn)隨著循環(huán)加載頻率的增加,超孔壓累積速率也增加,即抗液化性隨著加載頻率的增加而降低。Moayed等[13]通過循環(huán)三軸試驗,發(fā)現(xiàn)土工布在試樣中的排列對砂土動力特性起主導作用,當土工布靠近試樣頂部時砂土抗液化性能增加。Amini等[14]、Sonmezer[15]分別基于循環(huán)三軸試驗和循環(huán)剪切試驗,采用能量法分析了砂土的液化特性。Wang等[16]研究了不同波形的動荷載對砂土動力特性的影響。Yazarloo等[17]通過納米技術對砂土進行不同程度的加固,分析其抗液化性。

綜上所述,研究者在砂土液化領域做了很多室內(nèi)和現(xiàn)場試驗,但這些試驗存在成本高和周期長等局限性,尤其是大型試驗(如振動臺或離心機試驗)和參數(shù)分析。相比之下,在典型動力試驗結果的基礎上利用數(shù)值模擬標定現(xiàn)有的本構模型,進而進行參數(shù)分析,能夠有效避免上述問題。如凌賢長等[18]針對非自由可液化場地大型振動臺試驗,采用有效應力原理和應力循環(huán)孔壓增量計算方法建立數(shù)值計算模型,計算結果顯示孔壓比的增長與試驗結果基本一致。周健等[19]采用二維顆粒流程序PFC2D分析了砂土振動液化時的宏觀力學響應和細觀組構演化。唐亮等[20]通過有效應力原理建立土體液化分析模型,并基于該模型建立自由液化場地振動臺試驗的有限元模型,通過試驗值和計算值的對比驗證了土體液化模型的可靠性。Zhang等[21]提出了預測飽和砂土液化后變形的理論框架,能夠保證在有效應力消失時數(shù)值模擬的穩(wěn)定性。Shahir等[22]通過引入滲透系數(shù)和超孔壓比的關系對完全耦合的砂土動力學模型進行修正,結果表明修正后的模型能更準確地模擬砂土變形和孔壓的累積過程。Wang等[23-24]針對飽和砂土進行了離散元數(shù)值模擬,發(fā)現(xiàn)砂土細觀結構特征的變化會影響其抗液化性和液化后剪應變的發(fā)展,并建立了砂土細觀結構變化和抗液化性的關系。魏星等[25]基于離散元的數(shù)值模擬方法,發(fā)現(xiàn)由循環(huán)剪切引起的孔隙均勻化和土顆粒間孔隙增大是液化后循環(huán)剪應變逐漸增大的細觀原因。

目前已開發(fā)的用于模擬砂土液化的本構模型有很多,而對不同本構模型的合理選用是一個關鍵問題。王剛等[26]詳細分析了不同本構模型模擬砂土液化變形的機制,并進行了不同模型響應對比。Ramirez等[27]利用兩種數(shù)值計算方法,系統(tǒng)地比較了兩個砂土本構模型(PDMY02模型和SANISAND模型)在地震荷載作用下土層的沉降、加速度和超孔壓響應特性;通過對比離心機試驗結果發(fā)現(xiàn):兩個本構模型均能夠很好地預測土層的地震響應,但在模擬土體沉降方面仍需要改進。Sun等[28]基于修正的臨界狀態(tài)線發(fā)展了一種砂土邊界面塑性本構模型,并與Li等[29]提出的模型進行了比較,結果表明邊界面塑性本構模型能夠再現(xiàn)砂土典型的應力-應變響應。

因此,為了盡可能準確地模擬地震荷載下飽和砂土的響應特性,基于開源地震工程模擬系統(tǒng)OpenSees,對其材料庫中的4種典型砂土本構模型:PDMY(Pressure Depend Multi Yield)模型、Stress Density模型、Manzari Dafalias模型和CycLiqCPSP模型(砂土液化大變形統(tǒng)一本構模型)進行模擬測試,分析并對比每個模型的動力響應結果,以期為相關數(shù)值模擬提供參考和借鑒。

1 飽和砂土本構模型簡介

1.1 PDMY模型

PDMY模型是由Yang等[30]提出的一種彈塑性飽和砂土液化本構模型,主要用于模擬循環(huán)荷載條件下對圍壓敏感的土體材料的非線性響應特性。模型屈服函數(shù)遵循經(jīng)典的塑性約定。土體材料的彈性響應認為是線性和各向同性的,而塑性響應認為是非線性和各向異性的。在應力空間坐標系中,許多具有共同頂點和不同尺寸的相似屈服面形成了硬化區(qū),每個屈服面都與一個塑性剪切模量值相關聯(lián),最外側屈服面為破壞面。該模型采用偏分量相關聯(lián)和體積分量非關聯(lián)的流動規(guī)則,土體的剪脹和剪縮特性由非相關聯(lián)流動法則控制。

1.2 Stress Density模型

Stress Density模型是由Cubrinovsk等[31]基于初始狀態(tài)概念的框架構建的一種砂土應力-應變-剪脹模型。該模型基于孔隙比和法向應力狀態(tài)相關聯(lián)的思想建立,其應力-應變關系基于初始狀態(tài)概念建立,通過定義狀態(tài)指數(shù)變量來量化相對初始狀態(tài)。模型采用連續(xù)屈服假設定義了塑性應變增量方向與應力增量方向的依賴關系。

1.3 Manzari Dafalias模型

Manzari Dafalias模型[32]是由應力比控制、同時考慮臨界狀態(tài)的飽和砂土塑性模型。該模型基于一般土體本構模型的特征,對模型結構的簡單性與模型功能的復雜性進行了適當平衡。模型假定只有應力比的變化才能引起砂土塑形剪切變形和體積的變化,而在恒定的應力比下應力的增加僅由彈性應變引起,對于塑性應變則通過引入其他塑性加載機制來實現(xiàn)。模型采用剪脹張量模擬循環(huán)荷載下的剪脹和剪縮響應特性。

1.4 CycLiqCPSP模型

CycLiqCPSP模型是由Wang等[33]基于切割平面算法構建的一種飽和砂土液化本構模型,特別考慮了循環(huán)特性和液化后剪切大變形的累積。該模型基于物理方法計算零有效應力下剪切應變的累積,采用統(tǒng)一方法描述砂土液化前后的特性。模型基于剪脹引起的兩個體積應變分量分析砂土在單調(diào)和循環(huán)荷載下的剪脹特性,同時假定由平均有效應力引起的體積應變?yōu)閺椥?由剪切變形引起的體積應變?yōu)樗苄?。模型還引入了臨界狀態(tài)土力學概念,可以實現(xiàn)不同狀態(tài)下砂土的統(tǒng)一建模。

2 數(shù)值模型途徑

2.1 數(shù)值模型建立

為研究以上4種飽和砂土本構模型在相同荷載作用下響應的差異性,建立一個單元的二維飽和砂土流-固耦合測試模型(圖1)。模型尺寸為1 m×1 m,采用OpenSees中的9_4_QuadUP單元。該單元基于Biot多孔介質(zhì)理論,采用u-p公式的有限元形式(u表示土顆粒位移,p表示孔隙水壓力)[34],考慮了土骨架和孔隙水之間的相互影響,用于模擬固-液兩相完全耦合材料的動力響應。如圖1所示,土單元的4個角節(jié)點各具有3個自由度,其中前兩個自由度代表土顆粒位移,第3個自由度代表孔壓自由度。另外5個節(jié)點具有2個自由度,僅代表土顆粒位移。

圖1 二維飽和砂土流-固耦合模型Fig.1 Two-dimensional fluid-solid coupling model of saturated sand

單元建立過程中,首先對節(jié)點1、2、5在x和y方向施加固定約束,保證單元初始位移為0;對節(jié)點3、4施加孔壓約束,定義單元上邊界為水位線。其次分別對節(jié)點3、4、7和節(jié)點6、8、9施加相同位移約束,以保證模型在動荷載作用下維持純剪切狀態(tài)。最后在單元底部施加如圖2所示的正弦加速度,加速度頻率為1 Hz,幅值為0.2g,整個施加過程持續(xù)10 s。數(shù)值計算過程分成兩步:第一步計算土體在自重狀態(tài)下的靜力響應,分別實施土單元的彈性和塑性計算;第二步計算施加加速度后單元的動力響應。

圖2 基底輸入加速度Fig.2 Base input acceleration

動力分析中,為單元和節(jié)點分配阻尼。在OpenSees中使用瑞利阻尼時,單元或節(jié)點的阻尼矩陣D被指定為剛度和質(zhì)量比例阻尼矩陣的組合:

D=αM+βK

(1)

式中:M和K分別為質(zhì)量矩陣和剛度矩陣;α和β分別為質(zhì)量矩陣系數(shù)和剛度矩陣系數(shù),在對不同本構模型進行計算時,α和β統(tǒng)一取值,分別設置為0和0.002。

2.2 模型輸入?yún)?shù)確定

為了更好地對比不同砂土本構模型在相同荷載作用下響應的差異性,需要盡可能地保證定義的不同模型參數(shù)代表同一類砂土。因此,各個模型基準參數(shù)的選取參考了以往研究者對相似密度下中砂的模型標定結果[30-33,35]。

考慮到剪切模量是飽和砂土最重要的模型參數(shù),需要確保在循環(huán)荷載作用下不同本構模型具有類似的剪切模量。PDMY模型的剪切模量利用剪應力-剪應變骨干曲線描述,基于指定的參考圍壓與該參考圍壓下得到的小應變剪切模量參考值,建立本構模型中剪切模量隨平均有效應力的變化關系:

(2)

式中:Gr為參考剪切模量;p′r為參考有效圍壓;d為壓力系數(shù)。

Stress Density模型的剪切模量基于Iwasaki等[36]的研究成果,在建立剪切模量與有效應力關系時考慮了粒度分布對剪切模量的影響:

(3)

式中:A為彈性剪切模量參數(shù);pa為大氣壓(參考圍壓);ei為孔隙比;n為壓力系數(shù)。

Manzari Dafalias模型和CycLiqCPSP模型則根據(jù)Richart等[37]定義的經(jīng)驗公式建立剪切模量表達式,同時還考慮了初始孔隙比對剪切模量的影響:

(4)

式中:G0為剪切模量參數(shù)。

為消除數(shù)值模擬過程中剪切模量對響應結果的影響,對式(2)～(4)中的相關參數(shù)進行修正。圖3為修正后各本構模型剪切模量與平均有效應力的變化關系。由圖3可知,在這4個模型中,剪切模量與平均有效應力的變化關系曲線基本重合。表1～4分別列出了4個本構模型的輸入?yún)?shù),表5為數(shù)值模型中所采用單元的相關計算參數(shù)。

圖3 剪切模量與平均有效應力的變化關系Fig.3 Relationship between shear modulus and mean effective stress

參數(shù)數(shù)值密度ρ/(kg·m-3)1 900參考剪切模量Gr/kPa8.1×104參考體積模量Br/kPa2.0×105摩擦角φ/(°)33峰值剪應變γ0.1參考有效圍壓p′r/kPa101.3壓力系數(shù)d0.5相位轉換角φT/(°)27剪縮參數(shù)c0.07剪脹參數(shù)d1,d20.4,2液化參數(shù)l1, l2, l310,0.01,1

表2 Stress Density模型參數(shù)

表3 Manzari Dafalias模型參數(shù)

表4 CycLiqCPSP模型參數(shù)

表5 單元相關參數(shù)

3 計算結果分析

基于上述數(shù)值建模途徑,采用OpenSees有限元模擬程序進行數(shù)值計算,得到不同飽和砂土本構模型的動力響應結果。下面將從加速度、超孔隙水壓力、位移、剪應力-剪應變和平均有效應力路徑響應方面進行分析,研究各本構模型的動力響應規(guī)律,對比不同模型響應的差異性。

3.1 加速度響應

圖4給出了4種模型土單元節(jié)點3的水平加速度隨時間變化關系。由圖可知,PDMY和Cyc-LiqCPSP模型的節(jié)點加速度響應峰值和波形基本保持一致,與輸入的加速度相比,波形改變的幅度較小。另外注意到,Stress Density模型和Manzari Dafalias模型從第2個循環(huán)開始,加速度表現(xiàn)出較明顯的放大效應。對于不同模型,該放大效應的大小表現(xiàn)出一定的差別。

圖4 不同砂土本構模型節(jié)點3加速度時程Fig.4 Acceleration time histories of Node 3 for different sand constitutive models

3.2 超孔壓響應

圖5給出了循環(huán)動力荷載作用下不同本構模型飽和砂土單元在節(jié)點1處的超孔壓變化時程曲線。PDMY模型在起始階段,隨著加速度的施加超孔壓逐漸增加,但增加的速率逐漸降低;在第6 s,超孔壓增加到最大值,即達到初始有效應力,飽和砂土單元發(fā)生液化。Stress Density模型在加速度施加的第1 s出現(xiàn)了短暫的負超孔壓,這是因為該模型材料在循環(huán)荷載作用下產(chǎn)生了較明顯的剪脹和剪縮效果,即當加載時,砂土單元表現(xiàn)為加載膨脹,其內(nèi)部孔隙在一瞬間發(fā)生擴張,從而導致孔隙水壓力減小;當卸載時,砂土單元表現(xiàn)為卸載收縮,此時超孔壓開始上升[38]。

另外由圖5可知,超孔壓波動頻率比輸入加速度的頻率高,這與砂土在動荷載作用下的剪縮、剪脹交替變化有關。當剪應變增大時,孔壓出現(xiàn)回降(加載剪脹),而當剪應變回退時,孔壓又迅速上升(卸載剪縮)。由于每個完整的加速度周期有兩個加載段(正負兩個方向),所以超孔壓波動的頻率為輸入加速度頻率的2倍。不斷交替的剪脹和剪縮作用造成了Stress Density模型土單元超孔壓的急劇升高和降低,僅加載兩個周期的循環(huán)時間,超孔壓即在某些瞬間達到初始有效應力,砂土單元開始連續(xù)出現(xiàn)“瞬態(tài)液化”現(xiàn)象,且超孔壓持續(xù)在比較大的范圍內(nèi)波動。Manzari Dafalias模型同樣在加載兩個周期的循環(huán)時間出現(xiàn)了連續(xù)的“瞬態(tài)液化”現(xiàn)象,但相比Stress Density模型,其剪脹特性稍弱。該模型響應的初始階段未出現(xiàn)負超孔壓現(xiàn)象,在液化后超孔壓一直保持中等幅度的波動,持續(xù)到加載8 s左右,之后超孔壓波動幅度變小。結合圖4不難發(fā)現(xiàn),剪脹特性越強,其加速度的響應峰值也越高。Cyc-LiqCPSP模型在循環(huán)正弦荷載作用下,超孔壓先是呈穩(wěn)定的小幅值增長,隨后在第2 s增長到最大值,發(fā)生完全液化,緊接著發(fā)生剪脹作用,之后超孔壓保持最大值且呈直線發(fā)展。

圖5 節(jié)點1超孔壓響應Fig.5 Excess pore pressure response of Node 1

3.3 位移響應

在加速度施加過程中,記錄了土單元每個節(jié)點的位移響應。圖6對比了不同本構模型在節(jié)點3的水平位移。由圖可知,PDMY模型在0～6 s時間段內(nèi),位移響應幅值呈逐漸增加的趨勢;在6～8 s,位移響應幅值開始保持穩(wěn)定狀態(tài);在8 s之后,由于輸入加速度幅值下降,位移呈減小趨勢。Stress Density模型的位移時程響應顯示,該模型在第1個加速度循環(huán)即產(chǎn)生較大的位移變化,之后位移向反方向發(fā)展,并在此基礎上進行波動。這說明該土單元發(fā)生了永久變形,即該模型在循環(huán)荷載下更容易呈現(xiàn)永久變形的特性。結合圖5可知,在永久變形發(fā)生時,超孔壓剛好達到初始有效應力(即液化),因此這種永久變形是由土體液化引起的。Manzari Dafalias模型在第2個加速度循環(huán)達到最大位移,但位移幅值較小,之后幅值近似維持著恒定值,直到循環(huán)荷載結束。CycLiqCPSP模型的位移時程與PDMY模型相似,位移幅值在起始時段逐漸增加至穩(wěn)定狀態(tài),但CycLiqCPSP模型更早到達最大位移響應幅值。 這4種模型中,Stress Density模型由液化產(chǎn)生的累計變形最為顯著,其他3種模型最終產(chǎn)生的累計變形很小,這與動荷載強度、加載時間以及單元尺寸等因素有關。

圖6 節(jié)點3位移響應Fig.6 Displacement response of Node 3

3.4 剪應力-剪應變響應

圖7給出不同本構模型在循環(huán)荷載作用下的剪應力-剪應變響應。從圖7可知,PDMY模型隨加速度循環(huán)次數(shù)的增加,土單元剪應力逐漸減小,剪應變逐漸增大。當節(jié)點位移和超孔壓增長到最大值時,單元剪應力相應地減小到最小值,剪應變增長到最大值。值得注意的是,該模型的剪應力-剪應變響應呈階梯狀變化,即剪應力和剪應變分別以階梯的形式逐漸減小和增大。Stress Density模型在加載的第2個周期內(nèi)產(chǎn)生了較大的應變變化,此后滯回曲線一直停留在變化后的位置。這與位移響應結果相符合,說明該模型在動荷載作用下出現(xiàn)了累積永久變形。Manzari Dafalias模型在循環(huán)荷載作用下,其剪應力-剪應變響應始終穩(wěn)定在一個大致范圍內(nèi)。顯而易見,該模型的剪應變響應值較其他模型要小很多,這也是造成其位移響應值較小的直接原因。CycLiqCPSP模型隨著加速度的施加,土單元剪應力減小的同時剪應變也在增大;與此同時,剪切剛度隨著加速度循環(huán)次數(shù)的增加逐漸降低。當超孔壓增長到最大值時土單元發(fā)生液化,此時剪切剛度降為0。對比該模型與PDMY模型不難發(fā)現(xiàn),液化發(fā)生后兩個模型的剪應力-剪應變滯回曲線幾乎重合。

圖7 土單元剪應力-剪應變響應Fig.7 Shear stress-strain response of soil element

3.5 平均有效應力路徑響應

圖8給出了砂土單元平均有效應力路徑響應。由圖8可知,PDMY模型隨著動力循環(huán)荷載的施加,超孔壓不斷累積,從而導致平均有效應力不斷降低,直至降為0。砂土單元的強度不斷減小,位移幅值逐漸增大,最終達到完全液化狀態(tài)時,土體完全喪失強度,此時節(jié)點位移增加到最大值。Stress Density模型在初始階段,由于負超孔壓的出現(xiàn),平均有效應力有所增長;在之后的動荷載循環(huán)中,平均有效應力先是迅速衰減到最小值,隨后立刻恢復,連續(xù)在其強度包線內(nèi)進行大范圍的波動。這與超孔壓和單元剪應力的響應結果一致,都是由其剪脹特性所引起。Manzari Dafalias模型在第2個加載循環(huán)內(nèi),平均有效應力降到最小值,隨后在其強度包線內(nèi)波動,但波動的范圍明顯要小于Stress Density模型。CycLiqCPSP模型隨著循環(huán)荷載的施加,由于超孔壓的逐漸累積,其平均有效應力逐漸減小,在材料開始液化時,有效應力降低到最小值。需要注意的是,Manzari Dafalias模型和CycLiqCPSP模型的平均有效應力并沒有精確衰減到0點,這是因為材料在達到液化后尚存在殘余強度。

圖8 平均有效應力路徑Fig.8 Mean effective stress path

4 結論

針對飽和砂土液化問題,通過開源地震工程數(shù)值計算平臺OpenSees,對材料庫中PDMY模型、Stress Density模型、Manzari Dafalias模型和CycLiqCPSP模型進行單一單元的數(shù)值計算,分析相同循環(huán)動力荷載下,不同砂土本構模型動力響應的差異性,得到以下結論:

(1) 在模擬過程中,節(jié)點加速度表現(xiàn)出一定的放大效應,且這種放大效應與材料的剪脹特性存在緊密關系,剪脹特性越強,加速度放大效應越明顯。

(2) 在4種砂土本構模型中,Stress Density模型對循環(huán)荷載作用下土體永久變形的模擬效果最好,而PDMY模型和CycLiqCPSP模型在模擬土體可恢復變形時表現(xiàn)出一定優(yōu)勢。

(3) PDMY模型和CycLiqCPSP模型的超孔壓存在穩(wěn)定累積的過程,并且在增長到初始有效應力后保持不變。Stress Density模型和Manzari Dafalias模型的超孔壓表現(xiàn)出明顯的剪脹效應。

(4) PDMY模型和CycLiqCPSP模型在動荷載作用下剪切剛度逐漸降低,完全液化后剪切剛度降為0,且不再繼續(xù)變化。Stress Density模型和Manzari Dafalias模型在動荷載作用下剪切剛度時刻變化,沒有出現(xiàn)持續(xù)液化狀態(tài)。

(5) 在動荷載作用下,PDMY模型和Cyc-LiqCPSP模型的平均有效應力逐漸減小,達到最小值后保持不變。Stress Density模型和Manzari Dafalias模型在動荷載作用的某些瞬間平均有效應力降為最小值,而在另外大多數(shù)時刻始終處在波動狀態(tài)。

(6) 綜合考慮各個模型的響應結果,認為PDMY模型可以很好地模擬動荷載作用下飽和松砂超孔壓的累積效應和剪縮效應,而Stress Density模型和Manzari Dafalias模型則更適用于模擬中密和密砂的循環(huán)動力響應。