李志遠，鐘 紅，胡志強，林 皋

(1. 海岸與近海工程國家重點實驗室，大連理工大學，大連 116024；2. 工程抗震研究所，建設工程學部水利工程學院，大連理工大學，大連 116024；3. 中國水利水電科學研究院工程抗震研究中心，北京 100048)

多次震害調(diào)查表明，地下結(jié)構(gòu)對地震波有散射、折射、衍射等作用，顯著改變了地表的地震響應。水平巖層在構(gòu)造運動作用下會產(chǎn)生彎曲，形成局部褶皺。現(xiàn)有研究多假定場地為均質(zhì)或水平成層場地，忽略了局部復雜地形對地下結(jié)構(gòu)散射場的影響。因此研究局部褶皺對層狀地基中地下結(jié)構(gòu)散射作用的影響對地震安全性評價具有重要指導意義。

天然的工程場地常被假定為水平成層分布。Thomson 和Haskell[1 ? 2]提出的傳遞矩陣算法，可以提供分層介質(zhì)中彈性波動控制方程的系統(tǒng)解。為了解決求解時指數(shù)溢出的問題，Kausel 和Roesset[3]重新推導了力和位移關(guān)系，并建立了精確剛度法，Chen[4]定義了全局廣義的反射-透射矩陣。但是，這類算法只適用于求解具有規(guī)則邊界條件的水平成層地基，對包含不規(guī)則邊界的場地不能直接使用。Chen[5]利用模態(tài)綜合法研究了不規(guī)則層狀地基中Love 波的傳播特性，Lee 等[6 ? 7]采用加權(quán)余量法求解了體波和面波在不規(guī)則層狀地基的傳播特性。值得指出的是以上研究假定不規(guī)則土層均限于有限區(qū)域，遠場仍為水平成層場地，且不包含孔洞等局部構(gòu)造。Luco 和Apsel[8]基于傳遞矩陣法提出了一種層狀地基的格林函數(shù)，在此基礎(chǔ)上得到了移動荷載在層狀地基的動力響應[9]，利用層狀地基的移動荷載的格林函數(shù)，構(gòu)造了一種間接邊界元法，可以將求解三維結(jié)構(gòu)的動力響應。利用這種間接邊界元法，Barro 和Luco[10]分析了SV 波、P 波入射時，層狀地基中埋置孔洞的對地表幅值的響應。他們還分析了層狀半空間存在襯砌結(jié)構(gòu)的隧洞時，入射SV 波、P 波，孔洞襯砌的動應力集中系數(shù)[11]。Wuttke 等[12]結(jié)合離散波數(shù)法和邊界元法給出了一種兩步法用于分析層狀地基中波動傳播。國內(nèi)學者尤紅兵和梁建文[13]利用間接邊界元法求解了層狀地基中洞室對SV 波的散射。巴振寧等[14]利用間接邊界元法分析了場地動力特性對襯砌隧道地震反應的影響。趙密等[15]提出一種土-結(jié)構(gòu)相互作用的高效方法，提高了深厚土層地下結(jié)構(gòu)地震分析的計算效率。

值得指出的是以上文獻求解散射場是通過對散射波和自由場波動的疊加而得到，這種方法將地下結(jié)構(gòu)與周圍土體作為一個整體分析，散射波的求解通常需要對散射體施加邊界條件，當結(jié)構(gòu)復雜時處理起來會變的異常復雜，因此以上研究的散射體均假設為圓形、橢圓形等?？紤]多種復雜條件，如包含局部褶皺層狀地基中馬蹄形孔洞散射波的高精度解幾乎沒有。

本文基于Lin 等[16]提出散射場分析模型，通過一條人工邊界將層狀介質(zhì)中復雜邊界條件的波動散射問題轉(zhuǎn)化為子結(jié)構(gòu)波動響應求解的問題。這種轉(zhuǎn)化對于線彈性介質(zhì)是準確的，沒有引入任何假設，當輻射邊界條件和波動輸入為高精度解時，復雜邊界的散射問題也是高精度的?；诒疚姆椒ǎ懻摼植狂薨檶訝顖龅刂旭R蹄形孔洞散射場的影響。

