巴振寧，魯世斌，付繼賽，梁建文，蘆 燕

(1.中國(guó)地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點(diǎn)實(shí)驗(yàn)室，天津 300354；2.天津大學(xué)土木工程系，天津 300354)

準(zhǔn)確可靠的地震動(dòng)輸入是計(jì)算分析工程結(jié)構(gòu)地震反應(yīng)的前提[1]。目前工程應(yīng)用中大多采用地面運(yùn)動(dòng)預(yù)測(cè)方程(GMPEs)來(lái)獲取地震動(dòng)，但該方法基于遍歷假設(shè)，使得震源參數(shù)、傳播路徑和場(chǎng)地條件對(duì)地震動(dòng)時(shí)空分布的影響無(wú)法得到精確反應(yīng)[2]。因此，建立一種基于物理的全過(guò)程分析方法(包含震源破裂、地震波場(chǎng)傳播、土-結(jié)構(gòu)相互作用和結(jié)構(gòu)地震反應(yīng))，對(duì)精確結(jié)構(gòu)地震反應(yīng)分析具有重要意義。然而，由于全過(guò)程分析涉及地球物理、地震學(xué)、地震工程等多個(gè)學(xué)科，存在地殼層到工程結(jié)構(gòu)的多尺度跨越問(wèn)題，目前仍是一項(xiàng)難題。

近年來(lái)，諸多學(xué)者對(duì)全過(guò)程模擬進(jìn)行了有益探索，并提出一些全過(guò)程模擬方法，按照建模方法主要可分為3 類(lèi)：直接分析法、弱耦合兩步法和強(qiáng)耦合兩步法。其中，直接分析法是指將震源、傳播路徑和結(jié)構(gòu)建立在同一模型中直接進(jìn)行動(dòng)力反應(yīng)分析，例如：MAZZIERI 等[3]利用非連續(xù)伽遼金技術(shù)和譜元法，將區(qū)域尺度場(chǎng)地和工程結(jié)構(gòu)建立在同一模型中進(jìn)行分析，然而由于SEM中常用的顯式時(shí)間推進(jìn)的條件穩(wěn)定性，無(wú)法廣泛應(yīng)用于梁、板和墻等構(gòu)件的精細(xì)模擬；BRUN 等[4]利用混合異步時(shí)間積分器(HATI)，將SEM 和有限元(FEM)耦合在同一模型，并應(yīng)用于大壩非線(xiàn)性動(dòng)力分析，該方法可直接考慮從震源到結(jié)構(gòu)的全過(guò)程反應(yīng)，但由于數(shù)值方法的計(jì)算量隨頻率分辨率的16 倍增長(zhǎng)，需消耗大量計(jì)算資源；弱耦合兩步法是指先建立區(qū)域尺度模型求得地表地震動(dòng)，然后直接在結(jié)構(gòu)底部輸入所得地震動(dòng)進(jìn)行分析，例如：KRISHNAN 等[5]利用譜元法(SEM)求得南加州地區(qū)在兩次7.9 級(jí)設(shè)定地表作用下的地震波場(chǎng)，然后直接施加在一棟18 層鋼框架結(jié)構(gòu)底部進(jìn)行動(dòng)力分析，該方法操作流程較為簡(jiǎn)單，但未考慮土-結(jié)構(gòu)相互作用的影響；強(qiáng)耦合兩步法指將計(jì)算區(qū)域分為區(qū)域尺度模型(包含震源和傳播路徑)和局部區(qū)域模型(包含局部場(chǎng)地和結(jié)構(gòu))，首先對(duì)區(qū)域尺度模型進(jìn)行求解，然后利用區(qū)域縮減法將所求得的地震波場(chǎng)施加在局部區(qū)域的邊界，進(jìn)而對(duì)局部區(qū)域進(jìn)行動(dòng)力分析，例如：ICHINMURA 等[6-7]、MCCALLEN 等[8]、ZHANG 等[9]和 張 磊 等[10]學(xué)者，分別采用攝動(dòng)法、有限差分(FDM)、SEM、UCSB 等方法計(jì)算地震波場(chǎng)在區(qū)域尺度模型中的傳播，采用FEM 計(jì)算局部區(qū)域動(dòng)力反應(yīng)，該方法基于區(qū)域縮減[11-13]方法，相對(duì)于整體模擬法更為便捷，且具有可考慮土-結(jié)構(gòu)相互作用的優(yōu)點(diǎn)，應(yīng)用相對(duì)較廣，但計(jì)算效率仍受限于區(qū)域尺度模型的最高識(shí)別頻率和計(jì)算能力。

針對(duì)區(qū)域尺度模型中地震波場(chǎng)的傳播問(wèn)題，筆者前期建立了一種基于精確動(dòng)力剛度矩陣[14-15]的頻率波數(shù)域(FK)半解析方法[16-17]，該方法具有精度高、計(jì)算快和有效頻帶寬(可達(dá)0 Hz～20 Hz)等優(yōu)點(diǎn)，可有效克服上述計(jì)算瓶頸。因此基于以上現(xiàn)狀，本文以強(qiáng)耦合兩步法為框架，將FK 方法和FE 方法相結(jié)合，建立了從震源到工程結(jié)構(gòu)全過(guò)程地震反應(yīng)分析的FK-FE 混合方法。同時(shí)，為便于使用，筆者進(jìn)一步將FK 方法二次開(kāi)發(fā)至大型商用有限元分析軟件ABAQUS 中形成可視化插件，實(shí)現(xiàn)了僅在ABAQUS 平臺(tái)下即可完成位錯(cuò)點(diǎn)源作用下全過(guò)程地震反應(yīng)分析的功能，極大提高了建模和計(jì)算效率。

