寶 鑫 劉晶波， 李述濤 王 菲

1） 中國北京 100084 清華大學(xué)土木工程系

2） 中國北京 100036 軍事科學(xué)院國防工程研究院

引言

隨著技術(shù)進(jìn)步和經(jīng)濟(jì)發(fā)展以及“海洋強(qiáng)國”戰(zhàn)略的提出，我國對(duì)于海洋資源的開發(fā)和利用已步入快速發(fā)展階段，大量近海交通工程、海上風(fēng)機(jī)和海洋平臺(tái)的規(guī)劃和建設(shè)已提上日程.我國地處環(huán)太平洋地震帶西側(cè)，海域地質(zhì)構(gòu)造復(fù)雜，地震頻發(fā)，海洋工程的地震安全性問題不容忽視（李小軍，2006），海域場(chǎng)地的地震反應(yīng)分析作為海洋工程抗震設(shè)計(jì)和地震安全性評(píng)價(jià)的重要基礎(chǔ)應(yīng)予以高度關(guān)注.

早期對(duì)于海域場(chǎng)地地震動(dòng)的研究主要采用解析方法，將海域地震動(dòng)的產(chǎn)生和傳播簡化為數(shù)學(xué)模型，通過理論推導(dǎo)求解海域場(chǎng)地的應(yīng)力場(chǎng)或波場(chǎng)分布（Lindsay，1939；Brekhovskikh，1980；鄭天愉等，1985；朱鏡清，1988；馮啟民，1990）在揭示流固耦合開放系統(tǒng)的波動(dòng)傳播規(guī)律及其影響因素等方面取得了具有理論指導(dǎo)意義的研究成果.近年來，美國、日本等國家陸續(xù)在其臨近海域布設(shè)地震觀測(cè)系統(tǒng)，利用近海地震臺(tái)站記錄可分析海域地震動(dòng)的時(shí)頻特征和工程特性，Chen 等（2017）、陳蘇等（2018）和楊銘等（2020）分別利用美國加州和日本相模灣等地的海域地震記錄，揭示了近海地震動(dòng)具有長周期成分豐富，豎向分量強(qiáng)度低，反應(yīng)譜峰值大等特點(diǎn).同時(shí)，隨著計(jì)算機(jī)性能水平和數(shù)值模擬方法的成熟與發(fā)展，數(shù)值模擬方法憑借其在求解復(fù)雜模型和邊界條件等方面的優(yōu)勢(shì)，被更多地應(yīng)用于海域場(chǎng)地的地震反應(yīng)分析.采用譜元法（Linket al，2009）、有限差分法（Okamotoet al，2017）、傳遞函數(shù)模型（Liet al，2017）和有限元法（陳少林等，2019）建立了海水-海床流固耦合分析模型與方法，為復(fù)雜海域場(chǎng)地的地震動(dòng)場(chǎng)計(jì)算提供了實(shí)用分析工具.胡進(jìn)軍等（2017）通過數(shù)值模擬，比較了上覆海水層對(duì)海底地震動(dòng)的放大效應(yīng)在豎直方向和水平方向的差異，揭示了海底豎向地震動(dòng)，尤其是P 波引起的豎向震動(dòng)更易受海水層影響.榮棉水等（2013）分析了海底沉積層對(duì)海域地震動(dòng)場(chǎng)和地震動(dòng)參數(shù)的影響表明，當(dāng)考慮海底沉積層時(shí)，海域場(chǎng)地的地震動(dòng)場(chǎng)存在放大現(xiàn)象，且地震動(dòng)反應(yīng)譜周期也有所增大.