1 模型及計算方法

1.1 子結(jié)構(gòu)法控制方程

圖1 為典型的土-結(jié)構(gòu)相互作用體系，虛線ABCD 包圍的區(qū)域為近場復雜結(jié)構(gòu)，由虛線ABCD一直延伸到無窮遠處的部分為遠場，虛線ABCD交界面。本文只列出關(guān)鍵公式，詳細推導見文獻[17]。頻域中子結(jié)構(gòu)法的控制方程為：

圖1 土-結(jié)構(gòu)相互作用體系示意圖Fig.1 Soil-structure interaction system

式中： K 和 M分別為通過有限元計算得到剛度陣和質(zhì)量陣； ω為圓頻率。注意到當近場結(jié)構(gòu)為規(guī)則開挖土體時，整個體系退化為具有規(guī)則邊界的自由場，利用式(1)對自由場體系進行分析，可以得到如下關(guān)系：

將式(3)代入式(1)中可得：

式(4)即為子結(jié)構(gòu)法的控制方程。由式(4)可以看出，欲求復雜場地的動力響應，需要計算三部分內(nèi)容：自由場交界面處的動力剛度和波動位移；近場子結(jié)構(gòu)和開挖土體的動力剛度；控制方程的求解。

1.2 層狀地基輻射邊界

地下隧洞等結(jié)構(gòu)沿軸向長度遠大于截面尺寸，因此可以考慮為二維平面應變問題。建立直角坐標系，沿水平向右為x 軸正向，垂直向下為z 軸正向。層狀地基的材料參數(shù)為：拉梅常數(shù) λ和G，密度為 ρ，泊松比ν。忽略體力，用位移表示的波動方程為：

引入應力作為對偶變量可將方程化簡為一階常微分方程：

式中：V=[ uxux?σxz?σzz]T； H為水平波數(shù)κ、頻率 ω和地基材料參數(shù)的函數(shù)[18]。任選一土層分析，定義土層[za,zb]，上層面為 za，下層面為zb，土層的厚度為?，無外力作用時，土層兩端位移向量ub和應力向量 pa，滿足如下方程：

相鄰土層在層面處滿足位移連續(xù)和應力連續(xù)條件。相鄰土層1[za,zb]和土層2[zb,zc]都可表示成對偶方程式(7)的形式，利用層面間位移和應力的連續(xù)條件，消去層面b 處的位移和應力向量，可得到合并后的土層c[za,zc]的對偶方程，得到系數(shù)矩陣間的關(guān)系為：

采用精細積分計算，土層厚度可以取的非常小(如τ=?/220)，這時初始區(qū)段的系數(shù)矩陣可通過Taylor 級數(shù)展開獲得[19]。當土層厚度非常小時，F(xiàn)趨于單位陣，式(8)不再適用。因此針對微層間的合并時，式(8)做如下修改：

其中，L=(I+G1Q2)?1。

由初始土層進行合并，可得到任意層面間的對偶關(guān)系。下面采用子結(jié)構(gòu)法求解頻域-波數(shù)域中的整體剛度矩陣。首先對層狀地基在豎直方向進行離散，在計算點處設置土層層面，其次對相鄰計算點所在土層進行合并，以區(qū)段[za,zb]為例，可得到層面間的對偶方程(式(7))。為便于使用子結(jié)構(gòu)法，將對偶方程轉(zhuǎn)化成剛度形式：

式中：

或者寫成緊湊形式：

式中： P 表示外荷載的向量； K為區(qū)段的剛度矩陣；U 為位移向量。

