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

楊笑梅，陳鑫，游昊冉

（廣東工業(yè)大學(xué)土木與交通工程學(xué)院，廣東 廣州 510006）

引言

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

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

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

1 2D P－SV波傳播模擬在ABAQUS中的實(shí)現(xiàn)

2D場(chǎng)地地震反應(yīng)問(wèn)題理論上屬于無(wú)限域二維P－SV波傳播問(wèn)題，通用的有限元軟件無(wú)法直接模擬。需要解決人工邊界設(shè)置、土體的非線性本構(gòu)模型及土體的阻尼模型建立等相關(guān)問(wèn)題。為了更好地闡述文中構(gòu)建的計(jì)算過(guò)程，本節(jié)將簡(jiǎn)要介紹其相關(guān)理論及其在ABAQUS中的實(shí)現(xiàn)。

1.1 粘性人工邊界條件及地震輸入在ABAQUS中的實(shí)現(xiàn)

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

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

地震輸入以等效結(jié)點(diǎn)力F實(shí)現(xiàn)［18－19］，計(jì)算公式如式（2）：

式中：下標(biāo)x或y表示施加于結(jié)點(diǎn)處力的方向；上標(biāo)x或y表示外法線方向；Δt1及Δt2分別表示邊界處入射波和反射波的延遲時(shí)間。

1.2 多屈服面本構(gòu)模型理論簡(jiǎn)介

文中建立的本構(gòu)模型是基于Joyner［13］在Iwan［10］模型基礎(chǔ)上改進(jìn)的二維多屈面模型，由于該模型對(duì)應(yīng)力路徑能夠合理記憶，因此在任意復(fù)雜的往復(fù)荷載作用下模擬其加卸載曲線均能達(dá)到較好的結(jié)果。為理解文中計(jì)算過(guò)程，先簡(jiǎn)要介紹文獻(xiàn)［13］中的多屈面本構(gòu)模型。

1.2.1 相關(guān)應(yīng)力應(yīng)變量的定義

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

式中：Sij和Eij分別表示總應(yīng)力及總應(yīng)變；δij為Kronecher符號(hào)。由式（6）及式（7）可知，總應(yīng)力或總應(yīng)變量是通過(guò)計(jì)算出平均應(yīng)力或應(yīng)變及偏應(yīng)力或應(yīng)變后求得。模型假設(shè)塑性應(yīng)變?nèi)际怯善珣?yīng)力產(chǎn)生，平均應(yīng)力及平均應(yīng)變關(guān)系為線彈性。

1.2.2 多屈面模型的假設(shè)

文中彈塑性模型假設(shè)主要采用文獻(xiàn)［13］的理論確定偏應(yīng)力應(yīng)變?cè)隽筷P(guān)系，該模型的屈服準(zhǔn)則、加卸載準(zhǔn)則、移動(dòng)硬化定律以及流動(dòng)法則描述如下：

（1）屈服準(zhǔn)則

土體在偏應(yīng)力空間的應(yīng)力狀態(tài)下，假設(shè)有一系列的屈服面，滿足如下屈服函數(shù)：

式中：F是屈服函數(shù)；αij表示移動(dòng)硬化參數(shù)；k表示與屈服強(qiáng)度相關(guān)的常數(shù)；下標(biāo)n代表第n個(gè)屈服面。

（2）移動(dòng)硬化準(zhǔn)則

為了體現(xiàn)動(dòng)力作用，Joyner采用式（10）作為移動(dòng)硬化定律，從而保證應(yīng)力點(diǎn)和屈服面的一致性：

式中：上標(biāo)含“p”的物理量表示前一時(shí)刻的計(jì)算值，上標(biāo)含“p+1”的物理量表示當(dāng)前時(shí)間步的計(jì)算值。

（3）流動(dòng)法則

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

這里新出現(xiàn)的變量Ln與加卸載準(zhǔn)則相關(guān)，具體參數(shù)如下節(jié)，而hn為非負(fù)的塑性標(biāo)量因子，由式（12）計(jì)算：

式中，cn為第n個(gè)屈服面的屈服面系數(shù)。

（4）加卸載準(zhǔn)則

土體受到往復(fù)動(dòng)力荷載作用時(shí)，假設(shè)的加卸載準(zhǔn)則如下：

1.3 土體阻尼模型的選取

土體除了具有非線性特征還存在較強(qiáng)的吸能特征。選取在計(jì)算過(guò)程中保持定阻尼值的瑞利阻尼［20］考慮土體的耗能特征。參考文獻(xiàn)［21］給定的計(jì)算公式確定阻尼系數(shù)，即：

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

1.4 彈塑性模型VUMAT子程序的實(shí)現(xiàn)流程