然而，目前對(duì)海域場(chǎng)地地震反應(yīng)的理論分析與數(shù)值模擬研究主要針對(duì)平坦的海底場(chǎng)地，而近海工程場(chǎng)地存在向海洋傾斜的單體邊坡，且邊坡一側(cè)直接與海水耦合.此類橫向不均勻的特殊場(chǎng)地類型導(dǎo)致其自由場(chǎng)的求解存在天然的困難，對(duì)現(xiàn)有的基于自由波場(chǎng)的地震波動(dòng)輸入方法提出了挑戰(zhàn)；同時(shí)，場(chǎng)地與無限海水層和半無限海床的動(dòng)力相互作用也是開展近海場(chǎng)地地震反應(yīng)分析時(shí)需重點(diǎn)考慮的問題，其涉及的關(guān)鍵因素包括大范圍流-固耦合系統(tǒng)的數(shù)值模擬、無限固體域和流體域波動(dòng)輻射效應(yīng)的數(shù)值模擬等.

因此，綜合上述因素，本文擬采用基于聲流體單元的流固耦合算法（寶鑫，劉晶波，2017）模擬場(chǎng)地-海水動(dòng)力相互作用，利用流體介質(zhì)動(dòng)力人工邊界（劉晶波等，2017）和一致黏彈性人工邊界單元（劉晶波等，2006）模擬無限流、固介質(zhì)的輻射阻尼，提出一種基于人工邊界子結(jié)構(gòu)，利用混合波場(chǎng)實(shí)現(xiàn)近海場(chǎng)地地震動(dòng)輸入的方法，在此基礎(chǔ)上，建立近海場(chǎng)地地震反應(yīng)分析模型，并對(duì)其適用性及計(jì)算精度進(jìn)行驗(yàn)證.

1 近海場(chǎng)地地震反應(yīng)分析模型

近海場(chǎng)地與下部基巖、外部海水域和海床介質(zhì)共同構(gòu)成半無限的開放系統(tǒng)，在地震反應(yīng)分析中，需綜合考慮局部地形影響和場(chǎng)地與海水的耦合作用，截取有限的近場(chǎng)計(jì)算域進(jìn)行建模分析，并對(duì)流體和固體介質(zhì)截?cái)噙吔邕M(jìn)行處理，以模擬其波動(dòng)輻射效應(yīng).示意模型如圖1 所示，下文將對(duì)該數(shù)值模型所涉及的關(guān)鍵技術(shù)環(huán)節(jié)進(jìn)行簡要介紹.

圖1 近海場(chǎng)地地震反應(yīng)分析模型Fig. 1 Seismic response analysis model of offshore site

1.1 流體介質(zhì)動(dòng)力人工邊界

劉晶波等（2017）基于波動(dòng)理論，將流體介質(zhì)中的單側(cè)波動(dòng)方程轉(zhuǎn)化為截?cái)噙吔缣幍牡刃ЯW(xué)系統(tǒng)，從而提出一種離散化的流體介質(zhì)動(dòng)力人工邊界條件（圖2），可用于吸收流體域截?cái)噙吔缣幍耐庑胁▌?dòng).該流體介質(zhì)人工邊界由阻尼器與集中質(zhì)量構(gòu)成.二維模型中對(duì)應(yīng)的物理參數(shù)質(zhì)量M和阻尼C分別為

圖2 流體介質(zhì)動(dòng)力人工邊界示意圖Fig. 2 Schematic diagram of dynamic artificial boundaries of fluid medium

1.2 一致黏彈性人工邊界單元

式中，G，ρ，cS和cP分別為固體介質(zhì)的剪切模量、密度、橫波和縱波波速；h為人工邊界單元的厚度；R為波源至人工邊界點(diǎn)的距離；αT和αN為人工邊界參數(shù)，劉晶波等（2006）給出的推薦值分別為αT＝0.5 和αN＝1，α＝αN/αT.

1.3 近海場(chǎng)地-海水動(dòng)力相互作用有限元模型

基于通用有限元計(jì)算平臺(tái)Ansys 建立近海場(chǎng)地-海水動(dòng)力相互作用有限元模型.其中，固體介質(zhì)和流體介質(zhì)分別采用Solid45 和Fluid80 單元建模，并通過耦合流-固界面節(jié)點(diǎn)的法向自由度，模擬海水與場(chǎng)地的動(dòng)力相互作用（寶鑫，劉晶波，2017）.由于Solid45 和Fluid80 均為三維單元，應(yīng)用于二維模型時(shí)，可僅在平面內(nèi)建立一層單元，并約束全部節(jié)點(diǎn)的出平面自由度.分別在流、固介質(zhì)截?cái)噙吔缣幪砑忧皟晒?jié)中介紹的流體介質(zhì)動(dòng)力人工邊界和固體介質(zhì)一致黏彈性人工邊界單元.此外，為保證波動(dòng)問題的數(shù)值模擬精度，離散化網(wǎng)格的尺寸 Δx應(yīng)滿足（杜修力，2009）

