王競(jìng)雄 李鴻晶， 邢浩潔

1） 中國(guó)南京 211816 南京工業(yè)大學(xué)工程力學(xué)研究所

2） 中國(guó)北京 100081 中國(guó)地震局地球物理研究所

引言

覆蓋土層對(duì)地震波傳播特性具有重要影響（高武平等，2012；李平等，2012），通過土層地震反應(yīng)分析獲得工程場(chǎng)地的地面運(yùn)動(dòng)特征，為工程結(jié)構(gòu)抗震評(píng)估、設(shè)計(jì)和加固提供依據(jù).許多實(shí)際工程場(chǎng)地都可以被簡(jiǎn)化為水平成層場(chǎng)地模型，即假定覆蓋土層的力學(xué)性質(zhì)沿豎向呈水平成層變化，沿橫向均勻無限延伸，將地震激勵(lì)考慮為垂直向上入射的剪切波.按照此模型分析，土層地震反應(yīng)分析實(shí)際即為一維波動(dòng)問題.由于該模型形式簡(jiǎn)潔，物理意義明確，至今依然是場(chǎng)地地震反應(yīng)分析的一種重要途徑（Zalachoris，Rathje，2015）.

常用的一維土層地震反應(yīng)分析方法主要包括頻域和時(shí)域兩類方法.頻域法以等效線性化方法為代表，最早由Idriss 和Seed （1968）提出.該方法將不同應(yīng)變水平下的剪切模量和阻尼比以一個(gè)等效剪切模量和阻尼比代替，從而把非線性問題轉(zhuǎn)化為線性問題并在頻域內(nèi)求解.廖振鵬（1989）根據(jù)等效線性化方法開發(fā)出了土層地震反應(yīng)分析計(jì)算程序LSSRLI-1，在實(shí)踐中得到了廣泛認(rèn)可.近年來，袁曉銘等（2016）采用直頻法動(dòng)剪模量阻尼比求解技術(shù)，提出了新一代土層地震反應(yīng)分析方法，克服了傳統(tǒng)方法低估軟弱場(chǎng)地和深厚場(chǎng)地放大效應(yīng)的缺陷.孫銳和袁曉銘（2021）提出了全局等效剪應(yīng)變的概念和算法，建立了一種新的等效線性化分析方法.然而，等效線性化方法尚存在一些缺點(diǎn)，比如無法反映土層真實(shí)的受力狀態(tài)（王志良，韓清宇，1981；欒茂田，林皋，1992），且忽略了較多的高頻成分（齊文浩，薄景山，2007）.為了解決上述問題，發(fā)展出了土層地震反應(yīng)分析的時(shí)域方法，可以真實(shí)地反映土層在地震作用下的表現(xiàn).丁海平和周正華（1998）提出了一維土層地震反應(yīng)分析的時(shí)域有限元解法，并指出時(shí)域算法能夠達(dá)到與頻域算法相同的計(jì)算精度，但其計(jì)算耗時(shí)遠(yuǎn)少于頻域算法.然而傳統(tǒng)有限單元法由于形函數(shù)階次較低，通常需將每個(gè)波長(zhǎng)內(nèi)布置6—10 個(gè)單元（廖振鵬，2002）才能準(zhǔn)確刻畫出波場(chǎng)特征，在土層復(fù)雜或地震動(dòng)高頻成分較豐富的情況下可能會(huì)需要較大的計(jì)算工作量.對(duì)此，邢浩潔等（2017）提出了分析成層場(chǎng)地地震反應(yīng)的切比雪夫（Chebyshev）譜元法，僅需劃分少量單元即可獲得較高的計(jì)算精度.但其在計(jì)算切比雪夫譜單元質(zhì)量矩陣時(shí)沿用了傳統(tǒng)譜元法的做法，即利用切比雪夫多項(xiàng)式的性質(zhì)獲得單元質(zhì)量矩陣的解析解，由此導(dǎo)出的質(zhì)量矩陣為非對(duì)角形式，無法充分發(fā)揮顯式時(shí)間積分算法的高效率優(yōu)勢(shì).

