李會(huì)芳，趙 密，杜修力

(北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100124)

在地下結(jié)構(gòu)抗震分析中，有限元法因其靈活高效的特點(diǎn)，被廣泛用于有限域模擬。為了更有效地反映截?cái)酂o限土體的輻射阻尼效應(yīng)，需要引入人工邊界條件[1]。人工邊界的適用性和精度直接影響計(jì)算結(jié)果的精度，并通過影響計(jì)算區(qū)域的大小而影響計(jì)算效率。城市抗震防災(zāi)和重大工程項(xiàng)目建設(shè)中經(jīng)常遇到豎向成層介質(zhì)以及不規(guī)則地形場(chǎng)地。在我國西南部山嶺地區(qū)中這種地形地質(zhì)條件較為常見，如滇中引水工程中香爐山隧道穿越的含多條斷裂帶的復(fù)雜地形地質(zhì)條件[2]。廈門市軌道交通2 號(hào)線跨海盾構(gòu)隧道穿越了基巖、基巖強(qiáng)風(fēng)化、基巖全風(fēng)化等多個(gè)地層[3]。不均勻介質(zhì)和不規(guī)則地形場(chǎng)地條件引起復(fù)雜的波動(dòng)反射和透射，會(huì)造成與均勻半空間場(chǎng)地明顯不同的動(dòng)力反應(yīng)。復(fù)雜場(chǎng)地中結(jié)構(gòu)動(dòng)力反應(yīng)研究表明，不均勻地質(zhì)和不規(guī)則地形的存在會(huì)對(duì)場(chǎng)地和結(jié)構(gòu)的動(dòng)力反應(yīng)產(chǎn)生顯著影響[4-8]，在地下結(jié)構(gòu)抗震分析領(lǐng)域引起越來越多的關(guān)注。建立豎向成層介質(zhì)中的高精度人工邊界條件具有研究和實(shí)際意義。國內(nèi)外研究者已提出大量人工邊界條件[9]，其中包括：粘性邊界[10]、粘彈性邊界[11]、透射邊界[12]、無限元法[13-14]、邊界元法[15-16]、Engquist-Majda 邊界[17]、Higdon邊界[18]、Bayliss-Turkel 邊界[19]、Dirichlet-to-Neumann(DtN) 邊界[20-23]、一致邊界(薄層法)[24-26]、比例邊界有限元法[27]、吸收層法如完美匹配層法[28-29]等。以上人工邊界條件主要用于均勻介質(zhì)的半空間或全空間場(chǎng)地以及底部固定的單層或水平成層場(chǎng)地。

近年來適用于水平成層半空間場(chǎng)地的人工邊界條件得到進(jìn)一步發(fā)展。2009 年，蔣通和田治見宏[30]提出一致邊界結(jié)合粘性邊界來處理水平成層半空間介質(zhì)的波傳播問題，粘性邊界中阻尼器的吸能效果對(duì)整體精度影響較大。2011 年，Lee 和Tassoulas[31]提出一致邊界結(jié)合連分式吸收層的方法；2012 年，Jo?o 等[32]提出采用一致邊界結(jié)合完美匹配離散層來解決水平成層半空間介質(zhì)的波傳播問題，分析表明連分式吸收層可等效轉(zhuǎn)化為完美匹配離散層，其中匹配層厚度為基于頻率和外行波入射角的純虛數(shù)。2015 年，Hamdan 等[33]提出采用一致邊界結(jié)合旁軸邊界來模擬水平成層半空間介質(zhì)中波傳播問題，旁軸邊界采用二階泰勒展開代替半空間的解析剛度精度較低。2020 年，李會(huì)芳等[34]提出H 形高精度人工邊界處理水平成層半空間標(biāo)量波傳播問題，此邊界應(yīng)用簡單方便且能得到高精度、高效率的求解。

