流體飽和半空間二維地形三分量彈性波散射間接邊界元模擬

梁建文， 吳孟桃， 巴振寧*

1 天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室， 天津 300350 2 天津大學(xué)土木工程系， 天津 300350

0 引言

震害調(diào)查和強(qiáng)震觀測表明：局部地形(如凹陷、沉積和凸起)顯著影響地震動的幅值及空間分布(張建毅等, 2012).而近年來，眾多超高復(fù)雜建筑、重大交通工程和大型水利設(shè)施等不斷出現(xiàn)，這些位于場地條件復(fù)雜地區(qū)的防震減災(zāi)工作亟需精確可靠的地震動參數(shù).因此開展場地地震動局部地形效應(yīng)研究具有重要的理論價值和巨大的社會效應(yīng).

局部場地地震效應(yīng)，實(shí)質(zhì)上是場地對彈性波的散射問題，近40年來，引起了國內(nèi)外學(xué)者的重點(diǎn)關(guān)注.研究手段主要包括波函數(shù)展開解析解(Lee, 1984；袁曉銘和廖振鵬, 1993)和各種數(shù)值求解方法，例如有限元法(景立平等, 2005；宋貞霞和丁海平, 2013)、有限差分法(廖振鵬等, 1981；Zhang and Chen, 2006)、譜單元法(Komatitsch et al., 2004；周紅等, 2010)、邊界元法(Fu and Wu, 2001；Ge et al., 2005；梁建文和巴振寧, 2007)、廣義R/T遞推傳播矩陣法(Chen, 1990, 1995；Ge and Chen, 2007)、離散波數(shù)法(Bouchon et al., 1989；Fu and Bouchon, 2004)和各種逼近方法(Fu, 2005；Hu et al., 2009；Yu et al., 2010； 管西竹等, 2011；Yu and Fu, 2012, 2014；陳高祥等, 2017).

從地形散射機(jī)理來看，依據(jù)入射波方向與地形軸線的關(guān)系可分為二維散射(入射波場與地形垂直)和三維散射(入射波場與地形不垂直).前者自Trifunac(1971)開創(chuàng)性地給出半圓凹陷地形對均勻半空間SH波的散射后，Wong和Jennings(1975)、Sánchez-Sesma和Campillo(1991)、Yuan和Men(1992)、杜修力和熊建國(1993)、劉晶波(1996)、Fu和Wu(2001)、Fu等(2002)、李偉華和趙成剛(2004)、梁建文等(2006)、Tsaur(2011)、Gao等(2012)、陳少林等(2014)、Lee和Liu(2014)、巴振寧等(2017)、Shyu等(2018)、高玉峰(2019)等相繼對常規(guī)局部地形地震散射影響進(jìn)行了大量深入研究.相對于二維散射，地形對三維地震波散射的研究成果還較少，主要有Sánchez-Sesma等(1983)、Mossessian和Dravinski(1989)、Sánchez-Sesma和Luzón(1995)、Chávez-García(2003)、趙成剛和韓錚(2007)、梁建文等(2014)、王銘鋒等(2017)、Ba和An(2019)等采用各類方法求解了均勻或?qū)訝畎肟臻g中局部地形的三維響應(yīng)問題，獲得了二維模型無法反映的復(fù)雜場地三維動力響應(yīng)特征.三維模型與二維相比有著本質(zhì)差別，能夠全面反映地震動的空間分布特征，但對計算內(nèi)存需求大、求解耗時長，尤其對于時域響應(yīng)解.

