基于粘彈塑性模型的二維土體非線性地震反應分析

楊笑梅，陳鑫，游昊冉

（廣東工業大學土木與交通工程學院，廣東 廣州 510006）

引言

包含沉積土層的復雜場地地震反應及地表放大特征分析對防災減災及結構的抗震設計至關重要［1］。由于現有強震觀測臺站還不足以覆蓋所有城鎮區域，且強震觀測數據還遠不夠預測工程中需要的地表反應特征［2］，因此，地震波在復雜沉積土層傳播特征的研究仍主要依賴波傳播理論及數值模擬的手段［3］。目前通用有限元軟件在復雜場地數值模型建立上存在很大優勢［4］，但在覆蓋土層的非線性與吸能特征［5］及人工邊界設置上還很難達到要求。特別是對土體非線性特征的模擬，通用軟件內置的常規非線性本構模型如摩爾庫倫模型（Mohr-Coulomb）、D－P模型（Drucker-Prager）、劍橋模型等［6］并不適用。主要是因為這類模型均采用與動力分析不同的硬化規律，不能體現動力作用下應力路徑特征［7－8］。另外，工程中使用最廣泛的土動力本構模型是Idriss［9］提出的等效線性化模型。該模型雖然可以在一定程度上模擬土體的非線性及吸能特征，但其在軟土及強非線性反應中的模擬合理性一直受到質疑。因此尋求簡單有效的非線性本構模型，借助通用軟件使其適用于模擬土體高維非線性地震反應一直是復雜場地地震反應分析中重要的研究內容。

目前粘彈塑性本構仍是公認的能反映土體真實非線性特征的理論模型，其更適用于強震下的地震反應分析。Iwan［10］模型作為經典的土體動力非線性分析模型的代表，一直廣泛應用于土動力反應分析問題。很多學者在Iwan串聯模型基礎上進行改進，使其更適用于相關的研究問題。如為了合理考慮軟黏土的應變軟化特征，Rao等［11］借鑒Idriss等［12］軟化模型，進一步修正了Iwan串聯動力本構模型。在眾多描述土體非線性的粘彈塑性動力本構模型中，通常存在理論復雜，推導過程參數眾多等問題，且通過土動力學試驗確定的土體參數離散性較大，計算結果的可靠性常存在爭議。因此，很多模型還僅限于理論研究層面，難于應用到實際工程分析中。Joyner和Chen等［13－14］于1972年將Iwan串聯模型應用于高維非線性地震反應研究中，由于采用了合理硬化規律，不僅在模擬應變硬化土體時得到較好的結果，而且在應變較小時模擬應變軟化類土體也能達到較高的精度。特別是該模型相對于眾多的同類土體動力非線性模型來說，具有同等效線性化相類似的參數選取功能，較其他類似彈塑性本構模型具有參數簡單，易選取，模型穩定性較好等特征。

為了能借助通用軟件的強大數值模擬功能，并使其適合模擬包含土層的高維地震反應，文中從二維波動問題出發，基于ABAQUS的顯式計算過程：（1）推導了與其顯式計算相匹配的粘性邊界，解決了ABAQUS軟件無法直接模擬二維無限域波動問題的困難；（2）在Joyner［13］基于Iwan串聯模型發展的土體動力本構模型基礎上，考慮了土體的吸能特征，引入瑞利阻尼模型進行改進，使其適用于模擬土體的粘彈塑性分析；（3）借助通用軟件ABAQUS的良好開放性［15］，利用其中二次開發功能子程序VUMAT，將本構模型植入其中，克服ABAQUS自帶本構模型難以較好地模擬土體地震反應的劣勢。為了驗證文中方法，選取2個二維典型算例，進行地震反應分析并與參考解對比，吻合度較好。進一步選取文獻［16］的典型盆地為研究對象，通過與文獻［16］計算結果對比，探討了文中本構模型的參數選取原則，并闡述參數選取的靈活性及簡便性，為其能進一步應用于復雜工程場地地震反應分析提供指導。

1 2D P－SV波傳播模擬在ABAQUS中的實現