式中，λmin為離散網(wǎng)格模型中波動(dòng)傳播的最短波長，cmin為介質(zhì)中的最小波速，fmax為波動(dòng)問題數(shù)值模擬的截止頻率.

2 基于混合波場(chǎng)的近海場(chǎng)地地震動(dòng)輸入方法

在場(chǎng)地地震反應(yīng)分析中，近場(chǎng)計(jì)算域的截取和人工邊界條件的施加將直接影響地震動(dòng)的輸入過程.以域縮減法（Bielaket al，2003；Yoshimuraet al，2003）和波動(dòng)法（劉晶波，呂彥東，1998）為代表，目前常用的地震動(dòng)輸入方法一般以自由場(chǎng)地震動(dòng)為輸入地震動(dòng)場(chǎng)，通過將其轉(zhuǎn)化為等效地震荷載，在不影響人工邊界對(duì)于外行波動(dòng)吸收的前提下，將地震動(dòng)輸入至計(jì)算模型內(nèi)部.在此基礎(chǔ)上，為進(jìn)一步簡化計(jì)算流程，劉晶波等（2018）根據(jù)波動(dòng)法原理和有限元理論，提出一種土-結(jié)構(gòu)相互作用系統(tǒng)地震波動(dòng)輸入的人工邊界子結(jié)構(gòu)法.該方法不依賴于人工邊界條件的具體形式，將等效地震荷載的計(jì)算轉(zhuǎn)化為人工邊界子結(jié)構(gòu)模型的動(dòng)力分析，即保證了計(jì)算精度又有效地提高了分析效率.但上述地震動(dòng)輸入方法為了便于通過Shake91 和EERV 等一維化場(chǎng)地分析軟件計(jì)算得到自由場(chǎng)分布，主要分析橫向均勻、縱向成層的規(guī)則場(chǎng)地形式，而對(duì)于以圖1 中近海場(chǎng)地為代表的不規(guī)則場(chǎng)地，由于模型不同位置處的地形特征和地層分布存在明顯的差異，地震波入射下整體模型的自由波場(chǎng)難以預(yù)先求解，使得傳統(tǒng)的基于自由場(chǎng)地地震動(dòng)的地震波動(dòng)輸入方法難以實(shí)現(xiàn).為解決這一問題，本文提出一種基于人工邊界子結(jié)構(gòu)，利用混合波場(chǎng)實(shí)現(xiàn)近海場(chǎng)地中地震P 波和SV 波垂直輸入的方法，實(shí)施思路如下.

將近海場(chǎng)地模型中的輸入波場(chǎng)分成左側(cè)、右側(cè)和底部輸入波場(chǎng)三部分，如圖3所示.對(duì)于模型兩側(cè)的波動(dòng)輸入，可根據(jù)截?cái)辔恢锰幍牡貙犹匦越⒌刃У目v向成層、橫向均勻的場(chǎng)地模型，并根據(jù)輸入地震動(dòng)計(jì)算相應(yīng)的自由場(chǎng)分布，進(jìn)而利用人工邊界子結(jié)構(gòu)法，將該自由場(chǎng)分布轉(zhuǎn)化為等效地震荷載；而對(duì)于近場(chǎng)模型的底部截?cái)噙吔?，則僅考慮入射波場(chǎng)，并利用底部的人工邊界子結(jié)構(gòu)將其轉(zhuǎn)化為等效地震荷載.將計(jì)算得到的不同截?cái)噙吔缣幍牡刃У卣鸷奢d施加于近海場(chǎng)地模型，可完成地震動(dòng)的輸入過程.