基于水平成層半空間場(chǎng)地中人工邊界條件的研究工作，本文針對(duì)含豎向成層介質(zhì)以及不規(guī)則地形場(chǎng)地中標(biāo)量波傳播問題，建立了頻域下折線形高精度人工邊界條件。豎向成層介質(zhì)的底邊界采用廣義一致邊界來擬合，通過適當(dāng)?shù)目臻g變換并精確模擬兩側(cè)的半空間剛度，將原來適用于表面自由、底部固定的一致邊界擴(kuò)展應(yīng)用于兩側(cè)開放的豎向成層介質(zhì)底邊界；成層介質(zhì)側(cè)邊界采用基于連分式的高精度邊界。由于提出邊界的高精度特性，有限域可以取得盡量小，對(duì)于無內(nèi)域結(jié)構(gòu)僅關(guān)心地表反應(yīng)的情況，可以直接將人工邊界加在地表，極大地減少自由度，提高計(jì)算效率。通過在每個(gè)豎層邊界上引入斜角變換，廣義一致邊界可以為任意折線形，適用于豎向成層、不規(guī)則地表的場(chǎng)地反應(yīng)。

1 問題描述

地震作用下重大基礎(chǔ)設(shè)施結(jié)構(gòu)的動(dòng)力響應(yīng)需考慮與周圍介質(zhì)的相互作用，形成無限域介質(zhì)中的波傳播問題。圖1(a)所示為二維豎向成層半空間場(chǎng)地中的標(biāo)量波傳播問題，場(chǎng)地中包含任意不規(guī)則地形以及廣義結(jié)構(gòu)。為高效求解出平面波動(dòng)問題，引入人工邊界分別將豎向地層和兩側(cè)無限域截?cái)?，如圖1(b)所示計(jì)算模型1 由人工邊界條件與內(nèi)域有限元部分組成，其中有限域中可以包含任意不均勻、非線性材料以及不規(guī)則結(jié)構(gòu)，采用有限元法模擬。由于提出人工邊界條件的高精度，有限域尺寸可以盡量小。對(duì)于不含地下結(jié)構(gòu)僅關(guān)心地面反應(yīng)的情況如圖1(c)所示，可以直接將人工邊界加在地表，僅用折線形高精度人工邊界條件來模擬豎向成層半空間場(chǎng)地，如圖1(d)中計(jì)算模型2，最大限度減少計(jì)算自由度高效求解地表反應(yīng)。

圖1 豎向成層且地表起伏的半空間場(chǎng)地中標(biāo)量波傳播問題Fig. 1 Scalar wave propagation in half space including vertical stratified media and irregular topography

半空間場(chǎng)地中，左右兩側(cè)無限域以及豎向無限延伸的地層內(nèi)為均勻線彈性材料。頻域下標(biāo)量波傳播問題的控制方程為：

地表為自由邊界條件，各豎向地層和無限域之間的交界面以及人工邊界兩側(cè)均滿足位移和應(yīng)力連續(xù)條件。另外，滿足無限遠(yuǎn)處輻射條件以及初始靜止條件。

2 高精度人工邊界條件

以圖1(b)中計(jì)算模型1 為例給出本文高精度折線形人工邊界條件的實(shí)現(xiàn)過程。不含廣義結(jié)構(gòu)時(shí)，人工邊界條件可加在地表處，此時(shí)計(jì)算模型中無內(nèi)域和豎層側(cè)邊界，計(jì)算模型1 簡化為圖1(d)中計(jì)算模型2。另外，計(jì)算模型1 中豎層底邊界可以根據(jù)結(jié)構(gòu)形式、計(jì)算效率等需要設(shè)置為任意角度折線形邊界。

2.1 兩側(cè)無限域連分式邊界

則側(cè)邊無限域y方向的動(dòng)力剛度為：

以上為左側(cè)無限域的人工邊界條件，同理可以得到右側(cè)無限域的人工邊界條件。

2.2 豎層底部廣義一致邊界

2.3 有限域與人工邊界條件耦合

將有限域側(cè)邊力位移關(guān)系的連分式展開式(11)沿豎向x軸離散，通過線性插值得到側(cè)邊人工邊界條件有限元方程：