最后是整體剛度矩陣的組裝。在這種情況下，全局載荷矢量對應于界面處的外部作用力，整體剛度的組裝可應用結(jié)構(gòu)動力學子結(jié)構(gòu)法的定理和技術(shù)。整體的剛度矩陣通過在層狀體系的每個“節(jié)點”(層面)處疊加相鄰土層的貢獻來獲得：

式中： R∞為半空間層面處的動力剛度[19]；n 為土層個數(shù)。

1.3 自由場波動輸入

如圖2 所示，無限半空間上臥n?1 個均勻彈性層，有n 層層面，圖中 ρj、 λj、Gj表示第j 層的密度和拉梅常數(shù)。每層單獨建立坐標系，坐標原點設在該層的頂面上，水平向右為x 軸正向，豎直向下為z 軸正向。

圖2 層狀地基中波動傳播Fig.2 Wave propagation in a layered half-space

平面內(nèi)右行簡諧波的一般解可寫成：

式中： φ、 ψ分別為SV 波、P 波的位勢函數(shù)； κ1為水平向波數(shù)； κ2S、 κ2P分別為SV 波和P 波的豎直向波數(shù)； ES、 EP分別為SV 波、P 波上行波幅值；FS、 FP分別為SV 波、P 波下行波幅值。位移滿足如下關(guān)系：

將式(14)代入式(15)可得各層面處的位移表達式；由位移-應變關(guān)系式，以及應力-應變關(guān)系式可以得到各層面處的應力表達式。引入應力-位移狀態(tài)向量，整理后可表示為：

式中：

頂層的波幅向量可表示為：

式中：RSP、 RSS為SV 波的反射系數(shù)；RPS、 RPP為P 波的反射系數(shù)，可通過地表自由面的邊界條件求得[20]。

下面建立半空間和頂層波幅向量間的關(guān)系，并求解頂層和半空間的波幅向量。由第2 節(jié)可知，層狀體系可以通過精細積分進行高效合并。利用精細積分合并層狀體系，最終得到界面zn?1和界面 z0間的對偶方程：

式中，F(xiàn)1n、G1n、Q1n為合并層狀體系后得到的系數(shù)矩陣，結(jié)合式(16)可得半空間和頂層的波幅向量關(guān)系：

其中：

式(30)寫成緊湊形式：QEE=[amk]， (m,k=1,2,···,4)。將式(20)代入式(23)可得：

其中：

當P 波入射時，EPn=1， ESn=0；當SV 波入射時，ESn=1， EPn=0。求解式(24)可以得到頂層的上行波的幅值和半空間中下行波幅值。將頂層上行波系數(shù)和反射系數(shù)帶入式(20)可得到頂層的波幅向量HE1，將頂層波幅向量代入式(16)可得頂層的狀態(tài)向量，其他各層的狀態(tài)向量可由精細積分方法合并得到兩端邊值的對偶關(guān)系，最終可求解。本文提出的層狀自由場的求解方法，可稱為改進傳遞矩陣法，與傳統(tǒng)傳遞矩陣法的主要區(qū)別在微層傳遞矩陣，即式(23)，其優(yōu)點在于基于區(qū)段混合能的表示，使 Q和 G對于層厚很小時也不會發(fā)生數(shù)值病態(tài)[21]，保證了求解的順利進行。同時，改進傳遞矩陣法的便于推廣到各向異性層狀地基波動問題的求解。

2 方法驗證

計算均質(zhì)半空間中垂直入射SV 波時地表位移幅值的響應，與文獻[11]的結(jié)果對比驗證本文方法的精度。如圖3 所示，均質(zhì)半空間存在淺埋圓柱形孔洞，均質(zhì)半空間的密度ρ0=2200 kg/m3，泊松比ν=0.333，剪切波速cS0=600 m/s。圓柱形孔洞半徑r=6.0 m，埋置深度為H=1.5r=9.0 m。虛線ABCD 包括的部分為近場區(qū)域4r×3.5r。定義位移響應幅值與入射波幅值的比值為無量綱位移幅值。垂直入射SV 波，入射波頻率為η=ωr/πcS0=2.0，圖4 為圓柱形孔洞場地地表的位移幅值，計算結(jié)果與Luco 和Barros[11]采用間接邊界元法的計算結(jié)果對比，兩者吻合性良好。