本文在邢浩潔等（2017）的工作基礎(chǔ)上，擬通過節(jié)點(diǎn)積分法（Fried，Malkus，1975）建立集中質(zhì)量切比雪夫譜單元，解決傳統(tǒng)切比雪夫譜元法由于需要對(duì)質(zhì)量矩陣求逆而造成計(jì)算效率不高的問題，同時(shí)避免傳統(tǒng)的質(zhì)量集中方法的隨意性，如行和集中法（Zienkiewiczet al，2013）或?qū)窃胤糯蠓ǎ℉intonet al，1976）.在該集中質(zhì)量切比雪夫譜單元模型中嵌入多次透射人工邊界（multi-transmitting formula，縮寫為MTF）（Liao，Wong，1984），結(jié)合中心差分形成一種求解一維土層波動(dòng)問題的高階顯式算法.并利用日本Kik-net 強(qiáng)震觀測(cè)臺(tái)網(wǎng)提供的井下和地面實(shí)測(cè)記錄，以檢驗(yàn)本文方法的適用性.

1 水平成層場(chǎng)地的集中質(zhì)量切比雪夫譜元模型

1.1 水平成層場(chǎng)地簡(jiǎn)化模型

圖1 水平成層場(chǎng)地切比雪夫譜元模型Fig. 1 Chebyshev spectral element model of horizontal layered soil site

1.2 切比雪夫譜單元位移模式

一維土層分析的切比雪夫譜單元以 ［ ?1，1 ］ 區(qū)間內(nèi)不等間距分布的高斯-洛巴托-切比雪夫（Gauss-Lobatto-Chebyshev，縮寫為GLC）節(jié)點(diǎn)為參考單元節(jié)點(diǎn)，它們是第一類切比雪夫多項(xiàng)式的極值點(diǎn).對(duì)于n階切比雪夫譜單元，GLC 節(jié)點(diǎn)的位置為ξi＝?cos（iπ/n），i＝0，1，···，n.

設(shè)單元位移模式為

式中，Tk（ξ）為k階第一類切比雪夫多項(xiàng)式，其表達(dá)式為Tk（ξ）＝cos（kcos?1ξ），ci和ck為多項(xiàng)式系數(shù)，當(dāng)i＝0 或n時(shí)，取值為2，當(dāng)i＝1，···，n－1 時(shí)，取值為1.

1.3 切比雪夫譜單元特性矩陣

利用切比雪夫譜單元對(duì)水平成層場(chǎng)地模型進(jìn)行空間離散后，得到離散方程

式中ρ，η和G分別表示土體的質(zhì)量密度、阻尼系數(shù)和剪切模量， |J|為雅可比矩陣的行列式.對(duì)于一維譜單元，雅可比矩陣僅包含一個(gè)元素，即J＝Δz/2，Δz為物理單元的長(zhǎng)度.

為了建立集中質(zhì)量矩陣，本文采用節(jié)點(diǎn)積分法求解單元質(zhì)量矩陣，即以單元節(jié)點(diǎn)作為數(shù)值積分點(diǎn)，利用高斯-洛巴托（Gauss-Lobatto）數(shù)值積分計(jì)算單元質(zhì)量矩陣的積分（式（4））.由于拉格朗日形函數(shù)具有克羅內(nèi)克（Kronecker-δ）性質(zhì)，由此導(dǎo)出的質(zhì)量矩陣僅有主對(duì)角元素非零，其余非對(duì)角元素全部為零.單元GLC節(jié)點(diǎn)對(duì)應(yīng)的高斯-洛巴托積分權(quán)系數(shù)為

式中：wi為與節(jié)點(diǎn)ξi相對(duì)應(yīng)的積分權(quán)系數(shù)；φi（ξ）為式（2）中的單元形函數(shù).1—5 階切比雪夫譜單元的GLC 節(jié)點(diǎn)位置及對(duì)應(yīng)的積分權(quán)系數(shù)，列于表1.

表1 GLC 節(jié)點(diǎn)上的高斯?洛巴托積分權(quán)系數(shù)Table 1 Gauss-Lobatto quadrature weights based on GLC points

利用式（7）的數(shù)值積分格式計(jì)算式（4）單元質(zhì)量矩陣，則質(zhì)量矩陣的非對(duì)角元素全部為零，主對(duì)角元素可進(jìn)一步計(jì)算為

