基于粘彈性邊界的三維土體地震響應分析

吳 翔，丁海平

（蘇州科技大學 江蘇省結構工程重點實驗室，江蘇 蘇州 215011）

工程場地的土層地震反應分析是地震安全性評價工作中的重要組成部分，在進行場地地震反應分析時，土體的動力本構關系、模型尺寸和網格劃分、邊界條件，以及地震動輸入等都是需要考慮的問題。

在通常的土層地震響應分析中，常采用基于一維剪切土層模型等效線性化的波動頻域分析方法[1]；對地表地形起伏變化的場地，一般采用二維場地模型的時域線性分析。事實上，土體的非線性對地震響應影響很大。Drucker 和Prager 于1952 年提出了考慮靜水壓力影響的廣義Mises 屈服與破壞準則，被稱為Drucker-Prager（DP）準則[2]，后經Sandler、Krajcinovic 等人[3-4]修正后被廣泛使用。隨著對DP 模型研究的不斷深入[5-8]，更符合巖土材料特性的Extended Drucker-Prager（EDP）模型開始被廣泛使用[9]。宋宏偉、Luo Chuan 等[10-11]研究了DP 模型和EDP 模型的特點及相互轉換問題，證明了EDP 模型的可行性且更符合巖土材料的實際情況。事實上，土體的非線性對地震響應影響很大。巴振寧、梁建文和齊文浩等[12-15]分別對考慮土體線性與非線性、等效線性化與非線性在地震響應中的差異進行了研究分析。本文將依托于大型有限元分析軟件ANSYS，以實際工程為背景，建立三維成層土體非線性地震響應分析模型，借助DEEPSOIL 程序實現土體的等效線性化。將對線性模型、等效線性化模型和非線性（EDP）模型的地震響應進行分析，研究土體非線性特性對地震響應的影響。

1 粘彈性邊界

在考慮土體的地震響應分析中，核心問題之一是如何正確地模擬無限地基，如采用合適的人工邊界。常用的有：無窮邊界元法[16]、阻尼提取法[17]和人工透射邊界法[18]。本文選擇的是粘彈性人工邊界[19]。粘彈性人工邊界的實質是通過把彈簧和阻尼器設置在每個邊界節點上來模擬實際的邊界條件。Deeks、劉晶波等[19-20]首先在二維散射波為柱面波的假設下建立了二維粘彈性人工邊界，并分析了該邊界相對于粘性邊界的優勢。劉晶波等[21]基于球面波動理論，將粘彈性人工邊界由二維平面推廣到三維，并通過算例驗證了三維粘彈性人工邊界的準確性和適用性。本文采用的粘彈性邊界參數公式見式（1）。圖1 為邊界節點三向彈簧-阻尼系統示意圖。

圖1 邊界節點三向彈簧-阻尼系統示意圖

其中：KBT和KBN分別為切向彈簧和法向彈簧的剛度系數；CBT和CBN分別為切向阻尼器和法向阻尼器的阻尼系數；R 為波源至各個邊界節點的距離；cP和cS為介質的壓縮波速和剪切波速；G 和ρ 分別為土層的剪切模量和密度；αT與αN分別為彈簧切向與法向剛度修正系數，參考劉晶波等[22]研究，取值范圍及推薦值見表1，本文分別取值為αT=0.67，αN=1.33。

表1 粘彈性人工邊界修正系數

通過將地震荷載轉化為切向和法向等效節點力的方式實現粘彈性人工邊界的模擬，本文施加的等效節點力采用自由場應力法進行推導。具體公式如式（2）所列。

其中，σi（t）為邊界上一點i 在t 時刻的等效荷載；σf（x，y，t）為由原自由場產生的應力；k 和c 分別為彈簧剛度和阻尼系數；（x，y，t）和uf（x，y，t）分別為i 點的自由場速度與位移。邊界點上施加的等效力通過式（3）計算。

式中，A 為邊界點所代表的有效面積。采用公式（3）計算邊界節點的等效荷載力。

2 粘彈性邊界驗證

采用圖2 所示的土體模型為算例。模型尺寸為20 m×20 m×10 m，有限元網格尺寸為0.5 m，土體密度為2 000 kg/m3，剪切波速為200 m/s，壓縮波速為346.41 m/s，泊松比為0.25。入射波位移函數如式（4）所示；圖3為相應的位移時程。由波動理論可知，對于圖2 的彈性半空間模型，地表位移的理論值為輸入位移的兩倍。由圖4 可見，地表位移響應的理論解與數值解吻合很好，說明本文采用的粘彈性邊界精度良好，可以使用。

圖2 三維土體及粘彈性邊界有限元模型

圖3 入射波位移時程

圖4 地表位移響應的理論解與數值解對比

3 經典DP 模型與擴展DP 模型

3.1 經典DP 模型