將有限域與底部人工邊界條件式(17)以及兩側(cè)人工邊界條件式(18)組裝，得到頻域下耦合系統(tǒng)有限元方程：

式中，下標(biāo)I、B、C分別對(duì)應(yīng)內(nèi)域自由度、邊界自由度以及連分式引入的輔助自由度。

3 數(shù)值算例

本節(jié)給出不同地形、地質(zhì)條件下四種場(chǎng)地中波傳播問題的數(shù)值算例來驗(yàn)證提出高精度折線形人工邊界條件的有效性和精度。參考解為大區(qū)域有限元模型的時(shí)域解，模型尺寸足夠大能夠保證在觀測(cè)時(shí)間內(nèi)結(jié)果不受到截?cái)噙吔缣幏瓷洳ǖ挠绊?。同時(shí)采用了簡單常用的粘彈性人工邊界條件[11]進(jìn)行了對(duì)比分析，說明提出的折線形人工邊界條件適用于豎向成層不規(guī)則地形場(chǎng)地，且在計(jì)算精度上有顯著改進(jìn)。本文人工邊界條件中連分式階數(shù)統(tǒng)一取為5 階，階數(shù)為零時(shí)可退化為旁軸近似。采用Ricker 波作為荷載，其時(shí)程和頻譜曲線如圖2 所示。

圖2 Ricker 波的時(shí)程和頻譜曲線Fig. 2 Time history and Fourier spectrum of Ricker wavelet impulse

3.1 水平地表場(chǎng)地

平直地表半空間場(chǎng)地中含有一個(gè)豎向無限延伸地層，如圖3 所示。中間地層的剪切波速c1=100 m/s，兩側(cè)無限域介質(zhì)剪切波速c2= 200 m/s。坐標(biāo)原點(diǎn)取地表中心點(diǎn)A。計(jì)算中選取了兩種形式的人工邊界，如圖3 中點(diǎn)劃線和虛線所示。加載點(diǎn)為點(diǎn)A，觀察點(diǎn)為A(0, 0)、B(0, -100)、C(100, 0)、D(100, -100)。有限元網(wǎng)格尺寸為2 m，計(jì)算時(shí)間步長為0.01 s。

圖3 含豎直地層半空間場(chǎng)地中標(biāo)量波傳播問題Fig. 3 Scalar wave propagation in half space field with one vertical layer

觀察點(diǎn)的位移時(shí)程結(jié)果如圖4 所示。將采用新邊界的兩種計(jì)算模型和采用粘彈性邊界的計(jì)算模型1 與采用大區(qū)域有限元模型得到的參考解進(jìn)行比較。其中計(jì)算模型1 由人工邊界條件與內(nèi)域有限元部分組成，如圖3 中點(diǎn)劃線所示區(qū)域；計(jì)算模型2 為直接將人工邊界加在地表，如圖3 中虛線所示。下文中圖例意義相同。從圖中可以看出，無論是加載點(diǎn)還是邊界處節(jié)點(diǎn)以及材料交界面上節(jié)點(diǎn)，采用新人工邊界條件的兩種模型計(jì)算結(jié)果均與參考解吻合較好，而粘彈性邊界模型在邊界及材料交界面節(jié)點(diǎn)上計(jì)算誤差較大。驗(yàn)證了新提出人工邊界條件的有效性和精度。

圖4 水平地表場(chǎng)地中觀察點(diǎn)位移時(shí)程結(jié)果Fig. 4 Time histories of displacement solutions at observation points in a flat surface site

3.2 折線地表場(chǎng)地

折線形地表半空間場(chǎng)地中含有兩個(gè)豎向地層，如圖5 所示。左側(cè)地層的剪切波速c1=100 m/s，右側(cè)地層的剪切波速c2= 300 m/s，兩側(cè)無限域介質(zhì)剪切波速c3= 200 m/s。圖5 中各點(diǎn)坐標(biāo)為A(0, 0)、B(50, 50)、C(150, 0)、D(0, -50)、E(50, -50)、F(50,-50)，其中加載點(diǎn)為B，觀察點(diǎn)為B、C、E、F。有限元網(wǎng)格尺寸為1 m～3 m，計(jì)算時(shí)間步長為0.002 s。