2 時(shí)域逐步積分

對(duì)于土層地震反應(yīng)分析問題，一般將入射波作為邊界處的運(yùn)動(dòng)形式施加，因此需要將模型內(nèi)域節(jié)點(diǎn)和邊界節(jié)點(diǎn)分開處理.在本文建立的水平成層場(chǎng)地譜元模型中，邊界節(jié)點(diǎn)僅有一個(gè)，即模型底部人工邊界上的單元節(jié)點(diǎn)，其余節(jié)點(diǎn)均屬于內(nèi)域節(jié)點(diǎn).將式（3）中的矩陣和向量按照內(nèi)域和邊界進(jìn)行分塊，得

一般情況下質(zhì)量矩陣MI為對(duì)角陣，而阻尼矩陣CI通常不會(huì)呈對(duì)角矩陣形式，故由MI與CI組合的矩陣亦會(huì)是非對(duì)角陣.此情形下可考慮構(gòu)建僅同質(zhì)量矩陣呈比例的瑞雷阻尼矩陣，或者將阻尼矩陣對(duì)角化（Thomsonet al，1974）等措施.

3 基于透射邊界的地震動(dòng)輸入

圖2 切比雪夫譜單元的一階MTF 插值方案Fig. 2 Interpolation scheme for first-order MTF in Chebyshev spectral element

4 算例分析

4.1 均勻半空間土體地震反應(yīng)分析

為了檢驗(yàn)本文方法的正確性，首先對(duì)一個(gè)均勻半空間場(chǎng)地模型進(jìn)行分析，根據(jù)行波理論得到其解析解.該模型土體厚度為180 m，土體質(zhì)量密度為2 000 kg/m3，剪切波速為250 m/s.底部設(shè)置一階MTF 人工邊界，垂直向上輸入幅值為1 m、主頻為2 Hz 的雷克（Ricker）波位移脈沖.將整個(gè)模型劃分為8 個(gè)4 階譜單元，則單元尺寸約為輸入波最短波長(zhǎng)的一半.圖3 顯示了土層底部和地表的位移反應(yīng)時(shí)程，其中底部反應(yīng)時(shí)程的第一個(gè)波峰為入射波，第二個(gè)波峰為地表反射波.由圖可知，入射波從底部傳至地表時(shí)間約為0.7 s，與解析解吻合，同時(shí)，地表位移反應(yīng)峰值等于入射波峰值的2 倍，符合自由地表?xiàng)l件.此外，地表未出現(xiàn)模型底部反射而來的地震波，說明MTF 人工邊界條件成功實(shí)現(xiàn).本例表明集中質(zhì)量切比雪夫譜元法能夠處理土層地震反應(yīng)分析問題，而且在網(wǎng)格較為稀疏的情況下依然能夠獲得較高的計(jì)算精度.

圖3 均勻半空間模型在雷克波入射下的地表和底部位移反應(yīng)Fig. 3 Displacement responses of ground surface and bottom of a homogeneous half-space model under Ricker wave incidence

4.2 Kik-net 臺(tái)站記錄模擬

利用日本Kik-net 強(qiáng)震觀測(cè)臺(tái)網(wǎng)（NIED，2021）提供的實(shí)測(cè)地震動(dòng)記錄和鉆孔數(shù)據(jù)檢驗(yàn)本文方法處理實(shí)際場(chǎng)地地震反應(yīng)分析的能力.從Kik-net 臺(tái)網(wǎng)中隨機(jī)選取4 個(gè)具有代表性的臺(tái)站，分別對(duì)應(yīng)中國(guó)《GB 50011—2010 建筑抗震設(shè)計(jì)規(guī)范》（中華人民共和國(guó)住房和城鄉(xiāng)建設(shè)部，中華人民共和國(guó)國(guó)家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局，2010）所定義的Ⅰ1，Ⅱ ，Ⅲ和Ⅳ等四類場(chǎng)地.在每個(gè)臺(tái)站中選擇實(shí)測(cè)地面峰值加速度（peak ground acceleration，縮寫為PGA）等于0.05—0.1g，0.1—0.2g和0.2—0.4g的EW 向或NS 向地震記錄各一條，分別對(duì)應(yīng)較弱地震動(dòng)、中等強(qiáng)度地震動(dòng)和較強(qiáng)地震動(dòng).以井下實(shí)測(cè)記錄作為一維水平成層場(chǎng)地模型底部的基巖輸入波，計(jì)算地表加速度反應(yīng)并與地表實(shí)測(cè)加速度記錄進(jìn)行對(duì)比.由于Kik-net 數(shù)據(jù)庫(kù)中未提供土體質(zhì)量密度數(shù)據(jù)，本文利用Boore （2016）提出的公式根據(jù)P 波波速和剪切波速估算土體密度.各臺(tái)站的土層剖面剪切波速如圖4 所示，按照《GB 50011—2010 建筑抗震設(shè)計(jì)規(guī)范》計(jì)算得到的基本信息列于表2.

