王競雄 李鴻晶， 邢浩潔

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

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

引言

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

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

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

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

1.1 水平成層場地簡化模型

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

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

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

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

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

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

利用切比雪夫譜單元對水平成層場地模型進行空間離散后，得到離散方程

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

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

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

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

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

2 時域逐步積分

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

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

3 基于透射邊界的地震動輸入

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

4 算例分析

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

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

圖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 臺站記錄模擬

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

圖4 四個Kik-net 臺站的土層剪切波速圖Fig. 4 Shear wave velocity profiles at the four Kik-net stations

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

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

圖5 不同強度地震動作用下四個臺站地表加速度反應(yīng)時程（a） 較弱地震動；（b） 中等強度地震動；（c） 較強地震動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 給出了按本文方法計算得到的PGA 與地震中實際記錄PGA 的對比.由表中數(shù)據(jù)可知，對于Ⅱ類場地（MYGH10），本文預(yù)測PGA 最接近實際值，誤差小于2.1%.但對于Ⅲ類場地（KMMH14），本文方法計算得到的PGA 大于實測值.總體而言，本文方法計算得到的PGA 總是稍小于或偏大于實際值，而未出現(xiàn)嚴重偏小的情況.表4 列出了本文計算得到的PGA 放大倍數(shù)及實際測得的PGA 放大倍數(shù).所選臺站的實測放大倍數(shù)介于3.28—7.40 之間不等，而計算放大倍數(shù)在3.67—10.48 之間不等.對于Ⅰ1類場地（TCGH08）、 Ⅱ 類場地（MYGH10）和Ⅳ類場地（IKRH02）在小震和中震情況下本文方法得到的PGA 放大倍數(shù)均與實測值相差不大，但在大震情況下本文結(jié)果偏大.

表3 不同強度地震動作用下各臺站計算PGA 與實測PGA 對比Table 3 Comparison of computed PGA and recorded PGA for the stations under ground motions different intensities

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

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

圖6 不同強度地震動作用下四個臺站的地表加速度反應(yīng)譜（a） 較弱地震動；（b） 中強地震動；（c） 較強地震動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

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

5 結(jié)論

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

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