圖5 含豎向成層介質(zhì)且地表起伏的半空間場(chǎng)地中標(biāo)量波傳播問題Fig. 5 Scalar wave propagation in half space field with vertical stratified media and irregular topography

觀察點(diǎn)的位移時(shí)程結(jié)果如圖6 所示。觀察點(diǎn)包括地表節(jié)點(diǎn)、人工邊界上的節(jié)點(diǎn)以及材料交界面處的節(jié)點(diǎn)。從圖6 中可以看出，兩種新邊界計(jì)算模型的位移解與參考解吻合較好，而粘彈性邊界模型在邊界及材料交界面節(jié)點(diǎn)上計(jì)算誤差較大。表明提出的人工邊界條件可以很好的模擬邊界處和材料變化帶來的反射和散射波，以及隨時(shí)間空間變化的輻射效應(yīng)，驗(yàn)證了新人工邊界條件的有效性和精度。

圖6 折線地表場(chǎng)地中觀察點(diǎn)位移時(shí)程結(jié)果Fig. 6 Time histories of displacement solutions at observation points in a zigzag terrain site

3.3 階梯地形場(chǎng)地

階梯地形半空間場(chǎng)地中含有三個(gè)豎向地層，如圖7 所示。從左到右介質(zhì)剛度依次增大，左側(cè)無限域介質(zhì)剪切波速cl= 100 m/s，三個(gè)豎向地層的剪切波速分別為c2= 150 m/s、c3= 200 m/s、c4=250 m/s，右側(cè)無限域介質(zhì)剪切波速c5= 300 m/s。圖7 中定位點(diǎn)坐標(biāo)為A(0, 0)、B(50, 50/3)、C(100,100/3)、D(150, 50)、E(0, -50)、F(50, -50)、G(100,-50)、H(150,-50)，其中加載點(diǎn)為B，觀察點(diǎn)為A、C、E、G。有限元網(wǎng)格尺寸為1 m～2.2 m，計(jì)算時(shí)間步長為0.002 s。

圖7 含豎向成層介質(zhì)的階梯地形半空間場(chǎng)地中標(biāo)量波傳播問題Fig. 7 Scalar wave propagation in step-shaped half space field with vertical stratified media

觀察點(diǎn)的位移時(shí)程結(jié)果如圖8 所示。觀察點(diǎn)包括地表節(jié)點(diǎn)、人工邊界上的節(jié)點(diǎn)以及材料交界面處的節(jié)點(diǎn)。從圖8 中可以看出，兩種新邊界計(jì)算模型的位移解與參考解吻合較好，而粘彈性邊界模型在邊界及材料交界面節(jié)點(diǎn)上計(jì)算誤差較大。表明新人工邊界條件可以很好的吸收邊界以及材料變化產(chǎn)生的反射和散射波，模擬無限場(chǎng)地的輻射效應(yīng)，驗(yàn)證了新人工邊界條件的有效性和精度。

圖8 階梯地形場(chǎng)地中觀察點(diǎn)位移時(shí)程結(jié)果Fig. 8 Time histories of displacement solutions at observation points in a stepped terrain site

3.4 含地下結(jié)構(gòu)場(chǎng)地