圖4 四個(gè)Kik-net 臺(tái)站的土層剪切波速圖Fig. 4 Shear wave velocity profiles at the four Kik-net stations

表2 選用Kik-net 臺(tái)站的基本信息Table 2 Basic information of selected Kik-net stations

4 個(gè)臺(tái)站在不同強(qiáng)度地震作用下地表的加速度時(shí)程反應(yīng)，如圖5 所示.可見，在多數(shù)情況下，本文方法計(jì)算得到的地表加速度時(shí)程與實(shí)測(cè)的地表加速度記錄均能較好吻合.與輸入的基巖實(shí)測(cè)加速度時(shí)程相比，所有臺(tái)站的地表計(jì)算反應(yīng)均顯示出了土層的放大效應(yīng).即便對(duì)于IKRH02 臺(tái)站這樣覆蓋層厚度超過100 m 的深厚場(chǎng)地，也依然體現(xiàn)出了土層的放大效應(yīng)，未出現(xiàn)類似土層時(shí)域非線性軟件DEEPSOIL 嚴(yán)重低估地表反應(yīng)的問題（袁曉銘等，2016）.

圖5 不同強(qiáng)度地震動(dòng)作用下四個(gè)臺(tái)站地表加速度反應(yīng)時(shí)程（a） 較弱地震動(dòng)；（b） 中等強(qiáng)度地震動(dòng)；（c） 較強(qiáng)地震動(dòng)Fig. 5 Ground acceleration response histories of four stations under ground motions with different intensities（a） Weak ground motions；（b） Moderate ground motions； （c） Strong ground motions

表3 給出了按本文方法計(jì)算得到的PGA 與地震中實(shí)際記錄PGA 的對(duì)比.由表中數(shù)據(jù)可知，對(duì)于Ⅱ類場(chǎng)地（MYGH10），本文預(yù)測(cè)PGA 最接近實(shí)際值，誤差小于2.1%.但對(duì)于Ⅲ類場(chǎng)地（KMMH14），本文方法計(jì)算得到的PGA 大于實(shí)測(cè)值.總體而言，本文方法計(jì)算得到的PGA 總是稍小于或偏大于實(shí)際值，而未出現(xiàn)嚴(yán)重偏小的情況.表4 列出了本文計(jì)算得到的PGA 放大倍數(shù)及實(shí)際測(cè)得的PGA 放大倍數(shù).所選臺(tái)站的實(shí)測(cè)放大倍數(shù)介于3.28—7.40 之間不等，而計(jì)算放大倍數(shù)在3.67—10.48 之間不等.對(duì)于Ⅰ1類場(chǎng)地（TCGH08）、 Ⅱ 類場(chǎng)地（MYGH10）和Ⅳ類場(chǎng)地（IKRH02）在小震和中震情況下本文方法得到的PGA 放大倍數(shù)均與實(shí)測(cè)值相差不大，但在大震情況下本文結(jié)果偏大.

表3 不同強(qiáng)度地震動(dòng)作用下各臺(tái)站計(jì)算PGA 與實(shí)測(cè)PGA 對(duì)比Table 3 Comparison of computed PGA and recorded PGA for the stations under ground motions different intensities

表4 不同強(qiáng)度地震動(dòng)作用下各臺(tái)站PGA 放大倍數(shù)Table 4 Amplification factors of PGA for the stations under ground motions different intensities