圖3 基于混合波場(chǎng)的近海場(chǎng)地地震動(dòng)輸入方法Fig. 3 Seismic wave input method of offshore site based on hybrid wave field

從波動(dòng)傳播的角度分析該方法的實(shí)現(xiàn)過程：入射波場(chǎng)分別從側(cè)邊和底面邊界輸入至計(jì)算模型，與內(nèi)部局部地形相互作用產(chǎn)生反射波和散射波.對(duì)于側(cè)面邊界，由于輸入的自由波場(chǎng)同時(shí)包含了入射波和反射波成分，因此經(jīng)模型自由表面反射的外行波動(dòng)被與其對(duì)應(yīng)的等效地震荷載所抵消，而由局部地形產(chǎn)生的散射波則被人工邊界吸收；對(duì)于底面邊界，由于在等效地震荷載的計(jì)算中僅考慮了入射波場(chǎng)，因此反射波場(chǎng)和散射波場(chǎng)都將被人工邊界所吸收.此外，在兩種波場(chǎng)輸入位置的交界處，即模型兩側(cè)底部角點(diǎn)位置，可能存在輸入波場(chǎng)不協(xié)調(diào)的問題，由此引起的額外波動(dòng)最終也將被人工邊界所吸收.

根據(jù)以上思路提出的基于混合波場(chǎng)的近海場(chǎng)地地震反應(yīng)分析方法避免了預(yù)先求解不規(guī)則場(chǎng)地自由場(chǎng)分布的困難，通過構(gòu)建由模型兩側(cè)的自由波場(chǎng)和模型底部的入射波場(chǎng)組成的混合波場(chǎng)，并利用人工邊界子結(jié)構(gòu)法將其輸入至近場(chǎng)計(jì)算模型，實(shí)現(xiàn)此類不規(guī)則場(chǎng)地的地震反應(yīng)分析.該方法的具體實(shí)施步驟如下：

1） 建立近海場(chǎng)地有限元模型，在此基礎(chǔ)上截取對(duì)應(yīng)的人工邊界子結(jié)構(gòu)模型，該模型僅由人工邊界及與其相鄰的一層內(nèi)部介質(zhì)單元構(gòu)成，如圖4a 所示.

圖4 基于混合波場(chǎng)的近海場(chǎng)地地震波動(dòng)輸入方法實(shí)現(xiàn)步驟（a） 利用子結(jié)構(gòu)模型求解等效地震荷載；（b） 在近海場(chǎng)地模型中施加等效地震荷載Fig. 4 Implementation steps of seismic wave input method of offshore site based on hybrid wave field（a） Solving equivalent seismic loads through the substructure model；（b） Applying equivalent seismic loads on the offshore site model

2） 根據(jù)左右兩側(cè)邊界位置處的地層特性分別進(jìn)行等效自由場(chǎng)分析，獲得自由波場(chǎng)，其中右側(cè)截?cái)噙吔鐚?duì)應(yīng)的自由場(chǎng)模型為上覆海水層的半空間場(chǎng)地模型，Bao 等（2020）給出了不同類型波動(dòng)輸入下此類場(chǎng)地自由場(chǎng)分布的理論公式；固定人工邊界子結(jié)構(gòu)模型的最外層節(jié)點(diǎn)，對(duì)所有與側(cè)邊人工邊界相連的子結(jié)構(gòu)節(jié)點(diǎn)施加自由波場(chǎng)位移時(shí)程，而在其余的與底部人工邊界相連的子結(jié)構(gòu)節(jié)點(diǎn)上輸入入射波場(chǎng)位移時(shí)程（圖4a）；對(duì)子結(jié)構(gòu)模型進(jìn)行動(dòng)力分析，求得全部人工邊界節(jié)點(diǎn)上的反力，即為等效地震荷載.

3） 對(duì)近海場(chǎng)地有限元模型的人工邊界節(jié)點(diǎn)施加由上一步驟獲得的等效地震荷載（圖4b），進(jìn)行動(dòng)力時(shí)程計(jì)算，即可完成近海場(chǎng)地的地震反應(yīng)分析.

3 算例與分析