圖3 均質(zhì)半空間埋置圓形孔洞Fig.3 Schematic view of the shallow embedded circular cavity

圖4 淺埋圓柱形孔洞場地垂直入射SV 波地表位移幅值Fig.4 Surface displacement generated from scattering by a shallow embedded circular cavity for vertically incident SVwave with particle motion along x-axis

3 算例與分析

地質(zhì)運動作用下，褶皺是一種常見的工程地形，本節(jié)討論層狀半空間與上臥土層交界面為正弦波形狀時的對馬蹄形孔洞散射作用的影響。本文假定不規(guī)則土層交界面為有限區(qū)域，遠場仍為水平成層地基。

計算模型參數(shù)為：厚度為H=4r的單層土體，置于半空間之上；馬蹄形孔洞的穹頂半徑為r=10 m，邊墻高度為0.785r=7.85 m，埋置深度D=1.5r=15 m。半無限空間的材料參數(shù)：質(zhì)量密度ρ0=2000 kg/m3，泊松比ν0=0.35，剪切波速cS0=700 m/s，阻尼比ζ0=0.02；上臥土層的密度ρ1=2000 kg/m3，泊松比ν1=0.35，阻尼比ζ1=0.05，剪切波速cS1=350 m/s2。定義無量綱頻率η=ωr/(πcS1)=1.0。近場場地采用有限元計算近場模型，最大的單元尺寸為最小波長λmin=2r/ηmax=r的1/8。

如圖5 所示，討論兩種不同的工況，工況A 是褶皺埋入半空間中，工況B 是褶皺埋入土層中。實際工程中褶皺可能是各種形狀，本文為方便描述，討論正弦形褶皺，褶皺形狀分布滿足如下方程：

圖5 局部褶皺層狀地基中馬蹄形孔洞Fig.5 Horseshoe-shaped cavity embedded in a layered halfspace with local folds

褶皺的中心點與馬蹄形孔洞的中心在同一豎直線上，褶皺的分布范圍為x=?1.5r 和x=1.5r之間。上臥土層的厚度為H=4r，馬蹄形孔洞的埋置深度D=2.0r。入射SV 波，入射波頻率η=2，圖6 為工況A 和工況B 地表的位移幅值，作為參考，水平成層場地(即不考慮褶皺地形)的地表位移幅值也繪制在圖中。保持參數(shù)不變，入射頻率η=2的P 波，計算結(jié)果繪制在圖7 中。從圖6 中可以看出，SV 波入射時局部褶皺對地表位移幅值影響明顯；在地表?r

為了進一步研究局部褶皺對散射波的影響隨頻率變化的特性，圖8 中繪制了SV 波入射時觀測點P1 (x/r=1.5) 的位移頻譜。從圖8 可以看出，在低頻段(η<0.75)局部褶皺對地表水平向和豎直向位移幅值的影響均較弱；隨著頻率的增加，局部褶皺對地表位移幅值的影響增強，且不同的局部褶皺幾何構(gòu)造對地表位移幅值影響的程度不同。從整體上觀察位移頻譜的分布情況，可以看出：在低頻段(η<0.75)位移幅值有明顯的共振現(xiàn)象，可以認為，此時局部褶皺和馬蹄形孔洞對散射場影響均較弱，這時的波動響應接近自由場的波動響應；在0.75<η<1.5區(qū)段，位移幅值較小，也沒有明顯的共振現(xiàn)象，這時馬蹄形孔洞是影響散射場分布的主要因素；在1.5<η時，三種工況的位移幅值分布顯著不同，局部褶皺的影響明顯，此時局部褶皺已經(jīng)顯著的改變了位移幅值峰值和共振頻率。圖9 為P 波入射時觀測點P1 (x/r=1.5) 的位移頻譜。與SV 波入射時相比，具有如下特點：局部褶皺對地表位移幅值影響變顯著的起始頻率更低，SV 波出現(xiàn)在η=0.75附近，P 波出現(xiàn)在η=0.5附近；局部褶皺對地表位移幅值的影響不是簡單的隨著頻率的增加而增強，變化更加復雜，比如在1.5<η<1.75區(qū)間，三種工況的位移幅值比較接近，此頻率區(qū)段局部褶皺對位移幅值的影響要低于頻段1.0<η<1.5。