在大型有限元分析軟件ANSYS 中，DP 模型廣泛用于巖石、土體等材料分析中。在使用模擬單元Plane42進行二維實體模擬分析，以及單元Solid45 單元進行三維實體分析的情況下，可以使用經典DP 模型。ANSYS中設定經典DP 模型需要輸入3 個參數，即粘聚力、內摩擦角、膨脹角。其中的膨脹角是用來控制體積膨脹的大小的。在巖土工程中，一般密實的砂土和超強固結土在發生剪切的時候會出現體積膨脹，因為顆粒重新排列了；而一般的砂土或者正常固結的土體，只會發生剪縮。所以在使用DP 模型的時候，對于一般的土，膨脹角設置為0°是比較符合實際的。

3.2 擴展DP 模型

EDP 模型作為DP 模型的改進模型[10]，考慮材料的硬化特性，更符合巖土材料的實際情況。在新版本的ANSYS 中，經常使用Solid185 單元替代Solid45 單元，經典DP 模型不能應用在Solid185 單元中，而要采用擴展DP 模型，即EDP 模型。經典DP 模型到EDP 模型的參數需要轉換，即通過已知的粘聚力和內摩擦角等計算得到EDP 所需要的參數數值。通過對比經典DP 模型和EDP 模型的屈服函數，發現EDP 模型中的線性屈服函數與DP 模型的屈服函數形式上相似且屈服面形狀也相同，通過參數等效替換，EDP 參數的計算公式見下式。

其中，C1為EDP 模型第一個參數，即應力敏感度；C2為EDP 模型的第二個參數，即屈服強度；? 代表土體內摩擦角，C 代表土體粘聚力。

本文使用Solid185 三維實體單元模擬土體，故在考慮土體非線性時，對土體施加EDP 模型。

4 計算模型

4.1 場地條件及土體網格劃分

本文選用奧林匹克塔（位于北京市奧林匹克公園）所在位置的場地作為研究對象。土體1/4 模型如圖5 所示，尺寸為160 m×160 m×50 m。由下至上依次選取4 個監測點，監測點為每層土體平面中間點，具體坐標如下：A（80，80，0），B（80，80，20），C（80，80，30）和D（80，80，50）)。場地土參數如表2 所示。

圖5 土體1/4 模型

表2 場地土參數

在有限元數值模擬中，計算網格的劃分對計算結果有明顯的影響。網格尺寸過小將導致計算時間過長，同時還會導致計算結果不能收斂，而土體單元過大將導致計算精度不夠。因此，合理選擇是一個非常重要的問題。一般而言，土體單元高度應符合式（5）和式（6），土體單元的寬度通?？梢匀「叨鹊? 倍左右[23]

式中，λmin為地震波的最小波長（λmin=Vs/fmax；Vs為土體的剪切波速，fmax為最高頻率）。采用的有限元網格尺寸為2 m×2 m×2 m。

4.2 等效線性化模型

本文借助DEEPSOIL 程序，采用等效線性化方法建立等價非線性粘-彈塑性模型，用等效的剪切模量和阻尼比代替不同應變幅值下的剪切模量和阻尼比，將非線性問題轉化為線性問題，通過計算得到對應監測點的加速度時域及頻域結果。土類動剪切模量比和阻尼比采用袁曉銘[24]的推薦值，如表3 所示。

表3 土樣剪切模量比和阻尼比

4.3 地震波的選取

奧林匹克塔的設計基準期為50 年，所處地區的抗震設防烈度為8 度，設計地震分組為第一組，所在場地為Ⅱ類場地，基本地震加速度為0.20g[25]。

入射地震波選取EL-Centro 波，其加速度記錄如圖6。

圖6 EL-centro 地震波加速度記錄

5 計算結果及分析

5.1 位移響應分析

線性模型、非線性模型，以及等效線性化模型各個監測點分別在多遇地震作用下的水平位移響應見圖7，設防地震作用下位移響應見圖8，罕遇地震作用下位移響應見圖9。根據位移時程結果可以得出：（1）在多遇地震作用下，線性、非線性及等效線性化模型的位移響應結果較為接近。這是由于輸入地震動強度較小，土體沒有完全進入塑性階段。

圖7 多遇地震作用下各模型及監測點水平位移

圖8 設防地震作用下各模型及監測點水平位移

圖9 罕遇地震作用下各模型及監測點水平位移

（2）在設防地震和罕遇地震作用下，對于監測點A 點，線性模型、非線性模型和等效線性化模型的位移較為接近。監測點B 點，C 點和D 點均表現為：等效線性化的位移要明顯大于線性和非線性模型；線性模型的位移大于非線性模型的位移。隨著土的埋深變淺和地震動強度增加，非線性模型水平位移響應的延遲愈加明顯。這是由于土體非線性產生的變形及土體阻尼的增加耗散了一定的地震動能量。

（3）從圖7 至圖9 可以看出，在多遇地震作用下，線性、非線性和等效線性化的計算結果相似；在設防地震和罕遇地震作用下，等效線性化得到的位移響應顯著大于線性和非線性模型，這是等效線性化在地震動輸入較大的情況下過高估計了場地土的非線性特性所致。

5.2 加速度時程分析