建立近海場(chǎng)地有限元模型，對(duì)本文方法的有效性和計(jì)算精度進(jìn)行驗(yàn)證.假定近海場(chǎng)地、基巖和海床為均一介質(zhì)，材料參數(shù)列于表1.近場(chǎng)模型的整體尺寸為100 m×50 m，水深10 m，向海坡傾角為45° （圖5）.在截?cái)噙吔缣幏謩e添加流、固介質(zhì)人工邊界.采用本文方法將圖6所示的持時(shí)為0.2 s 的脈沖波分別以SV 波和P 波的形式垂直輸入至計(jì)算模型.模型中有限元網(wǎng)格的最大尺寸為1 m，滿足式（3）給出的計(jì)算精度要求.

圖5 計(jì)算模型示意圖Fig. 5 Diagram of calculation model

圖6 脈沖波時(shí)程Fig. 6 Time history of pulse wave

表1 介質(zhì)材料參數(shù)Table 1 Material parameters of media

由于本文方法誤差的主要來源可能是底部人工邊界難以充分吸收下行反射波，導(dǎo)致該波動(dòng)經(jīng)底部邊界再次反射后傳播至模型頂部的近海場(chǎng)地區(qū)域，對(duì)計(jì)算精度產(chǎn)生影響.可將模型底部邊界取至足夠遠(yuǎn)，使反射波在計(jì)算時(shí)長內(nèi)不傳回頂部場(chǎng)地觀測(cè)點(diǎn)，從而獲得擴(kuò)展網(wǎng)格解，作為本文方法的驗(yàn)證.

3.1 SV 波垂直入射

首先分析SV 波入射下近海場(chǎng)地地震動(dòng)場(chǎng)分布情況.以模型水平向總長度L＝100 m 對(duì)近海場(chǎng)地、向海坡以及海床位置的測(cè)點(diǎn)坐標(biāo)進(jìn)行歸一化處理，以該無量綱坐標(biāo)為縱坐標(biāo)，對(duì)比采用本文方法和通過擴(kuò)展網(wǎng)格方法計(jì)算得到的位移波形，如圖7 所示.同時(shí)，為衡量近海場(chǎng)地的地震動(dòng)放大效應(yīng)，定義地震反應(yīng)峰值比R為

圖7 SV 波入射下近海場(chǎng)地的水平（a）和垂直（b）位移波形Fig. 7 Displacement waveforms in horizontal （a） and vertical （b） directions on offshore site under incident SV wave

其中u為位移，下標(biāo)x和z分別代表水平和豎直方向，下標(biāo)0 代表輸入同一地震動(dòng)下均勻半空間平坦場(chǎng)地的地震反應(yīng)，其峰值為入射波峰值的二倍，x為測(cè)點(diǎn)水平坐標(biāo)，t為時(shí)間.

分析SV 波入射下近海場(chǎng)地位移峰值放大系數(shù)的空間分布（圖8）可以看出，本文方法計(jì)算得到的近海場(chǎng)地峰值地震反應(yīng)與擴(kuò)展網(wǎng)格解吻合良好.在靠近向海坡的近海場(chǎng)地區(qū)域，水平地震峰值放大系數(shù)大于1；向海坡由上至下，峰值比逐漸減小，至向海坡底部角點(diǎn)處峰值比達(dá)到最小值；在遠(yuǎn)離向海坡的區(qū)域，近海場(chǎng)地與海床處的峰值放大系數(shù)均接近于1，說明在此位置處的場(chǎng)地地震反應(yīng)受向海坡局部地形的影響較小.進(jìn)一步分析圖7 可知，由于采用混合波場(chǎng)進(jìn)行波動(dòng)輸入時(shí)，模型底面的人工邊界難以充分吸收下行的反射波，導(dǎo)致初至波峰過后的后續(xù)波動(dòng)存在較小的計(jì)算誤差.但鑒于峰值地震反應(yīng)是場(chǎng)地地震反應(yīng)分析和結(jié)構(gòu)抗震設(shè)計(jì)中最為重要的的地震動(dòng)參數(shù)，本文方法對(duì)于此類不規(guī)則的近海場(chǎng)地地震反應(yīng)分析具有較為良好的適用性.