利用了ABAQUS二次開(kāi)發(fā)編寫了對(duì)應(yīng)文中彈塑性本構(gòu)計(jì)算的VUMAT子程序，具體計(jì)算過(guò)程描述如下：（1）讀入由ABAQUS數(shù)據(jù)文件給定的初始總應(yīng)力、總應(yīng)變、初始狀態(tài)變量值，并根據(jù)式（6）及式（7）計(jì)算平均應(yīng)力及偏應(yīng)力，平均應(yīng)變及偏應(yīng)變初值；（2）調(diào)用VUMAT子程序，由ABAQUS中讀取的位移值，計(jì)算平均應(yīng)變?cè)隽考捌珣?yīng)變?cè)隽浚唬?）由平均應(yīng)變?cè)隽坑?jì)算平均應(yīng)力增量；（4）按彈性本構(gòu)模型計(jì)算試探應(yīng)力，并根據(jù)計(jì)算結(jié)果按屈服準(zhǔn)則即式（8）及式（9）判別是否屈服，若判別結(jié)果為彈性，則子程序自動(dòng)更新應(yīng)力及應(yīng)變值進(jìn)行下一個(gè)增量步的計(jì)算；（5）若判別為塑性，按塑性理論要求計(jì)算，遵循文中的移動(dòng)硬化準(zhǔn)則、流動(dòng)法則及加卸載準(zhǔn)則由偏應(yīng)變?cè)隽坑?jì)算偏應(yīng)力增量；（6）由計(jì)算所得的平均應(yīng)力增量及偏應(yīng)力增量，更新偏應(yīng)力及總應(yīng)力，并由平均應(yīng)變?cè)隽考捌珣?yīng)變?cè)隽扛缕珣?yīng)變及總應(yīng)變值；（7）重復(fù)（2）～（4）步直至計(jì)算停止。

2 基于ABAQUS數(shù)值模擬過(guò)程P－SV波傳播問(wèn)題的驗(yàn)證

為了驗(yàn)證文中計(jì)算過(guò)程的正確性，本節(jié)通過(guò)2個(gè)算例分別驗(yàn)證二維波傳播問(wèn)題邊界設(shè)置的合理性及多屈面本構(gòu)模型的可靠性。

2.1 半圓峽谷地形地表反應(yīng)

為了驗(yàn)證文中方法與計(jì)算的正確性，選取Kawase在文獻(xiàn)［22］中基于離散波速邊界元法分析的無(wú)限域半圓峽谷地表反應(yīng)問(wèn)題進(jìn)行分析，為了便于結(jié)果的比較，類似文獻(xiàn)［22］，計(jì)算時(shí)定義一個(gè)與頻率相關(guān)的參數(shù)：

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

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

從圖2結(jié)果對(duì)比可見(jiàn)。本文解與文獻(xiàn)［22］計(jì)算結(jié)果比較吻合，邊界反射的誤差波未對(duì)結(jié)果產(chǎn)生影響，也說(shuō)明地震波可以由邊界處外傳出計(jì)算域，x表示該點(diǎn)到谷地中央的距離，x/a為無(wú)量綱物理量。因此文中基于ABAQUS二維顯式有限元結(jié)合設(shè)置人工邊界及地震輸入的計(jì)算過(guò)程可以合理模擬二維無(wú)限域線性地震反應(yīng)。

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

2.2 二維土體非線性算例

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

圖3 土壩簡(jiǎn)化FEM模型［13］Fig.3 Simplified FEM model for earth dam［13］

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

表1 土體參數(shù)［13］Table 1 Soil parameters［13］

圖4 基底輸入地震動(dòng)時(shí)程［13］Fig.4 Time history of base-input ground motion［13］

3 非線性反應(yīng)分析參數(shù)設(shè)置問(wèn)題的討論

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

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

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

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

3.1 動(dòng)剪切模量的確定

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

式中：νs為土體的剪切波速；ρ為土體密度。根據(jù)文［16］給出的盆地土層剪切模量比曲線（圖8），當(dāng)前應(yīng)變狀態(tài)下的動(dòng)剪切模量G的計(jì)算公式為：

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

3.2 抗剪強(qiáng)度的選取原則

抗剪強(qiáng)度也是模擬土體非線性特征的重要參數(shù)之一，其值與動(dòng)剪切模量直接相關(guān)。由上節(jié)公式，抗剪強(qiáng)度τ可通過(guò)換算得到，計(jì)算公式如下：

根據(jù)式（17）將圖8的剪切模量比曲線轉(zhuǎn)化為圖9的應(yīng)力應(yīng)變曲線，再通過(guò)圖9的應(yīng)力應(yīng)變曲線可以直接獲得本模型每層土體的最大抗剪強(qiáng)度。

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

圖9 應(yīng)力應(yīng)變曲線Fig.9 Stress-stain curv

3.3 屈服面系數(shù)及屈服面數(shù)量的選取