本文主要工作如下：首先，對(duì)FK 方法和FK-FE混合方法的理論進(jìn)行闡述；然后，對(duì)FK-FE 混合方法在ABAQUS 軟件中的主要操作步驟、二次開(kāi)發(fā)流程和插件功能進(jìn)行詳細(xì)介紹；在驗(yàn)證該混合方法正確性的基礎(chǔ)上，進(jìn)一步以Benchmark 框架模型(美國(guó)SAC Steel Project 中的9 層抗彎鋼框架結(jié)構(gòu))為例建立全過(guò)程分析模型，探討了地殼層速度結(jié)構(gòu)、和近地表場(chǎng)地條件對(duì)鋼框架結(jié)構(gòu)地震反應(yīng)的影響。

1 FK-FE 混合方法

圖1 給出了全過(guò)程地震反應(yīng)求解分析的示意圖，包括圖1(a)所示的位錯(cuò)點(diǎn)源-傳播路徑-建筑結(jié)構(gòu)整體模型和圖1(b)所示的求解步驟。具體求解時(shí)將整體模型劃分為STEP Ⅰ 所示的不含建筑結(jié)構(gòu)的層狀半無(wú)限空間模型和STEP Ⅲ 所示的有限元模型。

圖1 全過(guò)程結(jié)構(gòu)地震反應(yīng)求解分析示意圖Fig.1 Schematic of end-to-end structural seismic analysis

層狀半無(wú)限空間模型自下而上包括半無(wú)限空間、地殼層及近地表覆蓋土層，用于求解位錯(cuò)點(diǎn)源作用下層狀半無(wú)限空間的地震波傳播；有限元分析模型包括人工邊界子結(jié)構(gòu)、內(nèi)部域土體和建筑結(jié)構(gòu)，用于求解STEP Ⅰ 求得的地震波場(chǎng)輸入下建筑結(jié)構(gòu)的地震反應(yīng)，該模型可考慮土-結(jié)構(gòu)相互作用的影響。此外，在人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)設(shè)置黏彈性邊界用于吸收有限元模型中建筑結(jié)構(gòu)產(chǎn)生的外行散射波

本文建立FK-FE 混合方法的總體步驟如下：

STEP Ⅰ：采用FK 方法計(jì)算層狀半無(wú)限空間(不含結(jié)構(gòu))人工邊界子結(jié)構(gòu)內(nèi)外層節(jié)點(diǎn)位置位移時(shí)程。

STEP Ⅱ：將 STEP I 獲得的位移時(shí)程輸入人工邊界子結(jié)構(gòu)(殺死有限元模型內(nèi)部單元)，采用有限元方法求得人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)等效節(jié)點(diǎn)力時(shí)程。

STEP Ⅲ：將 STEP II 獲得的等效節(jié)點(diǎn)力時(shí)程施加在人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)，采用有限元方法對(duì)整個(gè)有限元模型(激活內(nèi)部單元)精細(xì)分析求得建筑結(jié)構(gòu)的地震反應(yīng)。

以下分別對(duì)FK 方法模擬位錯(cuò)點(diǎn)源作用下地震波場(chǎng)和FK-FE 混合方法實(shí)現(xiàn)原理進(jìn)行介紹。

1.1 FK 方法：位錯(cuò)點(diǎn)源作用下地震波場(chǎng)模擬

位錯(cuò)點(diǎn)源作用下場(chǎng)地某處位移響應(yīng)表示為[18]：

式中：ui為x點(diǎn)沿i方向的位移分量；t為時(shí)間；?為卷積運(yùn)算符；g(x,t)為位錯(cuò)點(diǎn)源作用下動(dòng)力格林函數(shù)；STF(t)為震源時(shí)間函數(shù)；Mpq為地震矩張量的六種位態(tài)分量(p,q=x,y,z)，在笛卡爾坐標(biāo)系中為一對(duì)沿p方向作用相反、沿q方向分開(kāi)的力偶。

地震矩張量位態(tài)分量Mpq可表示為與走向、傾角和滑動(dòng)角相關(guān)的形式[19]：

式中：M0為位錯(cuò)點(diǎn)源地震矩； ? 、 δ 和 λ分別為位錯(cuò)點(diǎn)源的走向、傾角和滑動(dòng)角。

位錯(cuò)點(diǎn)源格林函數(shù)的表達(dá)式為[17,20]：

式中： (Rkm,Skm,Tmk)為面諧基矢量；i 為虛數(shù)單位；k為水平波數(shù)；ω為圓頻率；下標(biāo)r、 θ和z為柱坐標(biāo)系下x點(diǎn)位移分量方向；m為方位角模數(shù)。詳細(xì)公式推導(dǎo)可參考文獻(xiàn)[17，20]。

式(3)中位移幅值分量ui可表示為：

式中：Jm為第一類(lèi)m階貝塞爾函數(shù)；um和wm分別為平面內(nèi)水平和豎向位移幅值；vm為平面外位移幅值；um、vm和wm具體表達(dá)式見(jiàn)文獻(xiàn)[17]。

1.2 人工邊界子結(jié)構(gòu)法：FK-FE 混合方法實(shí)現(xiàn)

FK-FE 混合方法的建立主要包括兩個(gè)關(guān)鍵技術(shù)：有限元模型地震波場(chǎng)輸入和有限元模型邊界處理。在地震波場(chǎng)輸入方面，采用人工邊界子結(jié)構(gòu)法[21-22]輸入FK 方法得到的地震波場(chǎng)；在邊界處理方面，采用三維黏彈性動(dòng)力人工邊界[23-24]來(lái)吸收有限元模型中建筑結(jié)構(gòu)產(chǎn)生的散射波。以下對(duì)人工邊界子結(jié)構(gòu)法和三維黏彈性動(dòng)力人工邊界進(jìn)行介紹。