2D場地地震反應問題理論上屬于無限域二維P－SV波傳播問題，通用的有限元軟件無法直接模擬。需要解決人工邊界設置、土體的非線性本構模型及土體的阻尼模型建立等相關問題。為了更好地闡述文中構建的計算過程，本節將簡要介紹其相關理論及其在ABAQUS中的實現。

1.1 粘性人工邊界條件及地震輸入在ABAQUS中的實現

為了實現地震波在無限域傳播的合理模擬，需要在建立的有限域分析模型邊界設置特殊的邊界條件，以保證地震波傳播過程中不會在邊界處產生偽波。本計算過程采用由Lysmer等［17］提出的粘性人工邊界模型，通過在ABAQUS建立的有限元模型邊界處合理設置阻尼器實現。具體的阻尼器參數計算公式如式（1）：

式中：下標N代表法線方向；下標T代表切線方向；ρ為密度；vp為縱波波速；vs為剪切波波速；A為單元面積。

地震輸入以等效結點力F實現［18－19］，計算公式如式（2）：

式中：下標x或y表示施加于結點處力的方向；上標x或y表示外法線方向；Δt1及Δt2分別表示邊界處入射波和反射波的延遲時間。

1.2 多屈服面本構模型理論簡介

文中建立的本構模型是基于Joyner［13］在Iwan［10］模型基礎上改進的二維多屈面模型，由于該模型對應力路徑能夠合理記憶，因此在任意復雜的往復荷載作用下模擬其加卸載曲線均能達到較好的結果。為理解文中計算過程，先簡要介紹文獻［13］中的多屈面本構模型。

1.2.1 相關應力應變量的定義

文獻［13］所用的彈塑性模型將采用平均應力σM，平均應變eM，偏應力σij和偏應變eij的概念，定義分別如下：

式中：Sij和Eij分別表示總應力及總應變；δij為Kronecher符號。由式（6）及式（7）可知，總應力或總應變量是通過計算出平均應力或應變及偏應力或應變后求得。模型假設塑性應變全都是由偏應力產生，平均應力及平均應變關系為線彈性。

1.2.2 多屈面模型的假設

文中彈塑性模型假設主要采用文獻［13］的理論確定偏應力應變增量關系，該模型的屈服準則、加卸載準則、移動硬化定律以及流動法則描述如下：

（1）屈服準則

土體在偏應力空間的應力狀態下，假設有一系列的屈服面，滿足如下屈服函數：

式中：F是屈服函數；αij表示移動硬化參數；k表示與屈服強度相關的常數；下標n代表第n個屈服面。

（2）移動硬化準則

為了體現動力作用，Joyner采用式（10）作為移動硬化定律，從而保證應力點和屈服面的一致性：

式中：上標含“p”的物理量表示前一時刻的計算值，上標含“p+1”的物理量表示當前時間步的計算值。

（3）流動法則

文中采用如下的流動法則描述塑性變形增量的大小和方向：

這里新出現的變量Ln與加卸載準則相關，具體參數如下節，而hn為非負的塑性標量因子，由式（12）計算：

式中，cn為第n個屈服面的屈服面系數。

（4）加卸載準則

土體受到往復動力荷載作用時，假設的加卸載準則如下：

1.3 土體阻尼模型的選取

土體除了具有非線性特征還存在較強的吸能特征。選取在計算過程中保持定阻尼值的瑞利阻尼［20］考慮土體的耗能特征。參考文獻［21］給定的計算公式確定阻尼系數，即：

式中：α、β為瑞利阻尼系數；ξ為土體阻尼比系數；f1和f3分別為線性場地1，3階頻率。

1.4 彈塑性模型VUMAT子程序的實現流程

利用了ABAQUS二次開發編寫了對應文中彈塑性本構計算的VUMAT子程序，具體計算過程描述如下：（1）讀入由ABAQUS數據文件給定的初始總應力、總應變、初始狀態變量值，并根據式（6）及式（7）計算平均應力及偏應力，平均應變及偏應變初值；（2）調用VUMAT子程序，由ABAQUS中讀取的位移值，計算平均應變增量及偏應變增量；（3）由平均應變增量計算平均應力增量；（4）按彈性本構模型計算試探應力，并根據計算結果按屈服準則即式（8）及式（9）判別是否屈服，若判別結果為彈性，則子程序自動更新應力及應變值進行下一個增量步的計算；（5）若判別為塑性，按塑性理論要求計算，遵循文中的移動硬化準則、流動法則及加卸載準則由偏應變增量計算偏應力增量；（6）由計算所得的平均應力增量及偏應力增量，更新偏應力及總應力，并由平均應變增量及偏應變增量更新偏應變及總應變值；（7）重復（2）～（4）步直至計算停止。