如圖9 所示，地下結(jié)構(gòu)建立在含豎向成層介質(zhì)且兩側(cè)不等高地形的半空間場(chǎng)地中，其地表高低起伏可模擬河谷與山坡同時(shí)存在的復(fù)雜地形，本算例分析其中的標(biāo)量波傳播問題。從左到右5 種介質(zhì)剛度依次增大，左側(cè)無限域介質(zhì)剪切波速cl= 100 m/s，三個(gè)豎向?qū)拥募羟胁ㄋ俜謩e為c2=150 m/s，c3= 200 m/s，c4= 250 m/s，右側(cè)無限域介質(zhì)剪切波速c5= 300 m/s。因含有地下結(jié)構(gòu)，計(jì)算模型2 不再適用或者說誤差較大，本節(jié)僅采用計(jì)算模型1 與參考解比較，人工邊界位置如圖9中點(diǎn)劃線所示，分別施加新邊界條件和粘彈性邊界條件[11]。圖9 中定位點(diǎn)坐標(biāo)為A(0, 0)、B(50, -20)、C(150, 30)、D(200, 10)、E(0, -50)、F(50, -50)、G(150, -50)、H(200, -50)，地下結(jié)構(gòu)圓心坐標(biāo)為(100, -30)，半徑10 m。其中加載點(diǎn)為B，觀察點(diǎn)為B、C、E、H以及地下結(jié)構(gòu)上節(jié)點(diǎn)C1(90, -30)、C2(100, -20)。有限元網(wǎng)格最大尺寸小于4 m，計(jì)算時(shí)間步長為0.002 s。

圖9 含地下結(jié)構(gòu)的復(fù)雜地形地質(zhì)半空間場(chǎng)地中標(biāo)量波傳播問題Fig. 9 Scalar wave propagation in complex topographical geological half space field including an underground structure

觀察點(diǎn)的位移時(shí)程結(jié)果如圖10 所示。觀察點(diǎn)包括地表節(jié)點(diǎn)、人工邊界和材料界面處的節(jié)點(diǎn)以及地下結(jié)構(gòu)上的節(jié)點(diǎn)。從圖10 中可以看出，新人工邊界條件可以很好地模擬邊界處、材料界面以及地下結(jié)構(gòu)帶來的反射波和散射波，以及隨時(shí)間空間變化的輻射效應(yīng)。新邊界計(jì)算模型的場(chǎng)地反應(yīng)和結(jié)構(gòu)反應(yīng)均與參考解吻合較好，而粘彈性人工邊界誤差較大。驗(yàn)證了新人工邊界條件的有效性和精度。

圖10 含地下結(jié)構(gòu)場(chǎng)地中觀察點(diǎn)位移時(shí)程結(jié)果Fig. 10 Time histories of displacement solutions at observation points in a site containing underground structure

4 結(jié)論

本文針對(duì)豎向成層介質(zhì)中的標(biāo)量波傳播問題，基于連分式展開和擴(kuò)展的一致邊界，建立了含豎向成層介質(zhì)和地表不規(guī)則場(chǎng)地的折線形高精度人工邊界條件。半空間場(chǎng)地劃分為豎向成層介質(zhì)以及兩側(cè)無限域。兩側(cè)無限域的人工邊界條件基于動(dòng)力剛度的連分式展開，引入輔助變量來模擬無限域輻射效應(yīng)；豎向成層介質(zhì)為廣義的一致邊界，兩端開放與側(cè)邊連分式邊界耦合，并引入斜角變換使得提出的人工邊界條件可以為任意折線形。側(cè)邊界與底邊界耦合后的折線形人工邊界可用于豎向成層介質(zhì)中標(biāo)量波傳播分析。數(shù)值算例結(jié)果表明，提出的高精度人工邊界條件適用于多種豎向成層介質(zhì)中標(biāo)量波傳播問題，具有較高的計(jì)算精度，極大提高計(jì)算效率。

與水平成層半空間中標(biāo)量波傳播問題的H 形人工邊界[34]相比，本文折線形人工邊界條件主要是對(duì)原有方法的擴(kuò)展應(yīng)用。水平成層半空間模型轉(zhuǎn)換到豎向成層半空間模型的過程中，多層和半空間區(qū)域的無限延伸方向以及物理邊界均發(fā)生變化，人工邊界的截?cái)辔恢靡搽S之改變。因此針對(duì)水平成層半空間中的人工邊界條件不再直接適用，但具有重要參考價(jià)值。為進(jìn)一步擴(kuò)展人工邊界條件的適用范圍，適用于傾斜地層的人工邊界條件正在研究中。