1.2.1 人工邊界子結(jié)構(gòu)法

參考文獻(xiàn)[21]，輸入地震波場(chǎng)定義為自由波場(chǎng)，且可以轉(zhuǎn)化為在邊界上施加的等效節(jié)點(diǎn)力，等效節(jié)點(diǎn)力的確定只和自由波場(chǎng)(圖1 中STEP I)有關(guān)。因此，建立如圖2 所示的有限元模型和相對(duì)應(yīng)的自由場(chǎng)有限元模型。將模型按節(jié)點(diǎn)位置劃分為人工邊界子結(jié)構(gòu)內(nèi)、外層節(jié)點(diǎn)和剩余內(nèi)部節(jié)點(diǎn)。

圖2 人工邊界子結(jié)構(gòu)法求解示意Fig.2 Schematic for the substructure of artificial boundary method solution

按照節(jié)點(diǎn)位置，將兩模型的運(yùn)動(dòng)方程以分塊矩陣的形式表示：

式中：M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；u、u˙ 和u¨分別為位移、速度和加速度；Peff為施加在人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)的等效節(jié)點(diǎn)力；下標(biāo)i、b 和e 分別代表內(nèi)部域節(jié)點(diǎn)、人工邊界子結(jié)構(gòu)內(nèi)層節(jié)點(diǎn)和外層節(jié)點(diǎn)，上標(biāo)“0”表示自由場(chǎng)。

當(dāng)自由場(chǎng)有限元模型的人工邊界子結(jié)構(gòu)單元網(wǎng)格劃分與實(shí)際有限元模型一致時(shí)：

結(jié)合式(5)～式(7)可得等效節(jié)點(diǎn)力：

根據(jù)式(8)即可得到最終施加于人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)的等效節(jié)點(diǎn)力的表達(dá)式：

具體推導(dǎo)過(guò)程可見(jiàn)文獻(xiàn)[21 - 22]。由式(9)可知，等效節(jié)點(diǎn)力的求解只與人工邊界子結(jié)構(gòu)內(nèi)外層節(jié)點(diǎn)的運(yùn)動(dòng)有關(guān)，而與有限元模型內(nèi)部域節(jié)點(diǎn)無(wú)關(guān)。因此，采用1.1 節(jié)中描述的FK 方法得到層狀半無(wú)限空間中人工邊界子結(jié)構(gòu)位置處的節(jié)點(diǎn)(內(nèi)外層節(jié)點(diǎn))動(dòng)力響應(yīng)后，即可通過(guò)人工邊界子結(jié)構(gòu)法獲得相應(yīng)等效節(jié)點(diǎn)力(外層節(jié)點(diǎn))，進(jìn)而完成內(nèi)部域的FE 分析。

值得注意的是，運(yùn)用式(9)計(jì)算等效節(jié)點(diǎn)力時(shí)需要獲得人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)的位移、速度和加速度，但在實(shí)際應(yīng)用中有限元軟件可自動(dòng)求解速度和加速度，因此只需給定節(jié)點(diǎn)位移即可求得等效節(jié)點(diǎn)力。

1.2.2 三維黏彈性動(dòng)力人工邊界

人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)法向與切向彈簧和阻尼器參數(shù)參考式(10)和式(11)：

式中：CN和CT分別為阻尼器法向和切向阻尼系數(shù)；KN和KT分別為彈簧法向和切向剛度；G為剪切模量；ρ 為質(zhì)量密度；cP為壓縮波速；cS為剪切波速；R為散射波源到人工邊界的垂直距離；αN和αT分別表示法向和切向彈簧的修正系數(shù)，參考文獻(xiàn)[25]分別取為1.34 和0.67。

2 FK-FE 混合方法在ABAQUS 中的二次開(kāi)發(fā)

為實(shí)現(xiàn)FK-FE 混合方法在工程實(shí)際問(wèn)題中的可視化應(yīng)用，通過(guò)ABAQUS 提供的二次開(kāi)發(fā)接口將上述工作開(kāi)發(fā)為具有可視化圖形界面的程序插件。以下對(duì)FK-FE 混合方法在軟件中的主要步驟及插件二次開(kāi)發(fā)流程簡(jiǎn)要說(shuō)明。

2.1 FK-FE 混合方法在ABAQUS 軟件實(shí)現(xiàn)主要步驟

FK-FE 混合方法在ABAQUS 軟件中的實(shí)現(xiàn)主要步驟如圖3 所示，具體操作流程如下：

圖3 FK-FE 混合方法在軟件中的主要實(shí)現(xiàn)步驟Fig.3 Implement steps of the FK-FE hybrid method in software

Step 1：建立土-結(jié)構(gòu)相互作用有限元模型；

Step 2：采用生死單元法，殺死土-結(jié)構(gòu)相互作用有限元模型內(nèi)部單元，建立人工邊界子結(jié)構(gòu)；

Step 3：利用Python 腳本，提取人工邊界子結(jié)構(gòu)所有節(jié)點(diǎn)位置坐標(biāo)；

Step 4：將位置坐標(biāo)傳遞給FK 方法自編Fortran 程序，完成相應(yīng)節(jié)點(diǎn)位移時(shí)程計(jì)算；

Step 5：利用Python 腳本計(jì)算三維黏彈性動(dòng)力人工邊界條件參數(shù)，并賦予人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)；