2 基于ABAQUS數值模擬過程P－SV波傳播問題的驗證

為了驗證文中計算過程的正確性，本節通過2個算例分別驗證二維波傳播問題邊界設置的合理性及多屈面本構模型的可靠性。

2.1 半圓峽谷地形地表反應

為了驗證文中方法與計算的正確性，選取Kawase在文獻［22］中基于離散波速邊界元法分析的無限域半圓峽谷地表反應問題進行分析，為了便于結果的比較，類似文獻［22］，計算時定義一個與頻率相關的參數：

式中：λS為彈性介質中剪切波波長；νS表示彈性介質剪切波波速；a表示半圓半徑。文中η取2，介質計算參數如下：νs=200 m/s，密度ρ=1 800 kg/m3，泊松比υ=1/3；半圓峽谷模型計算參數取為：L×h=800 m×200 m，峽谷半圓半徑a=50 m。基于有限元軟件ABAQUS建立二維SV垂直入射線彈性有限元模型，如圖1（a）所示，邊界處網格尺寸為1 m×1 m。在邊界處設置阻尼器模擬外部無限域，如圖1（b）所示。根據式（1）計算邊界處相應阻尼器參數，并將地震輸入等效為結點力添加在邊界，邊界力值利用式（2）～式（5）計算，進行時域顯式分析，為了對比，地表水平與垂直方向的各點位移響應結果轉化為頻域結果，如圖2所示，圖2同時給出了Kawase計算結果。

圖1 半圓谷地有限元模型及局部粘性邊界Fig.1 The finite element（FE）model for semi-circular canyon and its local viscous boundary

從圖2結果對比可見。本文解與文獻［22］計算結果比較吻合，邊界反射的誤差波未對結果產生影響，也說明地震波可以由邊界處外傳出計算域，x表示該點到谷地中央的距離，x/a為無量綱物理量。因此文中基于ABAQUS二維顯式有限元結合設置人工邊界及地震輸入的計算過程可以合理模擬二維無限域線性地震反應。

圖2 半圓谷地地表位移幅值（η=2.0）Fig.2 The displacement amplitudes on surface of semi-circular valley for η=2.0

2.2 二維土體非線性算例

為了校驗本構模型，選取Willam B.Joyner文獻［13］的土壩模型算例進行模擬，計算模型如圖3所示，依據原文網格尺寸取為9 m×9 m，土壩土層主要由黏土與顆粒狀物質組成，土體動力參數請見表1，因文中本構模型參數主要包括抗剪強度、剪切波速、動剪切模量、屈服面數及阻尼參數，而文獻［13］未考慮阻尼故此例忽略該參數。地震動輸入為N21E方向的TAFT波，為了模擬強震輸入下土體的非線性特征，按原文的方法，將輸入地震波整體乘以4倍調整速度峰值為0.7 m/s，輸入地震動時程見圖4。參考原算例，在土壩的表面設置了5個輸出點（如圖3中1－5所示），模擬結果為速度時程，結果如圖5所示。為了比較，文獻［13］中對應5點的速度時程結果也繪于圖5。圖5顯示文中結果與文獻［13］結果吻合較好，驗證了文中本構模型的正確性。

圖3 土壩簡化FEM模型［13］Fig.3 Simplified FEM model for earth dam［13］

圖5 地表指定點（圖3所示）的水平速度時程Fig.5 The horizontal velocity time histories at the points on surface prescribed in Fig.3

表1 土體參數［13］Table 1 Soil parameters［13］

圖4 基底輸入地震動時程［13］Fig.4 Time history of base-input ground motion［13］

3 非線性反應分析參數設置問題的討論

