豎向成層介質(zhì)中標(biāo)量波傳播問(wèn)題的高精度人工邊界條件

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

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

在地下結(jié)構(gòu)抗震分析中，有限元法因其靈活高效的特點(diǎn)，被廣泛用于有限域模擬。為了更有效地反映截?cái)酂o(wú)限土體的輻射阻尼效應(yīng)，需要引入人工邊界條件[1]。人工邊界的適用性和精度直接影響計(jì)算結(jié)果的精度，并通過(guò)影響計(jì)算區(qū)域的大小而影響計(jì)算效率。城市抗震防災(zāi)和重大工程項(xiàng)目建設(shè)中經(jīng)常遇到豎向成層介質(zhì)以及不規(guī)則地形場(chǎng)地。在我國(guó)西南部山嶺地區(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)域引起越來(lái)越多的關(guān)注。建立豎向成層介質(zhì)中的高精度人工邊界條件具有研究和實(shí)際意義。國(guó)內(nèi)外研究者已提出大量人工邊界條件[9]，其中包括：粘性邊界[10]、粘彈性邊界[11]、透射邊界[12]、無(wú)限元法[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)地。

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

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

1 問(wèn)題描述

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

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

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

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

2 高精度人工邊界條件

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

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

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

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

2.2 豎層底部廣義一致邊界

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

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

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

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

3 數(shù)值算例

本節(jié)給出不同地形、地質(zhì)條件下四種場(chǎng)地中波傳播問(wèn)題的數(shù)值算例來(lái)驗(yàn)證提出高精度折線形人工邊界條件的有效性和精度。參考解為大區(qū)域有限元模型的時(shí)域解，模型尺寸足夠大能夠保證在觀測(cè)時(shí)間內(nèi)結(jié)果不受到截?cái)噙吔缣幏瓷洳ǖ挠绊?。同時(shí)采用了簡(jiǎn)單常用的粘彈性人工邊界條件[11]進(jìn)行了對(duì)比分析，說(shuō)明提出的折線形人工邊界條件適用于豎向成層不規(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è)豎向無(wú)限延伸地層，如圖3 所示。中間地層的剪切波速c1=100 m/s，兩側(cè)無(wú)限域介質(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í)間步長(zhǎng)為0.01 s。

圖3 含豎直地層半空間場(chǎng)地中標(biāo)量波傳播問(wèn)題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 中虛線所示。下文中圖例意義相同。從圖中可以看出，無(wú)論是加載點(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è)無(wú)限域介質(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í)間步長(zhǎng)為0.002 s。

圖5 含豎向成層介質(zhì)且地表起伏的半空間場(chǎng)地中標(biāo)量波傳播問(wèn)題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ì)算誤差較大。表明提出的人工邊界條件可以很好的模擬邊界處和材料變化帶來(lái)的反射和散射波，以及隨時(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è)無(wú)限域介質(zhì)剪切波速cl= 100 m/s，三個(gè)豎向地層的剪切波速分別為c2= 150 m/s、c3= 200 m/s、c4=250 m/s，右側(cè)無(wú)限域介質(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í)間步長(zhǎng)為0.002 s。

圖7 含豎向成層介質(zhì)的階梯地形半空間場(chǎng)地中標(biāo)量波傳播問(wèn)題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)生的反射和散射波，模擬無(wú)限場(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)量波傳播問(wèn)題。從左到右5 種介質(zhì)剛度依次增大，左側(cè)無(wú)限域介質(zhì)剪切波速cl= 100 m/s，三個(gè)豎向?qū)拥募羟胁ㄋ俜謩e為c2=150 m/s，c3= 200 m/s，c4= 250 m/s，右側(cè)無(wú)限域介質(zhì)剪切波速c5= 300 m/s。因含有地下結(jié)構(gòu)，計(jì)算模型2 不再適用或者說(shuō)誤差較大，本節(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í)間步長(zhǎng)為0.002 s。

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

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