圖6 SV 波入射時地表位移幅值Fig.6 Surface displacements caused by SV wave

圖7 P 波入射時地表位移幅值Fig.7 Surface displacements caused by P wave

圖8 SV 波入射時觀測點P1 的位移頻譜 (P1 (x/r=1.5))Fig.8 Displacement spectrums of observation point P1 caused by SV wave (P1 (x/r=1.5))

圖9 P 波入射時觀測點P1 的位移頻譜 (P1 (x/r =1.5))Fig.9 Displacement spectrums of observation point P1 caused by P wave (P1 (x/r=1.5))

場地地震安全性評價中，位移峰值是一個重要參數(shù)。因此，圖10 繪制了SV 波入射時地表各處位移峰值分布圖，縱坐標表示頻率在η ∈[0,2]區(qū)間中位移幅值的峰值。從圖10 可以看出，工況A 的褶皺幾何構(gòu)造增大了地表水平向位移峰值，工況B 的褶皺幾何構(gòu)造降低了地表水平向位移峰值；對于豎直向位移幅值，在近馬蹄形孔洞區(qū)域(?r

圖10 SV 波入射時地表位移峰值Fig.10 Peak surface displacements caused by SV wave

圖11 P 波入射時地表位移峰值Fig.11 Peak surface displacements caused by P wave

4 結(jié)論

本文提出了一種求解近區(qū)域包含復雜構(gòu)造的層狀場地波動響應的算法。基于子結(jié)構(gòu)法建立了復雜場地散射問題的控制方程，將復雜場地散射問題轉(zhuǎn)化具有規(guī)則邊界條件層狀地基的動力剛度求解和波動響應求解。這種轉(zhuǎn)化對于線彈性介質(zhì)是完全準確的，因此當具有規(guī)則邊界條件層狀地基的輻射問題和波動輸入問題得到高精度解答時，復雜場地散射問題的解答也是高精度的。通過與現(xiàn)有文獻結(jié)果對比，驗證了本文方法的精度，并利用該方法，分析了SV 波和P 波入射時，局部褶皺對層狀馬蹄形孔洞場地的散射場影響。

研究結(jié)果表明，局部褶皺對馬蹄形孔洞的散射場具有顯著影響。在低頻段局部褶皺對SV 波入射時馬蹄形孔洞散射場的影響不明顯，隨著頻率的增加，局部褶皺的影響增大；局部褶皺對P 波入射時馬蹄形孔洞散射場的影響更加復雜，在某些高頻段，也出現(xiàn)了局部褶皺對散射場影響不顯著的情況。在SV 波或P 波入射時，在近馬蹄形孔洞區(qū)域，局部褶皺對豎直向位移峰值影響明顯，在距離馬蹄形孔洞稍遠的位置，局部褶皺對豎直向位移峰值幾乎沒有影響；整體分析可得，地表豎直向位移峰值受局部褶皺和馬蹄形孔洞的共同影響；地表水平向位移峰值受褶皺地形的影響更顯著。本文僅分析了兩種簡單的局部褶皺構(gòu)造對散射場的影響，但本文方法可用于任意近區(qū)域復雜構(gòu)造對地下結(jié)構(gòu)散射場影響的分析。需要指出，本文假定孔洞和褶皺均滿足平面應變理論，若考慮場地三維散射效應，需將控制方程擴展到2.5 維。