彈塑性本構模型模擬土體非線性反應時，除上節給出的參數外，還需要設置很多與非線性分析相關的參數，參數選取時通常會遇到困難，導致很多本構模型很難運用于實際問題分析。相對多數彈塑性模型，文中本構模型參數相對簡單穩定，更適用復雜地質條件特別是含有復雜土介質的地震反應分析。本模型依賴的非線性分析參數主要包含抗剪強度、剪切波速、動剪切模量、屈服面參數以及粘性阻尼參數等，可以很容易通過實驗直接獲取或通過等效線性化模型的模量比曲線參數轉化，不需要增加額外的工作。為了更好地理解本模型參數選取的規則，文中以文獻［16］中的歐洲高山盆地簡化的非對稱場地模型為研究對象，采用文中方法模擬盆地地表反應進而闡述非線性相關參數選取的原則。如圖6所示，選取的計算模型尺寸總長約1 000 m，高約250 m，土層最大寬度約500 m，沉積土層最深處約225 m，共7層，場地左側邊界平直，右側邊界簡化為圓弧。輸入地震動為Ricker子波，見圖7。文獻［16］中土層統一了密度，具體土體物理參數見表2，土體動三軸試驗結果參見圖8。依據現有這些常規數據，下面將簡介本算例模擬過程中需要但未給定的相關參數選取原則。

圖6 盆地模型圖［16］（單位：m）Fig.6 Basin model diagram［16］（Unit：m）

圖7 地震輸入Ricker波時程［16］Fig.7 Time histories of seismic input of Ricker wave［16］

表2 土體的物理參數［16］Table 2 Physical parameters of the soil layers［16］

3.1 動剪切模量的確定

在非線性問題中剪切模量是隨著加載強度而變化的物理量。確定動剪切模量一般有2種方法：一是通過實驗直接得到，如文獻［13］的做法；另一種通過現場測量剪切波速及動三軸給出的剛度曲線參數換算得來［23］。因文中選取的盆地分析模型已提供了土體的物理參數及模量比曲線，故文中采用換算的方法計算動剪切模量。根據文獻［23］，土體在低應變情況下的最大剪切模量Gmax由式（15）計算：

式中：νs為土體的剪切波速；ρ為土體密度。根據文［16］給出的盆地土層剪切模量比曲線（圖8），當前應變狀態下的動剪切模量G的計算公式為：

式中：γ為剪應變；f(γ)為不同剪應變下G/Gmax的比值。

3.2 抗剪強度的選取原則

抗剪強度也是模擬土體非線性特征的重要參數之一，其值與動剪切模量直接相關。由上節公式，抗剪強度τ可通過換算得到，計算公式如下：

根據式（17）將圖8的剪切模量比曲線轉化為圖9的應力應變曲線，再通過圖9的應力應變曲線可以直接獲得本模型每層土體的最大抗剪強度。

圖8 模量比曲線Fig.8 Modulus ratio curve

圖9 應力應變曲線Fig.9 Stress-stain curv

3.3 屈服面系數及屈服面數量的選取

文中采用的非線性本構模型屬于彈塑性多屈面模型，因此計算中需要給出與屈服面相關的分析參數，主要包括屈服面系數和屈服面數量，其中屈服面系數為常數，可以通過室內實驗直接獲取，也可通過遞推公式計算［13］給出：

式中：cj表示當前待計算的屈服面系數；cn表示以往已計算的屈服面系數；e表示該屈服面上剪應變值的歸一值；k表示該屈服面上切應力值的歸一值，下標j及j+1分別表示當前及下一個屈服面。

因土體具有較強非線性特征，即在較弱荷載作用下土體仍會進入塑性，故其彈性域幾乎可以忽略不計，但文中本構模型仍然保留了較小的彈性域，通過增加屈服面數量達到減小彈性域的目的。文中假設僅取第一個屈服面的范圍為彈性域，隨著屈服面數量的增多、第一個屈服面的范圍越來越小，彈性區間也相應變小，塑性區間相應增大。為了討論屈服面數量合理取值范圍，文中通過將屈服面數量（以下用n表示）分別取為10、20、50、100，以Ricker子波為基底地震輸入（圖7），輸入加速度時程的峰值（PGA）分別取為0.5、2.5、5 m/s2進行模擬，分別計算地表3個監測點（圖6）處的加速度時程反應及其傅氏譜，以n=100計算結果為標準，由式（19）計算地表1～3監測點曲線的均方根誤差，計算所得的地表時程、傅氏譜誤差分別見圖10和圖11。