文中采用的非線性本構(gòu)模型屬于彈塑性多屈面模型，因此計(jì)算中需要給出與屈服面相關(guān)的分析參數(shù)，主要包括屈服面系數(shù)和屈服面數(shù)量，其中屈服面系數(shù)為常數(shù)，可以通過(guò)室內(nèi)實(shí)驗(yàn)直接獲取，也可通過(guò)遞推公式計(jì)算［13］給出：

式中：cj表示當(dāng)前待計(jì)算的屈服面系數(shù)；cn表示以往已計(jì)算的屈服面系數(shù)；e表示該屈服面上剪應(yīng)變值的歸一值；k表示該屈服面上切應(yīng)力值的歸一值，下標(biāo)j及j+1分別表示當(dāng)前及下一個(gè)屈服面。

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

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

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

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

從圖10及圖11的對(duì)比結(jié)果來(lái)看，屈服面數(shù)量的差異對(duì)地表地震反應(yīng)的模擬結(jié)果具有一定的影響?？傮w來(lái)說(shuō)屈服面數(shù)量越多，計(jì)算精度越好。從圖10及圖11反映的誤差分析結(jié)果可以看出n=50時(shí)計(jì)算結(jié)果良好，而n=10時(shí)模擬結(jié)果與n=100的結(jié)果還是有著有意義的差異，n=20相對(duì)于n=10的結(jié)果與n=100的時(shí)程結(jié)果吻合度尚可。從兩圖中可以看出3種強(qiáng)度下均體現(xiàn)類似的規(guī)律，即n=50誤差非常小，可以忽略。n=20的誤差也是可以接受的，而n=10時(shí)無(wú)論模擬的時(shí)程還是傅氏譜，均存在較大的差異，不建議選用。通過(guò)以上分析可以看出，屈服面數(shù)量選取應(yīng)大于20?？紤]到計(jì)算效率與計(jì)算精度的匹配及數(shù)值模擬結(jié)果，建議屈服面數(shù)量取為20，此時(shí)屈服面數(shù)量對(duì)地震反應(yīng)分析的差異影響變小。

3.4 數(shù)值對(duì)比分析

GéLIS在文獻(xiàn)［16］中采用相似的彈塑性本構(gòu)模型模擬了圖6描述的盆地，基于GéLIS在文獻(xiàn)［16］中給出的參數(shù)及計(jì)算模型（圖6），采用文中方法對(duì)該場(chǎng)地模型進(jìn)行模擬。通過(guò)與GéLIS的模擬結(jié)果進(jìn)行對(duì)比，驗(yàn)證文中分析過(guò)程的可靠性。

根據(jù)3.1及3.2節(jié)的說(shuō)明，計(jì)算本模型的動(dòng)剪切模量，抗剪強(qiáng)度。根據(jù)3.3節(jié)的討論，本計(jì)算過(guò)程屈服面數(shù)量取為20，屈服面系數(shù)由式（18）計(jì)算，阻尼系數(shù)由式（13）計(jì)算。圖12及圖13分別給出了地表3個(gè)監(jiān)測(cè)點(diǎn)的速度時(shí)程及其反應(yīng)譜的模擬結(jié)果，為了對(duì)比，圖中還給出了文獻(xiàn)［16］的對(duì)應(yīng)解。

圖12 監(jiān)測(cè)點(diǎn)（圖6所示）加速度時(shí)程Fig.12 Acceleration time histories of the output points prescribed in Fig.6

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

圖13 監(jiān)測(cè)點(diǎn)（圖6所示）加速度時(shí)程傅氏譜Fig.13 Fourier spectra of accelerations on the points prescribed in Fig.6

4 結(jié)論

文中在通用軟件ABAQUS的基礎(chǔ)上，建立了適用于包含復(fù)雜土層場(chǎng)地的二維地震反應(yīng)分析過(guò)程，使其可以同時(shí)考慮土體結(jié)構(gòu)空間的復(fù)雜性及土體的非線性及粘性特征。通過(guò)對(duì)本分析過(guò)程的論證，可以得出以下一些結(jié)論：

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

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

（3）相較于現(xiàn)有二維動(dòng)力本構(gòu)模型，文中選取的動(dòng)力非線性本構(gòu)模型具有參數(shù)簡(jiǎn)單，易選取的特征，可以利用等效線性化模型參數(shù)直接轉(zhuǎn)化，也可以通過(guò)土樣試驗(yàn)直接獲得。

（4）文中計(jì)算過(guò)程是基于ABAQUS的顯式計(jì)算過(guò)程，具有良好的計(jì)算效率，并可以實(shí)現(xiàn)較大的計(jì)算規(guī)模，為大尺度非線性地震反應(yīng)分析過(guò)程的實(shí)現(xiàn)提供了借鑒。