Step 6：修改.inp 文件，輸入人工邊界子結(jié)構(gòu)兩層節(jié)點(diǎn)位移時(shí)程，計(jì)算得到等效節(jié)點(diǎn)力Peff；

Step 7：修改.inp 文件，施加等效節(jié)點(diǎn)力Peff于人工邊界子結(jié)構(gòu)外層節(jié)點(diǎn)；

Step 8：激活內(nèi)部單元，完成位錯(cuò)點(diǎn)源作用下土-結(jié)構(gòu)相互作用有限元分析，求得建筑結(jié)構(gòu)地震響應(yīng)。

2.2 二次開(kāi)發(fā)流程

本插件開(kāi)發(fā)主要分為FK 方法Fortran 程序在ABAQUS 中的接入和ABAQUS-GUI 程序代碼編寫(xiě) 兩 部 分。 其 中， FK 方 法Fortran 程 序 在ABAQUS 中的接入實(shí)現(xiàn)了層狀半無(wú)限空間地震波場(chǎng)在有限元模型中的輸入，ABAQUS-GUI 程序代碼實(shí)現(xiàn)了FK-FE 混合方法在ABAQUS 中的手動(dòng)操作到內(nèi)核自動(dòng)執(zhí)行程序和可視化界面的轉(zhuǎn)化。

需要注意的是，由于不同版本ABAQUS 所對(duì)應(yīng)的內(nèi)置Python 解釋器不同，進(jìn)行編譯和相應(yīng)numpy 庫(kù)下載時(shí)，需處理好版本兼容問(wèn)題。本文所用ABAQUS 版本為6.14-4，Python 解釋器版本為2.7，NumPy 庫(kù)版本為NumPy 1.11.2。

二次開(kāi)發(fā)設(shè)計(jì)思路如圖4 所示，具體流程如下：

圖4 二次開(kāi)發(fā)設(shè)計(jì)思路Fig.4 Design ideas of secondary development

Step 1：將FK 方法Fortran 子程序?qū)憺?F90 格式，其中.F90 代碼中Module 封裝的變量用Common重新定義為全局變量，合并修改后的.F90 子程序，封裝為子函數(shù)(Subroutine)；

Step 2：利用Python 語(yǔ)言開(kāi)源庫(kù)中的f2py.exe工具(f2py.exe 可通過(guò)下載與Python 對(duì)應(yīng)版本的numpy 庫(kù)得到，是一個(gè)Fortran 程序到Python 程序的接口轉(zhuǎn)換工具)以及GFortran 和mingw32 編譯器對(duì).F90 程序進(jìn)行編譯，形成.pyd 文件(編譯指令為“f2py -m Dislocation_Source -c Dislocation_Source _Abaqus.f90 –fcompiler = gFortran –compiler=mingw32”，封裝的.pyd 模塊名稱(chēng)為Dislocation_Source)；

Step 3：在內(nèi)核執(zhí)行程序中寫(xiě)入.pyd 模塊的調(diào)用命令，關(guān)聯(lián)FK 方法Fortran 程序與二次開(kāi)發(fā)插件(調(diào)用命令為模塊導(dǎo)入指令“import Dislocation_Source”和子函數(shù)調(diào)用指令“Dislocation _Source.Sub”，模塊內(nèi)子函數(shù)名稱(chēng)為Sub)；

Step 4：利用ABAQUS 軟件自帶的RSG 構(gòu)造器進(jìn)行插件程序設(shè)計(jì)，主要內(nèi)容包括設(shè)計(jì)交互界面的“GUI”(圖形界面文件)和負(fù)責(zé)關(guān)聯(lián)腳本函數(shù)的“Kernel”(注冊(cè)文件)兩部分。

2.3 插件簡(jiǎn)介

開(kāi)發(fā)完成的插件由圖形界面文件、注冊(cè)文件、內(nèi)核執(zhí)行程序、.pyd 文件以及.txt 文件組成。各組件功能說(shuō)明見(jiàn)表1，插件界面展示如圖5～圖7 所示。

表1 二次開(kāi)發(fā)過(guò)程各組件功能Table 1 The function of each component in the secondary development process

圖5 建立人工邊界子結(jié)構(gòu)參數(shù)設(shè)置界面Fig.5 Parameter setting interface for building substructure of artificial boundary

圖6 FK 法輸入波場(chǎng)參數(shù)設(shè)置界面Fig.6 Parameter setting interface of calculating the input waves field with the FK method

圖7 施加等效節(jié)點(diǎn)力參數(shù)設(shè)置界面Fig.7 Parameter setting interface of applying equivalent nodal force

開(kāi)發(fā)完成后，將插件放入ABAQUS 的起始工作路徑，即可在 “Plugin-ins”菜單中找到相應(yīng)的插件啟動(dòng)按鈕。插件控件功能說(shuō)明詳見(jiàn)表2。

表2 FK-FE 混合方法控件功能說(shuō)明Table 2 Control function description of FK-FE hybrid method

3 方法驗(yàn)證

本節(jié)分別采用FK 方法和FK-FE 混合方法計(jì)算同一層狀半空間模型在位錯(cuò)點(diǎn)源作用下的地震反應(yīng)來(lái)驗(yàn)證FK-FE 混合方法的精度，采用FK-FE混合方法模擬2021 年漾濞6.4 級(jí)地震動(dòng)，并與強(qiáng)震記錄對(duì)比來(lái)驗(yàn)證方法的有效性。

3.1 FK-FE 混合方法精度驗(yàn)證