在不失一般性的情況下，建立二維地形的三分量散射模型(即2.5維模型)是一種計算量遠(yuǎn)小于三維，并能在一定程度上反映場地三維動力響應(yīng)特征的有效方法.Pedersen等(1994)采用基本解為全空間格林函數(shù)的間接邊界元方法(IBEM)研究了二維沉積谷地的三維散射問題，但該方法需要離散自由地表.Stamos和Beskos(1996)采用直接邊界元方法(DBEM)求解了二維凹陷地形的三分量彈性波散射，其解也僅適用于均勻半空間情況.Ba和Liang(2010a, b)采用基本解為半空間格林函數(shù)的IBEM求解了層狀半空間凹陷地形的2.5維散射問題.以上研究僅限于單相介質(zhì)，而實(shí)際地表土層很多是飽和多孔介質(zhì)，尤其是濱海地區(qū)的沉積谷地中沉積層多為飽和土.根據(jù)Biot理論(Biot, 1956)，波在流體飽和介質(zhì)中的傳播與彈性介質(zhì)中的傳播特征截然不同.目前關(guān)于流體飽和半空間二維地形的三維地震散射問題研究，仍鮮有報道.

本文將Ba和Liang(2010a,b)給出的方法進(jìn)一步拓展到層狀流體飽和半空間情況，推導(dǎo)建立了飽和土層及下臥基巖的三維精確動力剛度矩陣和移動線源下的動力格林函數(shù)，進(jìn)而求解了流體飽和半空間中二維地形三分量彈性波散射問題.文中通過與已發(fā)表結(jié)果的比較驗(yàn)證了方法的可靠性和求解精度，并以梯形凹陷和半橢圓沉積地形為算例，在空間域內(nèi)給出了一些時頻數(shù)值結(jié)果，以期為揭示層狀介質(zhì)中局部地形地震波的三維散射機(jī)理提供一定參考.

1 理論模型及求解

如圖1所示，流體飽和層狀場地中存在一無限長任意截面形狀局部地形，場地模型由水平成層的飽和土及其下的基巖半空間組成，局部地形與土層的交界面為S.地震體波(P、SV或SH波)由下臥面入射，入射波波陣面的法向量在水平面上的投影與y軸的夾角為θh(水平入射角)，在豎直面內(nèi)與z軸的夾角為θv(豎向入射角).

圖1 流體飽和半空間二維地形三分量彈性波散射模型示意圖Fig.1 Sketch of three-component scattering of elastic waves by 2-D topography in a fluid-saturated half-space

由于彈性波斜入射下，半空間中局部地形沿其y軸的兩個相鄰截面的響應(yīng)只存在相位差(因截面位置不同產(chǎn)生)，那么可僅選取任一截面進(jìn)行離散和求解，然后通過沿地形軸向的波數(shù)拓展求得其他截面位置的動力響應(yīng).本文采用以移動斜線荷載(孔隙水壓)動力格林函數(shù)為基本解的2.5維IBEM進(jìn)行問題求解.具體求解時，通過波場分離將整體模型分解為開口的層狀半空間(波場1)和閉合的局部地形(波場2)，并假定波場1中存在著自由波場和散射波場，波場2中僅存在散射波場.首先提出層狀飽和場地自由波場(無地形存在)求解的精確動力剛度矩陣方法，然后通過在各離散斜線單元上施加一組沿軸線方向(y軸)移動的三坐標(biāo)方向虛擬均布荷載和虛擬孔隙水壓，求解位移和應(yīng)力的格林影響函數(shù)來模擬散射波場.虛擬移動均布荷載和孔隙水壓的密度，可通自由地表和地形的邊界條件(位移、應(yīng)力連續(xù)條件和透水條件)求得.對于時域響應(yīng)則可通過對頻域結(jié)果進(jìn)行快速Fourier逆變換得到.

1.1 Biot理論及其一般解

根據(jù)文獻(xiàn)(Biot, 1962)，流體飽和多孔介質(zhì)體系的土體本構(gòu)方程和滲流連續(xù)方程可表示為

σij=λ*eδij+2μ*εij-αδijpf,

(1a)

-pf=αMe+Mξ,

(1b)

以土骨架位移U和流體相對位移w表示的飽和土動力控制方程為

(2a)

(2b)