圖6 為4 個(gè)臺(tái)站在不同強(qiáng)度的地震作用下的地表加速度反應(yīng)譜，阻尼比按5%計(jì)算.總體而言，本文計(jì)算反應(yīng)譜在低頻段（周期＞0.5 s）與實(shí)測(cè)記錄的反應(yīng)譜較為接近.觀察圖6 可知，在較弱的地震作用下，Ⅰ1和 Ⅱ 類場(chǎng)地的反應(yīng)譜峰值與實(shí)測(cè)記錄相近，而Ⅲ和Ⅳ類場(chǎng)地的反應(yīng)譜形狀與實(shí)測(cè)值吻合，其中Ⅳ類場(chǎng)地的計(jì)算反應(yīng)譜曲線與實(shí)測(cè)反應(yīng)譜十分接近.在中等強(qiáng)度地震作用下，Ⅰ1和Ⅳ類場(chǎng)地的反應(yīng)譜峰值與實(shí)測(cè)記錄較為吻合.在較強(qiáng)地震作用下，本文計(jì)算反應(yīng)譜曲線的形狀大致與實(shí)測(cè)記錄得到的反應(yīng)譜一致，但對(duì)于Ⅲ和Ⅳ類場(chǎng)地計(jì)算得到的反應(yīng)譜峰值比觀測(cè)結(jié)果偏大.

圖6 不同強(qiáng)度地震動(dòng)作用下四個(gè)臺(tái)站的地表加速度反應(yīng)譜（a） 較弱地震動(dòng)；（b） 中強(qiáng)地震動(dòng)；（c） 較強(qiáng)地震動(dòng)Fig. 6 Ground acceleration response spectra of four stations under ground motions with different intensities（a） Weak ground motions；（b） Moderate ground motion； （c） Strong ground motions

造成本文計(jì)算結(jié)果與實(shí)際臺(tái)陣記錄存在一定差異的原因可能有：① 本文采用的一維水平成層場(chǎng)地模型是一種高度簡(jiǎn)化的計(jì)算模型，無法反映出實(shí)際三維場(chǎng)地的地形特征，因此一維波動(dòng)的數(shù)值計(jì)算結(jié)果一般無法體現(xiàn)出對(duì)實(shí)際場(chǎng)地地震動(dòng)有重要影響的地形放大效應(yīng)；② 實(shí)際場(chǎng)地并不一定僅受剪切波作用，還有可能受到縱波、面波等多種地震波的作用，并且地震波的傳播方向也不一定恰好垂直向上的；③ 本文方法未考慮土體的非線性特性.

5 結(jié)論

本文提出了一種求解水平成層場(chǎng)地地震反應(yīng)的時(shí)域集中質(zhì)量切比雪夫譜元法.通過節(jié)點(diǎn)積分法嚴(yán)格地導(dǎo)出集中質(zhì)量矩陣，克服了傳統(tǒng)切比雪夫譜元法由于質(zhì)量矩陣為非對(duì)角矩陣形式所帶來的計(jì)算效率不高的問題.采用中心差分法進(jìn)行時(shí)間積分，并嵌入多次透射人工邊界，形成了高效的時(shí)域顯式算法.與傳統(tǒng)有限單元法相比，該方法在每個(gè)波長(zhǎng)內(nèi)僅需布置一個(gè)譜單元即可獲得令人滿意的計(jì)算精度，顯著降低了對(duì)空間網(wǎng)格尺寸的要求.

選擇日本Kik-net 強(qiáng)震臺(tái)網(wǎng)中分屬四種不同場(chǎng)地類型的臺(tái)站記錄檢驗(yàn)了本文方法處理實(shí)際場(chǎng)地的能力.計(jì)算結(jié)果表明，本文方法對(duì)于Ⅰ1，Ⅱ和Ⅳ類場(chǎng)地在較弱地震和中等強(qiáng)度地震作用下能夠較好地預(yù)測(cè)地面運(yùn)動(dòng)特征，但對(duì)于Ⅲ類場(chǎng)地及強(qiáng)震作用下的地表反應(yīng)計(jì)算結(jié)果與觀測(cè)結(jié)果相比仍存在一定誤差，后續(xù)工作中可考慮添加土體的時(shí)域非線性本構(gòu)關(guān)系.