驗(yàn)證模型如圖8 所示，場(chǎng)地的速度結(jié)構(gòu)參數(shù)見(jiàn)表3[25]。內(nèi)部域有限元模型尺寸為400 m×400 m×400 m，網(wǎng)格大小10 m。設(shè)置3 個(gè)觀測(cè)臺(tái)站位置坐標(biāo)分別為A(0, 0, 0)、B(0, 100, 200)和C(100, 100, 300)。位錯(cuò)點(diǎn)源埋置在地表下4 km 處，地震矩M0=1022dyne.cm，走向、傾角、滑動(dòng)角分別為135°、70°和30°。震源時(shí)間函數(shù)為Ricker 子波，設(shè)定中心頻率f0為5 Hz，波峰位置為0.3 s，計(jì)算時(shí)間間隔為0.025 s，Ricker 子波時(shí)域和頻域波形如圖9 所示。

表3 驗(yàn)證模型場(chǎng)地速度結(jié)構(gòu)[26]Table 3 Site velocity structure of validate model

圖8 FK-FE 混合方法驗(yàn)證模型Fig.8 Verification analysis model of the FK-FE hybrid method

圖9 震源時(shí)間函數(shù)Fig.9 Source time function

圖10～圖11 分別給出了采用兩種方法得到的觀測(cè)臺(tái)站處的位移、加速度時(shí)程的相對(duì)誤差分布，兩種結(jié)果十分吻合，位移時(shí)程誤差均在±0.5%以?xún)?nèi)，加速度時(shí)程誤差均在±1%以?xún)?nèi)，檢驗(yàn)了本文所建立的FK-FE 混合方法精度。

圖10 FK 方法和FK-FE 混合方法位移時(shí)程對(duì)比Fig.10 Comparison of displacement time history for the FK method and the FK-FE hybrid method

圖11 FK 方法和FK-FE 混合方法加速度時(shí)程對(duì)比Fig.11 Comparison of acceleration time history for the FK method and the FK-FE hybrid method

3.2 FK-FE 混合方法有效性驗(yàn)證

有效性驗(yàn)證模型與精度驗(yàn)證模型類(lèi)似，以觀測(cè)臺(tái)站為原點(diǎn)建立內(nèi)部域有限元模型，模型尺寸為100 m×100 m×100 m，網(wǎng)格大小2 m。參考漾濞6.4 級(jí)地震相關(guān)信息[27]設(shè)定震源參數(shù)，位錯(cuò)點(diǎn)源埋深7.5 km，地震矩M0=4.47×1025dyne.cm，走向、傾角、滑動(dòng)角分別為135°、82°和-165°。震源時(shí)間函數(shù)取鐘形函數(shù)(式12)，τ表示上升時(shí)間，取τ=0.5 s，計(jì)算時(shí)間間隔為0.025 s。漾濞地區(qū)一維地殼波速結(jié)構(gòu)可參考文獻(xiàn)[28]。

圖12 給出了FK-FE 混合方法模擬的2021 年6.4 級(jí)漾濞地震中2 個(gè)遠(yuǎn)場(chǎng)臺(tái)站處地震動(dòng)加速度時(shí)程與實(shí)測(cè)記錄的對(duì)比結(jié)果。53SDX 為施甸臺(tái)站，震中距127.88 km，53YRH 為仁和鎮(zhèn)臺(tái)站，震中距142.33 km。計(jì)算結(jié)果均為50 s 的加速度時(shí)程，曲線(xiàn)末端所標(biāo)數(shù)值為地震時(shí)地面最大峰值加速度PGA，單位為cm/s2，黑色與紅色時(shí)程曲線(xiàn)分別為強(qiáng)震記錄和模擬結(jié)果。比較發(fā)現(xiàn)，本文模擬的地震動(dòng)加速度時(shí)程在幅值、持時(shí)、波形上與強(qiáng)震記錄對(duì)比良好，論證了FK-FE 混合方法的有效性。

圖12 模擬地震動(dòng)時(shí)程與實(shí)測(cè)記錄對(duì)比Fig.12 The time-history comparison between the simulations and the measured records

4 算例分析——以9 層Benchmark鋼框架模型為分析對(duì)象

本節(jié)將建立的FK-FE 混合方法應(yīng)用于美國(guó)SAC Steel Project 的 9 層Benchmark 抗彎鋼框架的全過(guò)程地震反應(yīng)分析，探討了地殼層速度結(jié)構(gòu)與場(chǎng)地類(lèi)別兩個(gè)因素對(duì)該結(jié)構(gòu)地震反應(yīng)的影響。該Benchmark 鋼框架已被國(guó)內(nèi)外多位學(xué)者使用[29-30]，具有一定代表性與普適性，模型具體資料可參考文獻(xiàn)[29]。

4.1 層狀半無(wú)限空間模型

層狀半無(wú)限空間模型場(chǎng)地速度結(jié)構(gòu)參數(shù)見(jiàn)表4[31]。位錯(cuò)點(diǎn)源埋置在地表下18 km 處，全局坐標(biāo)為(20 km, 20 km, 18 km)。震源參數(shù)設(shè)置如下：走向、傾角和滑動(dòng)角分別為135°、70°和30°，地震矩M0=3.55×1026dyn.cm。震源時(shí)間函數(shù)同驗(yàn)證模型，中心頻率f0為5 Hz，截止頻率fmax約為15 Hz，波形如圖9 所示。

表4 層狀半無(wú)限空間模型場(chǎng)地速度結(jié)構(gòu)[31]Table 4 Site velocity structure of layered infinite half-space model

4.2 有限元模型