其中，ρ=(1-n)ρs+nρf，n為孔隙率，ρ、ρs和ρf分別為兩相土、土骨架和流體的質(zhì)量密度；b、m、α和M均為飽和相關(guān)參數(shù).

函數(shù)f(x,y,z,t)對時間t以及水平坐標(biāo)x和y的Fourier變換及其逆變換可表示為

(3a)

f(x,y,z)=

(3b)

其中，上標(biāo)“=”表示對空間分量(x,y)的Fourier變換，上標(biāo)“～”表示對時間t的Fourier變換.對式(2)進(jìn)行式(3a)的Fourier變換并進(jìn)行求解，可得到土骨架位移以及流體相對位移在頻率-波數(shù)域中解為

-kSzkxβ(ASVeikSzz-BSVe-ikSzz)-kyβ(ASHeikSzz+BSHe-ikSzz),

(4a)

-kSzkyβ(ASVeikSzz-BSVe-ikSzz)+kxβ(ASHeikSzz+ASHe-ikSzz),

(4b)

-(ASVeikSzz+BSVe-ikSzz)/β,

(4c)

-D3(ASVeikSzz+BSVe-ikSzz)/β,

(4d)

D3=ρfω2/(ibω-mω2)，AP1、BP1、AP2、BP2、ASV、BSV、ASH和BSH分別為上行和下行P1、P2、SV和SH波的幅值.kP1、kP2和kS分別為P1、P2和S波的波數(shù)，可由下式確定：

(5)

其中，

β1=[(mω2-ibω)(λ*+2μ*+αM)+ρω2M

-2αω2MρfB]/[M(λ*+2μ*)],

(6a)

(6b)

1.2 飽和介質(zhì)三維精確動力剛度矩陣及自由波場求解

考慮上行和下行的波幅值，每一飽和土層上下表面處的土骨架位移和流體相對于土骨架位移可由式(7a)表示.結(jié)合式(1)和式(4)，同時令飽和土層上表面處的外力

(7a)

[S]{AP1,BP1,AP2,BP2,ASV,BSV,ASH,BSH}T,

(7b)

(8)

(9)

1.3 移動線荷載(孔隙水壓)動力格林函數(shù)及散射波場模擬

流體飽和半空間層內(nèi)部作用的移動斜線荷載(孔隙水壓)如圖2所示，荷載作用層厚度為d，px0、py0和pz0分別為沿坐標(biāo)x、y和z方向的外荷載幅值，pf0為孔隙水壓幅值.該移動線源在時間-空間域中展開形式為：

p(x,y,z,t)=p0δ(z-xtanθ)δ(y-ct),

(10)

(11)

如圖2，在節(jié)點(diǎn)1、2處引入附加交界面以使荷載作用介質(zhì)層固定，總響應(yīng)涉及的特解、齊解和反力解可由直接剛度法方便求得.流體飽和半空間動力格林函數(shù)的詳細(xì)推導(dǎo)過程可參考層狀彈性介質(zhì)的思路(Ba and Liang, 2010a)，不再贅述.

圖2 流體飽和層中移動斜線均布荷載(孔隙水壓)Fig.2 Moving distributed forces (pore pressure) acting on an inclined line in fluid-saturated layers

×exp(iωt-ikxx-ikyy)dkxdky

(12)

利用該基本解，我們得到散射波場在地形邊界S各單元上產(chǎn)生的土骨架位移以及流體相對位移、應(yīng)力和孔隙水壓計算公式為

(13a)

(13b)

1.4 邊界條件及求解

求解中涉及到的邊界條件，按地形類型可分為：(1)零應(yīng)力邊界條件；(2)位移和應(yīng)力連續(xù)邊界條件.按是否透水則有：(a)排水邊界(透水情況)；(b)不排水邊界(不透水情況).注意到，對于透水情況，孔隙流體自由流動，孔隙水壓為零；對于不透水情況，孔隙流體無法自由流動，流體較土骨架的相對位移為零.而無論是否透水，求解步驟均相同，以下給出一般性邊界條件：