圖8 SV 波入射下近海場(chǎng)地位移峰值放大系數(shù)的空間分布Fig. 8 Spatial distributions of peak displacement amplification coefficient on offshore site under incident SV wave

3.2 P 波垂直入射

進(jìn)一步將圖6 所示的脈沖波以P 波的形式垂直輸入至計(jì)算模型中，比較近海場(chǎng)地、斜坡以及海床位置處的位移波形（圖9），同時(shí)計(jì)算P 波入射下近海場(chǎng)地位移峰值放大系數(shù)的空間分布（圖10）.

圖9 P 波入射下近海場(chǎng)地的水平（a）和垂直（b）位移波形Fig. 9 Displacement waveforms in horizontal （a） and vertical （b） directions on offshore site under incident P wave

圖10 P 波入射下近海場(chǎng)地位移峰值放大系數(shù)的空間分布Fig. 10 Spatial distributions of peak displacement amplification coefficient on offshore site under incident P wave

與SV 波入射時(shí)的計(jì)算結(jié)果（圖8）類似，P 波垂直入射下，采用本文方法與擴(kuò)展網(wǎng)格方法計(jì)算得到的峰值地震反應(yīng)也具有良好的一致性，且地震動(dòng)場(chǎng)的分布規(guī)律與SV 波入射時(shí)的結(jié)果相似，均體現(xiàn)為靠近向海坡的近海場(chǎng)地區(qū)域地震動(dòng)場(chǎng)被放大，而向海坡由上至下峰值地震反應(yīng)逐漸減小.觀察圖9 可以看出，采用本文方法計(jì)算得到的初至波峰過后的后續(xù)波動(dòng)（0.6—1.0 s）與擴(kuò)展網(wǎng)格解相比也存在一定的誤差，這同樣是由于底部邊界對(duì)下行反射波吸收不充分導(dǎo)致的.整體而言，采用本文方法能夠較為合理準(zhǔn)確的實(shí)現(xiàn)P 波入射下近海場(chǎng)地的地震動(dòng)輸入和地震反應(yīng)分析.

4 討論與結(jié)論

針對(duì)橫向不規(guī)則的近海場(chǎng)地地震反應(yīng)問題，本文采用由場(chǎng)地近場(chǎng)截?cái)嗄Ｐ蛢蓚?cè)的自由波場(chǎng)和模型底部的入射波場(chǎng)構(gòu)成的混合波場(chǎng)作為輸入波場(chǎng)，改進(jìn)基于人工邊界子結(jié)構(gòu)的地震波從輸入方法以實(shí)現(xiàn)近海場(chǎng)地的地震動(dòng)輸入.同時(shí)，綜合考慮場(chǎng)地-海水動(dòng)力相互作用以及遠(yuǎn)場(chǎng)流、固介質(zhì)的輻射阻尼，構(gòu)建了一種近海場(chǎng)地地震反應(yīng)分析模型與時(shí)域整體分析方法.

數(shù)值算例的結(jié)果表明，采用本文模型與方法能較為合理可靠地計(jì)算得到的近海場(chǎng)地的峰值地震反應(yīng).鑒于峰值地震反應(yīng)是場(chǎng)地地震反應(yīng)分析和結(jié)構(gòu)抗震設(shè)計(jì)中最為研究人員所關(guān)心的地震動(dòng)參數(shù)，本文方法對(duì)于此類不規(guī)則的近海場(chǎng)地地震反應(yīng)分析具有較為良好的適用性.

采用本文提出的混合波場(chǎng)方法進(jìn)行波動(dòng)輸入時(shí)，由于模型底部的人工邊界難以充分吸收下行反射波，導(dǎo)致計(jì)算得到的場(chǎng)地地震動(dòng)的后續(xù)波動(dòng)存在一定的誤差.進(jìn)一步合理優(yōu)化近海場(chǎng)地模型的底部輸入波場(chǎng)，是該領(lǐng)域的后續(xù)研究中有待解決的關(guān)鍵問題.