土結(jié)構(gòu)相互作用會(huì)顯著影響框架結(jié)構(gòu)的地震反應(yīng)[32]。因此，本文運(yùn)用ABAQUS 建立內(nèi)部域土-結(jié)構(gòu)相互作用有限元模型，如圖13 所示。

圖13 土-結(jié)構(gòu)相互作用有限元模型 /mFig.13 Finite element model of soil-structure interaction

模型布置方面，內(nèi)部域土體模型尺寸為200 m×200 m×23.65 m，網(wǎng)格大小為3 m～5 m，滿(mǎn)足15 Hz波場(chǎng)計(jì)算精度；土體模型頂面為自由面，其余邊界面設(shè)置黏彈性邊界；上部Benchmark 鋼框架外層灰色一榀框架采用連接器以MPC-pin 的形式鉸接于主體結(jié)構(gòu)；框架柱底和周?chē)鷪?chǎng)地土進(jìn)行綁定約束。

單元類(lèi)型選取方面，場(chǎng)地土體采用三維實(shí)體單元(C3D8)，Benchmark 鋼框結(jié)構(gòu)采用梁?jiǎn)卧?B31)。

材料選取方面，場(chǎng)地土為彈性，其各項(xiàng)參數(shù)根據(jù)實(shí)際工況中場(chǎng)地速度結(jié)構(gòu)進(jìn)行設(shè)置。鋼材采用雙線(xiàn)性等向強(qiáng)化本構(gòu)模型，塑性階段(屈服后)模量取為1/100 倍的彈性階段彈性模量。有限元模型采用Raiyleigh 阻尼，本模型取場(chǎng)地和結(jié)構(gòu)的第一階和第二階振型自振頻率計(jì)算Raiyleigh 阻尼系數(shù)α 和β，其中土層阻尼比見(jiàn)表4，鋼材阻尼比取0.03。

4.3 結(jié)構(gòu)模型合理性分析

對(duì)9 層Benchmark 鋼框架有限元模型(不含場(chǎng)地)進(jìn)行模態(tài)分析，如表5 所示，將所得前四階振型自振頻率與文獻(xiàn)[29]結(jié)果對(duì)比，自振頻率誤差均控制在±5%以?xún)?nèi)，驗(yàn)證了結(jié)構(gòu)模型的合理性。

表5 鋼框架模型自振頻率與文獻(xiàn)[29]結(jié)果對(duì)比Table 5 Results comparison of this paper and literature [29]for Natural frequency of steel frame model

4.4 地殼層速度結(jié)構(gòu)對(duì)建筑結(jié)構(gòu)地震反應(yīng)的影響

地震波地殼層傳播過(guò)程中受到行波效應(yīng)以及非均勻地形的影響，會(huì)發(fā)生反射和折射，導(dǎo)致幅值和方向發(fā)生變化[33]。為探討地殼速度結(jié)構(gòu)對(duì)建筑結(jié)構(gòu)地震反應(yīng)的影響，以表4 速度結(jié)構(gòu)為基準(zhǔn)(表6 中速度結(jié)構(gòu)1)，以莫霍面壓縮波速(8.10 km/s)為上限，控制地殼層介質(zhì)泊松比與密度不變，遞增地殼各層剪切波速，得到表6 所示三種地殼層介質(zhì)。本文重點(diǎn)關(guān)注了結(jié)構(gòu)一側(cè)一榀框架的動(dòng)力響應(yīng)，各響應(yīng)輸出點(diǎn)位置示意圖14 如所示。

表6 三種地殼層速度結(jié)構(gòu)Table 6 Three crust velocity structures

圖14 結(jié)果展示位置示意圖Fig.14 Schematic plots of position for result display

圖15～圖16 分別給出了FK-FE 自由場(chǎng)模型中地表P點(diǎn)三分量加速度時(shí)程和反應(yīng)譜結(jié)果。為便于觀察，將反應(yīng)譜中所有y向分量縮小5 倍、z向分量縮小25 倍?？梢园l(fā)現(xiàn)，地殼層速度結(jié)構(gòu)軟硬程度的變化對(duì)地震動(dòng)加速度響應(yīng)的影響較為明顯，其中尤以水平向分量最為顯著。隨地殼層介質(zhì)變硬，地震動(dòng)響應(yīng)出現(xiàn)時(shí)刻提前，幅值減小，持時(shí)更短，衰減更加迅速。且反應(yīng)譜以短周期成分為主，幅值隨介質(zhì)變硬而變小，但包含的頻率成分基本不變。例如：相較于速度結(jié)構(gòu)1，速度結(jié)構(gòu)3 工況下地表P點(diǎn)三分量加速度響應(yīng)分別減小約72.6%、86.3%和43.8%；衰減為峰值幅值1/10所用時(shí)間分別減小2.30 s、1.95 s 和1.15 s；反應(yīng)譜幅值減小約52.1%、45.0%和36.6%。

圖15 自由場(chǎng)模型中地表P點(diǎn)加速度時(shí)程Fig.15 Acceleration time history of pointPon the surface for the model of free field

圖16 自由場(chǎng)模型中地表P點(diǎn)加速度三分量加速度反應(yīng)譜Fig.16 Three-component acceleration response-spectrum of pointPon the surface for the model of free field