圖10 地表加速度時程誤差Fig.10 Errors of acceleration time histories at the points prescribed

圖11 地表傅氏譜誤差Fig.11 Errors of Fourier spectra at the points prescribed

式中：RMES為所求均方根誤差；n代表屈服面數；m代表加速度時程或傅氏譜的數據點數。

從圖10及圖11的對比結果來看，屈服面數量的差異對地表地震反應的模擬結果具有一定的影響。總體來說屈服面數量越多，計算精度越好。從圖10及圖11反映的誤差分析結果可以看出n=50時計算結果良好，而n=10時模擬結果與n=100的結果還是有著有意義的差異，n=20相對于n=10的結果與n=100的時程結果吻合度尚可。從兩圖中可以看出3種強度下均體現類似的規律，即n=50誤差非常小，可以忽略。n=20的誤差也是可以接受的，而n=10時無論模擬的時程還是傅氏譜，均存在較大的差異，不建議選用。通過以上分析可以看出，屈服面數量選取應大于20。考慮到計算效率與計算精度的匹配及數值模擬結果，建議屈服面數量取為20，此時屈服面數量對地震反應分析的差異影響變小。

3.4 數值對比分析

GéLIS在文獻［16］中采用相似的彈塑性本構模型模擬了圖6描述的盆地，基于GéLIS在文獻［16］中給出的參數及計算模型（圖6），采用文中方法對該場地模型進行模擬。通過與GéLIS的模擬結果進行對比，驗證文中分析過程的可靠性。

根據3.1及3.2節的說明，計算本模型的動剪切模量，抗剪強度。根據3.3節的討論，本計算過程屈服面數量取為20，屈服面系數由式（18）計算，阻尼系數由式（13）計算。圖12及圖13分別給出了地表3個監測點的速度時程及其反應譜的模擬結果，為了對比，圖中還給出了文獻［16］的對應解。

圖12 監測點（圖6所示）加速度時程Fig.12 Acceleration time histories of the output points prescribed in Fig.6

從圖12中的結果可以看出文中給出的3個輸出點的時程反應與文獻［16］的結果基本吻合，無論是時程曲線波形還是峰值均比較接近參考解，只是在監測點1處的時程在1.5～2 s之間與參考解稍微不同，這個特征主要在于對盆地產生的面波反應分析存在稍許差異，可能是阻尼參數選取差異的影響。圖13計算的傅氏譜存在稍許差異，監測點2及監測點3的特征周期值與參考解基本吻合，特征周期處的傅氏譜峰值稍有差異，這可能也源于文中采用了與文獻［16］不同的阻尼模型假設。對監測點1，文中解與參考解存在一定的差異，主要體現在傅氏譜峰值對應的特征周期相對于文獻［16］解略小。原因可能是分析模型介質參數存在一定差異或輸出位置不同而產生的影響。通過對時程反應和傅氏譜模擬結果總體分析，文中解與文獻參考解的計算規律基本吻合，說明文中計算過程及分析參數選取方法是基本合理的，應用文中計算過程分析二維復雜土層地震反應是可靠的。

圖13 監測點（圖6所示）加速度時程傅氏譜Fig.13 Fourier spectra of accelerations on the points prescribed in Fig.6

4 結論

文中在通用軟件ABAQUS的基礎上，建立了適用于包含復雜土層場地的二維地震反應分析過程，使其可以同時考慮土體結構空間的復雜性及土體的非線性及粘性特征。通過對本分析過程的論證，可以得出以下一些結論：

（1）文中設置的粘性邊界，可以使外傳波成功通過人工邊界，保證計算的合理性。

（2）文中選用的非線性動力本構模型適用于土體二維地震反應分析，能合理模擬土體的非線性特征。

（3）相較于現有二維動力本構模型，文中選取的動力非線性本構模型具有參數簡單，易選取的特征，可以利用等效線性化模型參數直接轉化，也可以通過土樣試驗直接獲得。

（4）文中計算過程是基于ABAQUS的顯式計算過程，具有良好的計算效率，并可以實現較大的計算規模，為大尺度非線性地震反應分析過程的實現提供了借鑒。