多遇地震、設防地震和罕遇地震各模型及監測點的加速度時域分析見圖10 至圖12；為便于對比分析，加速度時域分析取前10 s 的計算結果。從加速度時域結果可以得出：

圖10 多遇地震作用下各模型及監測點水平加速度

圖11 設防地震作用下各模型及監測點水平加速度

圖12 罕遇地震作用下各模型及監測點水平加速度

（1）在多遇地震作用下，線性模型、非線性模型及等效線性化模型的加速度響應總體趨勢相同；在加速度峰值方面，非線性模型和等效線性化模型較為接近，二者均小于線性工況下的加速度峰值，這是由于考慮土體非線性后，土體產生的變形及土體阻尼的增加耗散了一定的地震動能量。由于輸入地震動強度較小，土體沒有完全進入塑性階段，三種工況下整體區別不大。

（2）在設防地震和罕遇地震作用下，線性條件下的加速度峰值大于非線性條件下的加速度峰值，且隨著土的埋深變淺和地震動強度的增大，二者的差值愈加明顯，等效線性化的加速度峰值大于非線性條件下加速度峰值。對比多遇地震可見：非線性和線性地震響應存在的差異隨著入射地震動幅值的增大而增大。這是由于土體非線性耗散了地震動能量。

（3）從圖10 至圖12 可以看出，等效線性化得到的加速度峰值與非線性較為接近，但隨著輸入地震動強度的增加，等效線性化的加速度峰值大于非線性。在輸入地震動強度較小的情況下，非線性和等效線性化的計算結果波形相似，在地震動輸入強度變大的情況下，二者的差別較大，等效線性化所得的波形變得稀疏。該結論與文獻[15]所得結論相同。

5.3 加速度傅里葉譜分析

加速度傅里葉譜見圖13 至圖15。從加速度頻域結果分析可以得出：

圖13 多遇地震作用下各模型及監測點水平加速度傅里葉譜

圖14 設防地震作用下各模型及監測點水平加速度傅里葉譜

圖15 罕遇地震作用下各模型及監測點水平加速度傅里葉譜

（1）在多遇地震作用下，由于土體的濾波效應，非線性和等效線性化的傅里葉譜幅值小于線性模型。等效線性化條件下土體基頻降低，圖形向左移。由于地震動強度較小，土體沒有完全進入非線性階段，線性、非線性和等效線性化的加速度傅里葉譜總體相差不大。

（2）在設防地震作用下，等效線性化的傅氏譜峰值小于非線性模型，非線性模型的傅里葉譜幅值小于線性模型。在罕遇地震作用下，等效線性化和非線性模型的傅里葉譜幅值遠小于線性模型。這是由于土體非線性時模量衰減，基頻降低。由于土體的濾波效應，隨著埋深變淺，非線性時中高頻成分被土體吸收。監測點C和監測點D 非線性和等效線性化條件下的低頻顯著降低。這一現象是由于地震動輸入增大，場地土非線性特性影響增強的結果。這現象與文獻[15]所描述規律相似。

（3）等效線性化條件下，土體模量衰減，剛度降低，基頻降低，傅氏譜圖形向左移了。隨著地震動動輸入增大，中高頻部分快速降低，這是等效線性化在地震動輸入較大的情況下過高估計了場地土的非線性特性所致。

6 結論

基于大型有限元分析軟件ANSYS 平臺及其APDL 編程語言編制了相應的程序，使用EDP 模型考慮土體的非線性特性，施加粘彈性人工邊界并驗證了其正確性，然后以某實際工程場地為背景對其進行了計算分析，通過與DEEPSOIL 程序得出的等效線性化計算結果進行對比分析，探討了土體非線性特性對其自身地震響應的影響?；诒疚牡挠嬎惴治?，可得如下結論:

（1）當輸入地震動強度較小且土體埋深較大時，土的非線性影響不夠明顯，線性、非線性及等效線性化計算結果相差較小。隨著埋深變淺和地震動強度的增大，非線性的位移響應大于線性情況；等效線性化的位移響應顯著大于線性和非線性模型。

（2）考慮土體非線性后，土體產生的變形及土體阻尼的增加耗散了地震動能量，隨著埋深變淺和輸入地震動強度的增加，線性條件下的加速度峰值大于非線性條件下的加速度峰值，非線性條件下的加速度峰值大于等效線性化條件下的加速度峰值，線性計算結果的加速度放大倍數保持不變，非線性計算結果的加速度放大倍數呈遞減趨勢。能反應“隨輸入地震動強度增大，土層非線性效應表現越明顯”這一基本規律[25]。

（3）在EL-Centro 波作用下，隨著埋深變淺和輸入地震動強度的增加，土體非線性時模量衰減，剛度降低，低頻部分顯著降低，同時由于土體的濾波效應，非線性時中高頻成分被土體吸收。等效線性化方法雖簡單方便，但其往往會過高估計場地土的非線性，因而在地震動輸入較大的情況下須謹慎使用。

由此可見，土體非線性特性對其自身地震響應的影響十分顯著，在實際工程中使用非線性模型應當考慮土體非線性的影響。