圖17 給出了FK-FE 土-結(jié)構(gòu)相互作用模型中地表P點(diǎn)三分量加速度時(shí)程。比較土-結(jié)構(gòu)相互作用模型與自由場(chǎng)模型地表加速度時(shí)程發(fā)現(xiàn)，考慮土-結(jié)構(gòu)相互作用的影響后，地表加速度響應(yīng)在三種地殼層速度結(jié)構(gòu)下均有不同程度的放大。例如：速度結(jié)構(gòu)1 中，地表三分量加速度分別放大1.16 倍、1.65 倍和2.19 倍，速度結(jié)構(gòu)2 為1.18 倍、2.23 倍和2.54 倍，速度結(jié)構(gòu)3 則是1.4 倍、2.18 倍和2.61 倍。由此可見(jiàn)，土-結(jié)構(gòu)相互作用對(duì)較硬地殼介質(zhì)下結(jié)構(gòu)地震反應(yīng)的影響更為顯著。

圖17 土-結(jié)構(gòu)相互作用模型中地表P點(diǎn)加速度時(shí)程Fig.17 Acceleration time history of pointPon the surface for the model of soil-structure interaction

圖18 分別給出了不含結(jié)構(gòu)的FK-FE 自由場(chǎng)模型和含有結(jié)構(gòu)的FK-FE 土-結(jié)構(gòu)相互作用模型在不同地殼層速度結(jié)構(gòu)下地表P點(diǎn)的三分量加速度反應(yīng)譜，對(duì)應(yīng)阻尼比為5%。

圖18 兩種模型地表P點(diǎn)三分量加速度反應(yīng)譜Fig.18 Three-component acceleration response-spectrum of pointPon the surface for the two model

在土-結(jié)構(gòu)相互作用的模型中，反應(yīng)譜幅值同樣會(huì)隨地殼介質(zhì)變硬而變小，這一現(xiàn)象在短周期成分較為明顯，而對(duì)長(zhǎng)周期成分影響較小。值得注意的是，考慮土-結(jié)構(gòu)相互作用的影響后，加速度反應(yīng)譜峰值對(duì)應(yīng)周期增大，反應(yīng)譜低頻成分增多，且在土-結(jié)構(gòu)相互作用模型中，峰值對(duì)應(yīng)周期會(huì)隨地殼介質(zhì)的變硬而增大。相較速度結(jié)構(gòu)1，速度結(jié)構(gòu)3 工況下三分量加速度反應(yīng)譜峰值對(duì)應(yīng)周期分別提高了0.1 s、0.075 s 和0.05 s。

圖19 給出了框架角柱柱頂處(圖14)在三種地殼層速度結(jié)構(gòu)條件下x與y方向的位移時(shí)程；圖20給出了結(jié)構(gòu)角柱柱頂產(chǎn)生最大位移時(shí)三種地殼層速度結(jié)構(gòu)下鋼框架變形云圖，其中圖20(a)～圖20(c)為一榀框架(圖14 所示位置框架)x向變形云圖，圖20(d)～圖20(f)為框架整體變形云圖。

圖19 不同地殼層速度結(jié)構(gòu)下頂層角柱位移時(shí)程Fig.19 Displacement time history of top corner pillar in different crust velocity structures

結(jié)果顯示，隨著地殼層剪切波速和壓縮波速的增大，位移峰值點(diǎn)出現(xiàn)時(shí)刻逐漸前移，且峰值大小逐漸減小。例如：三種速度結(jié)構(gòu)x向位移峰值分別出現(xiàn)在11.875 s、10.375 s 和9.25 s，最大相對(duì)位移分別為21.6 cm、19.8 cm 和16.7 cm，和地殼層速度結(jié)構(gòu)1 相比，速度結(jié)構(gòu)3 工況下建筑結(jié)構(gòu)x向最大位移減小了22.7%；y向位移峰值分別出現(xiàn)在11.875 s、10.375 s 和9.275 s，最大相對(duì)位移分別為21.1 cm、19.4 cm 和18.1 cm，和地殼層速度結(jié)構(gòu)1 相比，速度結(jié)構(gòu)3 建筑結(jié)構(gòu)y向最大位移減小了14.3%。

上述現(xiàn)象是由于在較硬的地殼層介質(zhì)中地震波傳播速度較快，P 波初至?xí)r刻提前，且不同地殼層速度結(jié)構(gòu)對(duì)地震波放大作用也不同導(dǎo)致的。

此外，不同地殼層速度結(jié)構(gòu)下，頂層角柱峰值位移時(shí)刻結(jié)構(gòu)最大形變分布沿層高也有所不同。對(duì)于本算例所示的一榀框架，速度結(jié)構(gòu)1 工況下結(jié)構(gòu)最大形變出現(xiàn)在第3 層～4 層的位置；速度結(jié)構(gòu)2 中出現(xiàn)在第4 層～5 層及第8 層～9 層位置；速度結(jié)構(gòu)3 則在第2 層～3 層及第9 層位置。

圖21 給出了三種地殼層速度結(jié)構(gòu)下結(jié)構(gòu)層間位移角的分布情況。由圖可見(jiàn)，隨著地殼層剪切波速和壓縮波速的增大，結(jié)構(gòu)各層層間位移角逐漸減小，相較速度結(jié)構(gòu)1 的位移響應(yīng)，結(jié)構(gòu)在地殼層速度結(jié)構(gòu)3 場(chǎng)地條件下x向?qū)娱g位移角減小為約30.5%～61.3%，y向?qū)娱g位移角減小為約26.9%～62.0%?？梢?jiàn)，地殼層速度結(jié)構(gòu)的變化會(huì)顯著影響結(jié)構(gòu)層間位移角的大小。

圖21 不同地殼層速度結(jié)構(gòu)下最大層間位移角Fig.21 The maximum inter-story drift angle under different cases of crust velocity structures