(1)零應(yīng)力邊界條件(如凹陷地形，以“Sc”表示)

(14)

其中[W(s)]為權(quán)函數(shù)，一般取為單位矩陣.將式(13)代入式(14)得：

(15)

其中，

(16a)

(16b)

(2)位移和應(yīng)力連續(xù)邊界條件(如沉積地形，以“Sv”表示)

(17a)

(17b)

同樣，將式(13)代入式(17)得：

(18a)

(18b)

其中，

(19a)

(19b)

(19c)

(19d)

求解式(15)和(18)，可得層狀飽和半空間中局部地形附近地表任意位置的位移和應(yīng)力(孔隙水壓)幅值，對于時域動力響應(yīng)可通過快速Fourier逆變換求得.

2 退化驗(yàn)證和精度比較

圖3 本文結(jié)果與Liang等(2006)給出結(jié)果的比較Fig.3 Comparison of the solutions given by the present study and Liang et al.(2006)

圖4 本文結(jié)果與Chen 等(2008)給出結(jié)果的比較Fig.4 Comparison of the solutions given by the present study and Chen et al.(2008)

de Barros和Luco(1995)采用邊界積分方程法首次給出彈性半空間中二維半圓沉積谷地三分量彈性波散射結(jié)果的數(shù)值解.在本文方法退化為單相介質(zhì)情況的基礎(chǔ)上，分別以均勻地形和層狀地形對P波散射的求解模型為例，圖5比較了本文方法和文獻(xiàn)方法求得的三個坐標(biāo)方向地表位移幅值在不同斜入射角度下的計算結(jié)果.相關(guān)參數(shù)取值如下：αv=2βv=βL，αL=2βL，αR=2βR=4βL，ρv=(2/3)ρL，ρR=(4/3)ρL，ζR=ζL=ζv=0.005，其中α和β分別為壓縮波波速和剪切波波速，上標(biāo)“v”、“L”和“R”分別表示與沉積內(nèi)、沉積外和下臥基巖有關(guān)的變量.無量綱頻率取為η=ωa/πβL=0.5，地震波入射角度取為θh=0°和θv=30°.由圖5可見，無論在均勻介質(zhì)還是在層狀介質(zhì)中，兩種方法的計算結(jié)果均具有良好的一致性，較好地驗(yàn)證了本文方法的求解精度.

圖5 本文結(jié)果與De Barros和Luco(1995)給出結(jié)果的比較Fig.5 Comparison of the solutions given by the present study and Barros & Luco (1995)

3 數(shù)值算例

3.1 凹陷地形對彈性波的三維散射

圖6 地表透水時梯形凹陷附近的孔壓變化(a) 透水：η=0.5, n=0.3, θh=0°, θv=30°; (b) 透水：η=0.5, n=0.3, θh=90°, θv=30°; (c) 透水：η=1.0, n=0.3, θh=0°, θv=30°; (d) 透水：η=1.0, n=0.3, θh=90°, θv=30°; (e) 透水：η=2.0, n=0.3, θh=0°, θv=30°; (f) 透水：η=2.0, n=0.3, θh=90°, θv=30°.Fig.6 Variation of pore pressure near the trapezoidal canyon with drained condition

圖7 地表不透水時梯形凹陷附近的孔壓變化(a) 不透水：η=0.5, n=0.3, θh=0°, θv=30°; (b) 不透水：η=0.5, n=0.3, θh=90°, θv=30°; (c) 不透水：η=1.0, n=0.3, θh=0°, θv=30°; (d) 不透水：η=1.0 n=0.3, θh=90°, θv=30°; (e) 不透水：η=2.0, n=0.3, θh=0°, θv=30°; (f) 不透水：η=2.0, n=0.3, θh=90°, θv=30°.Fig.7 Variation of pore pressure near the trapezoidal canyon with undrained condition