以上數(shù)據(jù)表明，在抗震設(shè)計(jì)中，應(yīng)考慮所在地區(qū)地殼層速度結(jié)構(gòu)對(duì)建筑地震反應(yīng)的影響，并對(duì)建筑結(jié)構(gòu)進(jìn)行合理規(guī)劃。

4.5 不同場(chǎng)地類(lèi)別對(duì)結(jié)構(gòu)地震反應(yīng)的影響

依據(jù)我國(guó)《建筑抗震設(shè)計(jì)規(guī)范》，進(jìn)行建筑抗震設(shè)計(jì)時(shí)應(yīng)考慮建筑場(chǎng)地對(duì)建筑抗震利害的影響，建筑場(chǎng)地類(lèi)別的劃分是較為重要的參考因素[33]。因此，本文以土層等效剪切波速為設(shè)計(jì)依據(jù)，設(shè)計(jì)了Ⅱ、Ⅲ、Ⅳ三種場(chǎng)地類(lèi)別探究其對(duì)結(jié)構(gòu)地震反應(yīng)的影響。以表4 場(chǎng)地速度結(jié)構(gòu)為基礎(chǔ)，調(diào)整頂部覆蓋層厚度為100 m，覆蓋層參數(shù)見(jiàn)表7。震源參數(shù)設(shè)置同4.1 節(jié)。

表7 場(chǎng)地覆蓋層介質(zhì)物理參數(shù)Table 7 Physical parameter of site overburden layers

圖22～圖23 給出了鋼框架頂層角柱上端在不同場(chǎng)地類(lèi)別下x與y方向的位移時(shí)程結(jié)果以及層間位移角沿高度的分布情況?？梢园l(fā)現(xiàn)，當(dāng)建筑結(jié)構(gòu)位于較硬的場(chǎng)地覆蓋層介質(zhì)上時(shí)，產(chǎn)生的動(dòng)力響應(yīng)較小，隨著場(chǎng)地覆蓋層介質(zhì)剪切波速和壓縮波速的減小，結(jié)構(gòu)的整體位移響應(yīng)及結(jié)構(gòu)各層層間位移角均逐漸增大。例如：Ⅱ類(lèi)場(chǎng)地相比于Ⅳ類(lèi)場(chǎng)地，x向最大位移由24.9 cm 增加至30.1 cm，幅值變化約17.2%；y向最大位移由25.6 cm 增加至32.9 cm，幅值變化約22.3%；x向最大層間位移角約增大16.9%～44.7%，y向最大層間位移角約增大20.8%～30.0%。

圖22 不同場(chǎng)地類(lèi)別下頂層角柱位移時(shí)程Fig.22 Displacement time history of top corner pillar under different cases of different site classes

圖23 不同場(chǎng)地類(lèi)別下最大層間位移角Fig.23 The maximum inter-story drift angle under different cases of different site classes

由此可見(jiàn)，結(jié)構(gòu)在位錯(cuò)點(diǎn)源地震作用下，場(chǎng)地覆蓋層介質(zhì)的軟硬程度會(huì)顯著影響結(jié)構(gòu)的地震反應(yīng)。此外，觀察位移時(shí)程可發(fā)現(xiàn)，場(chǎng)地類(lèi)別對(duì)位移峰值出現(xiàn)時(shí)刻基本無(wú)影響。這是由于，雖然近地表覆蓋層介質(zhì)不同，但覆蓋層厚度較薄(100 m)，因此對(duì)P 波初至?xí)r刻影響不大。

5 結(jié)論

本文基于區(qū)域縮減理論提出了一種FK-FE 混合方法。該方法采用FK 法計(jì)算位錯(cuò)點(diǎn)源作用下整體模型地震波場(chǎng)，運(yùn)用FE 法模擬有限元分析域中土-結(jié)構(gòu)相互作用，實(shí)現(xiàn)了從位錯(cuò)點(diǎn)源到結(jié)構(gòu)的全過(guò)程物理模型分析。以9 層Benchmark 抗彎鋼框架模型為應(yīng)用實(shí)例，研究了地殼層速度結(jié)構(gòu)與場(chǎng)地類(lèi)別對(duì)結(jié)構(gòu)地震反應(yīng)的影響。得到主要結(jié)論如下：

(1) 本文建立的FK-FE 混合方法可模擬震源到建筑結(jié)構(gòu)的全過(guò)程地震反應(yīng)，能夠有效解決寬頻地震動(dòng)作用下地殼層到工程結(jié)構(gòu)的多尺度跨越問(wèn)題，具有較高的求解精度。

(2) 所形成的FK-FE 混合方法ABAQUS 插件具有操作簡(jiǎn)便、應(yīng)用性強(qiáng)的特點(diǎn)，可綜合考慮震源和場(chǎng)地參數(shù)信息。該插件不僅可以對(duì)本文Benchmark抗彎鋼框架進(jìn)行全過(guò)程地震反應(yīng)分析求解，對(duì)其他各種結(jié)構(gòu)同樣適用，例如：高層結(jié)構(gòu)、大跨空間結(jié)構(gòu)、地下結(jié)構(gòu)和橋梁結(jié)構(gòu)等。

(3) 位錯(cuò)點(diǎn)源作用下，地殼層速度結(jié)構(gòu)和場(chǎng)地類(lèi)別對(duì)上部結(jié)構(gòu)的變形和層間位移角有顯著影響。本文算例中，地殼層速度結(jié)構(gòu)3(較硬)相比速度結(jié)構(gòu)1(較軟)，框架最大位移和層間位移角分別減小了22.7%和62.0%；Ⅱ類(lèi)場(chǎng)地相較于Ⅳ類(lèi)場(chǎng)地，框架最大位移和層間位移角分別減小了22.3%和44.7%。