結(jié)果顯示，飽和排水情況(地表透水)和飽和不排水情況(地表不透水)的場地動力響應(yīng)顯著不同，且與入射波的頻率、角度及觀測點(diǎn)位置密切相關(guān).總體上，飽和不排水時凹陷表面的孔隙水壓幅值更大，原因在于孔隙流體的存在對波的傳播有著耗能的作用，而隨著觀測點(diǎn)逐漸遠(yuǎn)離凹陷地形(x/a逐漸增大)，兩者的幅值差異逐漸減小，地表是否透水的影響減弱.隨著入射頻率的增大，半空間中孔隙水壓幅值的變化及分布逐漸變得復(fù)雜.η=0.5(較低)時，幅值響應(yīng)主要由其自由波場決定，η=2.0(較高)時，起主導(dǎo)作用的則為凹陷地形引起的附加散射波場.隨著水平入射角度的變化，凹陷周圍土層孔隙水壓幅值的空間分布變化較為明顯.θh=0°時，孔隙水壓幅值關(guān)于凹陷軸線對稱，θh=90°時，孔隙水壓幅值在波入射側(cè)較大，且分布較為復(fù)雜.

表1 基巖上單一飽和土層中半橢圓沉積谷地材料參數(shù)Table 1 Material parameters of semi-elliptical alluvial valleys in a single saturated soil layer overlying the bedrock half-space

事實(shí)上，與經(jīng)典的半圓和半橢圓凹陷地形相比(Wong and Jennings, 1975；袁曉銘和廖振鵬, 1993；Liang et al., 2006)，梯形凹陷的土層和地形放大效應(yīng)研究還較少，本文僅給出了一個均勻凹陷三維散射問題的簡單結(jié)果，后續(xù)可針對成層場地中更多復(fù)雜截面形狀凹陷地形開展研究.已有結(jié)果表明，對于給定的入射頻率和入射角度，層狀飽和半空間中凹陷附近動力響應(yīng)由層狀場地自身的動力特性(成層場地自由波場)和凹陷截面形狀(凹陷地形形成的散射波場)共同決定(巴振寧和梁建文, 2013)，而均勻飽和半空間情況僅取決于凹陷截面形狀.因此為更為準(zhǔn)確的確定凹陷地形附近的地震響應(yīng)特征，需綜合考慮凹陷附近的土層分布情況及其截面形狀.

3.2 沉積地形對彈性波的三維散射

結(jié)果顯示，沉積地形的存在(圖中虛線位置)顯著改變了飽和層狀半空間中孔隙水壓振幅的分布.對于入射P1波，孔隙水壓幅值的空間分布十分依賴于豎向入射角θv，但對水平入射角θh相對不敏感，其幅值大小似乎隨著豎向入射角的增大而逐漸增大.孔隙水壓幅值的最大值出現(xiàn)在靠近谷地與土層交界面的右側(cè)，這意味著這些區(qū)域具有較大的應(yīng)變響應(yīng).當(dāng)θv=90°(P1波垂直入射)時，幅值的分布關(guān)于谷地軸線對稱，響應(yīng)與水平角度無關(guān).對于入射SV波，水平和豎向角度都對孔隙水壓的變化(尤其是幅值)有著顯著影響.比較圖8和圖9，可以明顯看出SV波入射下谷地內(nèi)部孔隙壓力的空間分布更復(fù)雜，幅值也更大.其原因在于不同類型入射波幅值的差異主要是自由波場的差異，孔隙流體承擔(dān)傳遞來的激勵能量改變，孔隙水壓自然不同.同樣，當(dāng)θv=90°(SV波垂直入射)時，幅值的分布具有對稱性，這也為我們的結(jié)果提供了額外的信心.

圖8 P1波入射下半橢圓沉積附近的孔壓變化(a) θv=45°, θh=30°； (b) θh=45°, θv=30°； (c) θv=45°, θh=60°； (d) θh=45°, θv=60°； (e) θv=45°, θh=90°； (f) θh=45°, θv=90°.Fig.8 Pore pressure amplitudes around the semi-elliptical valley under P1 waves incidence

圖9 SV波入射下半橢圓沉積附近的孔壓變化(a) θv=45°, θh=30°； (b) θh=45°, θv=30°； ( c) θv=45°, θh=60°； (d) θh=45°, θv=60°； (e) θv=45°, θh=90°； (f) θh=45°, θv=90°.Fig.9 Pore pressure amplitudes around the semi-elliptical valley under SV waves incidence

圖10 均勻飽和半空間中沉積附近地表位移幅值時域結(jié)果Fig.10 Time domain solutions for surface displacement amplitude around the valley in a uniform half-space

圖10展示了均勻半空間情況沉積附近地表位移時域響應(yīng)，圖中|Ux/ASV|、|Uy/ASV|和|Uz/ASV|分別表示三個坐標(biāo)方向的無量綱位移.可以看到，沉積地形的存在使得谷地附近的位移幅值與自由波場(無沉積地形存在)存在明顯差異，且由于局部不規(guī)則形成的反射波和散射波與入射波相互干涉，使得谷地附近位移幅值的持續(xù)時間明顯長于自由波場情況.另外，SV波在整個半空間的傳播過程在時域內(nèi)可更為清楚的觀察到，如谷地左側(cè)(x/a≤-1.0)首先接收到直達(dá)波，隨后接收到谷地左角點(diǎn)(x/a=-1.0)產(chǎn)生的散射波，最后接收到谷地右角點(diǎn)(x/a=1.0)產(chǎn)生的散射波.同時，入射波在較軟的沉積介質(zhì)中經(jīng)過多次反射，導(dǎo)致谷地內(nèi)部位移相較于外部的持續(xù)時間明顯更長，且幅值更大.

4 結(jié)論

基于Biot理論，建立了層狀多孔介質(zhì)內(nèi)部移動線荷載(孔隙水壓)格林函數(shù)，進(jìn)而提出了求解流體飽和半空間中二維局部地形三分量彈性波散射問題的間接邊界元方法(IBEM).該方法優(yōu)勢在于離散僅限于地形邊界(無須離散自由地表)，且格林函數(shù)計算不存在奇異性(荷載直接加在邊界上)，因而計算精度容易控制，且對復(fù)雜邊界條件具有很強(qiáng)的適應(yīng)性.

文中以給出的方法具體研究了梯形凹陷地形和半橢圓沉積地形的三維地震效應(yīng)，分別在頻域和時域內(nèi)進(jìn)行了計算分析.頻域結(jié)果表明，局部地形的存在顯著改變了飽和層狀半空間中孔隙水壓的幅值及空間分布，且依賴于地震波類型、入射頻率、入射角度和邊界透水條件等影響因素.地形附近動力響應(yīng)由半空間動力特性和地形截面形狀共同決定.時域結(jié)果表明由于局部不規(guī)則形成的反射波和散射波與入射波相互干涉，使得地形附近位移幅值的持續(xù)時間明顯長于自由波場情況.與均勻半空間相比，基巖的存在的位移時程的空間分布更為復(fù)雜，振動時長增加.層狀半空間，土層自身特性(土層剛度、厚度等)對波的傳播路徑和持續(xù)時間有重要影響.

本文建立的以移動斜線荷載(孔隙水壓)格林函數(shù)為基本解的IBEM在近地表復(fù)雜場地對地震波三維散射和土-結(jié)構(gòu)相互作用方面有著重要應(yīng)用前景.后續(xù)研究將進(jìn)一步考慮有實(shí)際觀測數(shù)據(jù)的場地，以期為區(qū)域地震危險性分析、烈度區(qū)劃以及局部地形附近建筑物的抗震設(shè)防等工作提供一定理論基礎(chǔ).