巴振寧 ，張家瑋 ，梁建文 ，吳孟桃

(1.中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點實驗室(天津大學)，天津 300350；2.濱海土木工程結構與安全教育部重點實驗室，天津300350；3.天津大學土木工程系，天津 300350)

場地自由場地震反應是研究地震波散射和土-結構動力相互作用的基礎[1－2]，因此，該課題一直是巖土工程、地震工程和地震學等領域研究的重點。目前，國內外學者針對彈性波在不同介質中的傳播問題已開展了較多研究。對于各向同性彈性介質，Haskell[3]和 Thomson[4]最先給出了層狀介質中波傳播問題的傳遞矩陣解；隨后Kausel和Ro?sset[5]采用Haskell-Thomson傳遞矩陣的方法進一步推導了層狀半空間動力剛度矩陣，并通過改變波型和控制點的位置，用剛度矩陣方法研究了層狀彈性半空間場地豎直和水平運動分量變化對自由場的影響；Wolf和Obernhuber[6]推導了層狀半空間精確動力剛度矩陣，并采用剛度矩陣的方法求解了不同波入射時層狀半空間時的動力響應；于國友等[7]應用Haskell-Thomson傳遞矩陣方法，給出了不同性質場地對地震波的放大作用曲線；梁建文和巴振寧[8]將Wolf動力剛度矩陣推廣至三維，分析了基巖上單一土層場地的動力響應。對于各向同性飽和多孔介質，Deresiewicz和Rice[9]推導了含粘性液體飽和多孔介質 Biot動力平衡方程的通解，并給出了相速度、衰減系數(shù)、反射角和振幅比的解析表達式；Deresiewicz[10]推導了含有液體的多孔半空間中表面波的速度衰減系數(shù)表達式；Yang和 Sato[11]基于Biot兩相介質理論，模擬了飽和場地的地震反應，分析了飽和效應對單土層豎向地震激勵的影響；Lin[12]研究了孔隙率、固結狀態(tài)、邊界排水等多種因素對無粘流體飽和多孔半空間動力響應的影響；李偉華等[13]應用顯式有限元方法，分析了含有飽和軟弱土層場地的動力響應；巴振寧和梁建文[14]在頻域內求解了流體飽和半空間中埋置球面波源(球面膨脹P1和P2波源以及剪切SV波源)的動力響應。

值得指出的是，上述研究均基于各向同性介質(彈性或飽和)的假定。然而，天然土體具有固有各向異性和應力各向異性已逐漸成為巖土工程界的一種共識。由于風化沉積作用，天然土體的水平彈性模量往往大于其豎直彈性模量[15]，因此，采用橫觀各向同性(transversely isotropic，簡稱 TI)介質力學模型來描述土體的各向異性更為合理。對于彈性波在TI介質的傳播問題，陳镕等[15]建立了層狀TI場地模型，利用剛度矩陣的方法計算了其動力響應；薛松濤等[16]分析了 TI場地在三種不同底部邊界條件下的動力響應?；?Biot兩相介質理論，Schmitt[17]推導了柱坐標系下 TI飽和介質的波動方程，研究了柱面波的頻散和衰減曲線及與參數(shù)相關的靈敏度系數(shù)；Sharma[18]推導了直角坐標系下平面諧波的特征方程，并給出了三種準波波速的解析表達式；Liu等[19]考慮了流體粘度，研究了TI飽和多孔介質中的波傳播問題；陳勝立和張建民[20]建立TI飽和半空間模型，并給出了受荷載作用時場地的動力響應；陸建飛等[21]建立了軸對稱情況下 TI飽和土層的反射、透射矩陣，求解了集中力作用下TI飽和場地的頻域響應。然而，目前很少有針對地震波斜入射下TI飽和場地時域響應分析的研究。

本文將Haskell-Thomson傳遞矩陣方法拓展到層狀 TI飽和半空間，通過層間位移、應力連續(xù)條件及地表邊界條件建立了傳遞矩陣方程，進而結合快速傅里葉變換，求解了地震波斜入射下層狀 TI飽和場地的時域自由場反應。文中通過將退化后場地的計算結果與已有文獻結果對比驗證了方法的正確性，并進行了不同工況下的計算分析，研究了土體 TI性質及飽和特性對自由場響應的影響，得到了一些有益結論。

1 層狀TI飽和半空間模型及其求解

層狀TI飽和半空間模型如圖1所示。模型中各土層的厚度為dj(j=1,2,…,N)，相對應的土層下表面深度為zj(j=1,2,…,N)。土層和半空間性質由五個工程常數(shù)(Eh、Ev、Gv、νh、vvh)和飽和相關參數(shù)(m1、m3、α1、α3、k1、k3、η、?)以及質量密度ρ和阻尼比ζ進行描述。設波入射位置在基巖露頭，入射方向與豎直方向夾角為θ。

具體求解時：首先，根據(jù) TI飽和多孔介質動力平衡方程及孔隙流體運動方程求得其基本解；然后，根據(jù)層間連續(xù)條件建立傳遞矩陣，聯(lián)立地表邊界條件求得各層土中來波和去波的幅值；最后，代入基本解即可得到層狀TI飽和半空間自由場響應。

圖1 層狀TI飽和半空間模型Fig.1 The model of multi-layered TI saturated half-space

1.1 TI飽和多孔介質波動方程及其基本解

根據(jù)文獻[19]，TI飽和多孔介質平面內應力-應變關系為：

式中：σx、σz、τzx=τxz為總的應力分量；p為孔壓；ux和uz為土骨架的位移；wx和wz為流體相對位移；A11、A13、A33、A44和M1、M3均為與兩相材料有關的參數(shù)，其具體表達詳見附錄A。

直角坐標系下[22]，以u-w形式表示的Biot動力平衡方程為：

式中：ρm為 TI飽和多孔介質密度ρ=(1-?)ρs+?ρf，ρs和ρf分別為土體密度和孔隙流體密度；t為時間變量。

討論簡諧彈性波的傳播，不失一般性可取解的形式為：

式中：l1=sinθ，l3=cosθ；ax、bx、az和bz為波幅值；k=ω/c為水平波數(shù)，ω為振動圓頻率，c為視速度。

將式(4)代入式(2)和式(3)并考慮 qP1、qP2和qSV波所引起的位移之間的耦合關系，可得位移表達式如下：

式中：χ1~χ6、δ1~δ6具體表達式詳見文獻[23]；AqP1、AqP2和AqSV可視為在z軸負向傳播的qP1波、qP2波和qSV波的幅值；BqP1、BqP2和BqSV可視為在z軸正向傳播的qP1波、qP2波和qSV波的幅值；λqP1、λqP2和λqSV分別對應qP1波、qP2波和qSV波的波數(shù)，其具體表達式如下：

式中，Δ、Λ、?、d1及d2具體表達式見附錄B。

將式(5)代入式(1)中，可得到平面內應力及孔壓表達式如下：

1.2 層狀TI飽和半空間平面內傳遞矩陣

如圖1所示，結合式(5)以及式(7)并忽略e-ikx，可得第j層平面內部分位移、應力及孔壓與上、下行波向量的關系如下：

式中，S為系數(shù)矩陣，矩陣元素詳見附錄C。

如圖2所示，在土層交界面處，應力、位移及孔壓有以下連續(xù)條件：

圖2 TI飽和土層層間連續(xù)條件Fig.2 Continuity conditions between layers of TI saturated soil

將式(9)代入式(8)可得第j層上、下行波向量與第j+1層上、下行波向量的關系：

進一步變換式(10)并簡化表達，可得第j層上、下行波向量幅值與第j+1層上、下行波向量幅值之間的關系：

式中，E(hj)(j)矩陣元素詳見附錄D。

由式(11)可得第j層上、下行波向量幅值與第j+1層上、下行波向量幅值之間的關系：

式中，T(j)=[S(j+1)]-1S(j)E(hj)(j)。

進一步，可以推得第N+1層上、下行波向量幅值與第1層上、下行波向量幅值之間的關系：

1.3 邊界條件及地表動力響應

對于波在飽和層狀半空間中傳播問題的研究，由于傳播介質為固液兩相材料，建立邊界條件時，不但要考慮土體的性質，還要考慮土體與孔隙流體的相互作用[24－25]。

假定半空間表面完全透水，邊界條件可表示為：

假定半空間表面不透水時，邊界條件可表示為：

分別將式(5)和式(7)代入透水邊界條件和不透水邊界條件，可得地表土層中入射波幅值與反射波幅值的關系：

式中：當邊界透水時，M1~M9具體表達式見附錄E；當邊界不透水時；M1~M9具體表達式見附錄F。

聯(lián)立式(13)及式(16)，可求得地表處入射波及反射波幅值AqP1、AqP2、AqSV、BqP1、BqP2和BqSV，將其代入式(5)和式(7)可分別得到地表處位移及應力響應。

2 方法驗證

為了驗證本文頻域結果的正確性，將退化結果與文獻[26]結果進行比較。圖3(a)和圖3(b)分別給出了P波和SV波入射下各向同性彈性半空間中位移隨深度的變化規(guī)律。計算參數(shù)取固體顆粒密度ρs=2000 kg?m-3，剪切模量G=2.0 GPa，泊松比v=0.33，阻尼比ζ=0.05，圓頻率ω=100 rad/s，與孔隙流體相關的參數(shù)取為零。從圖3可以看出，本文方法的結果與文獻[27]結果吻合良好。其次，將TI飽和土層退化為TI彈性土層，與文獻[27]計算結果進行比較。計算中，取土層厚度d=20 m，土層剪切模量Gh=26.7 MPa，Gv=13.3 MPa，基巖半空間剪切模量Gh=Gv=200 MPa，土層阻尼比取ζ=0.05，基巖阻尼比取ζ=0.02，密度均取ρs=ρR=2000 kg?m-3。從圖4可以看出，本文結果與文獻[27]結果一致，驗證了本文方法的正確性。

圖3 各向同性半空間中位移幅值與深度的關系Fig.3 The relationship between displacement amplitude and depth in an isotropic half-space

圖4 層狀TI場地地表位移幅值譜放大曲線Fig.4 Amplification curve of the surface displacement for the surface of a layered TI site

為了驗證本文時域結果的正確性，將退化結果與EERA程序計算結果進行比較。以基巖半空間上單層各向同性場地為計算模型，計算參數(shù)取土層厚度d=10 m，土層剪切模量Gh=Gv=20 MPa，基巖半空間剪切模量Gh=Gv=200 MPa；密度均取ρs=ρR=2000 kg?m-3，阻尼比均取ζ=0.02。圖5給出了 ElCentro波(0.1g)作用時場地地表的水平加速度時程曲線。從圖5可以看出，結果吻合良好。其次，將TI飽和土層退化為TI彈性土層，與文獻[27]計算結果進行比較。計算模型及計算參數(shù)取值為土層厚度d=100 m，土層剪切模量Gh=26.7 MPa，Gv=13.3 MPa，基巖半空間剪切模量Gh=Gv=200 MPa，土層阻尼比取ζ=0.05，基巖阻尼比取ζ=0.02，密度均取ρs=ρR=2000 kg?m-3。圖6 給出了 El Centro 波(0.1g)作用時場地地表的加速度時程曲線。對比表明，兩者結果一致，再次驗證了本文方法的正確性。

圖5 層狀各向同性場地地表加速度時程曲線Fig.5 Acceleration time-history curve for the surface of a layered isotropic site

圖6 層狀TI場地地表加速度時程曲線Fig.6 Acceleration time-history curve for the surface of a layered TI site

3 算例和分析

算例模型為基巖上單一 TI飽和土層場地，研究了土體TI參數(shù)及飽和特性對場地自由場的影響，以歸一化峰值為0.3g的CNTEWGXE地震波(2007年 8月 16日四川珙縣地震長寧臺記錄，I0類基巖波)[28]作為輸入地震動作用于露頭處，獲得了不同工況下場地地表的加速度時程曲線及其反應譜等動力響應結果。文中選取了3種具有不同模量比的TI飽和介質模型進行計算，具體參數(shù)見表1，其中，材料1為各向同性飽和介質。三種場地模型中，土層厚度d=60 m，阻尼比ζ=0.05，孔隙率?=0.3?；鶐r半空間考慮為各向同性材料，剪切模量G=185 MPa，彈性模量Eh=Ev=462.5 MPa，阻尼比ζ=0.02，固體顆粒體積模量ks=1.8 GPa，其余參數(shù)與TI飽和土層參數(shù)相同。

圖7和圖8分別給出了qP1和qSV波以不同角度(θ=0°,30°,60°)入射時，層狀 TI飽和場地地表的水平和豎向加速度時程曲線。從圖中可以看出：

1)TI飽和場地與各向同性飽和場地地表的時程曲線存在一定差異?？傮w上，TI飽和場地(Gh/Gv=2.0和 3.0)的加速度峰值小于各向同性飽和場地(Gh/Gv=1.0)的加速度峰值。對于qP1波，TI飽和場地和各向同性飽和場地差異不大且波形相近；但對于qSV波，隨著場地剪切模量比的增加，波的到時逐漸向右移動。這是由于土體的 TI性質改變了場地的動力特性及波在半空間中的傳播速度，進而改變了加速度峰值及峰值時刻。

2)地震波垂直入射與斜入射時場地地表的時程曲線也有所不同，且主要體現(xiàn)在時程峰值方面。對于qP1波，其垂直入射下的豎向加速度峰值明顯大于其斜入射下的峰值；對于qSV波，隨著入射角的增加，峰值逐漸減小且土體 TI性質的影響逐漸減小。這主要是因為qP1和qSV波斜入射時，改變了波在水平和豎直兩個方向的輸入分量，同時，在土層交界面及地表處存在波型轉換，進而改變了不同類型波的相互作用機制。

表1 TI飽和場地的材料參數(shù)Table 1 Material parameters for TI saturated site

圖7 qP1波入射下場地地表豎向加速度時程曲線Fig.7 Time-history curves of vertical acceleration for the surface under qP1 wave incidence

圖9和圖10進一步給出了相應的場地地表加速度反應譜。從圖中可以看出：對于qP1波，TI飽和場地(Gh/Gv=2.0和 3.0)的反應譜幅值大于各向同性飽和場地(Gh/Gv=1.0)的反應譜幅值，且其在垂直入射下的幅值明顯大于其斜入射下的幅值。對于 qSV波，隨著Gh/Gv比值的增大，反應譜幅值逐漸減??；隨著入射角度的增加，反應譜走勢變得不再規(guī)律且其峰值變化不如qP1波顯著。這是由于Gh/Gv比值的變化改變了場地的動力特性，不同TI參數(shù)場地對地震波產生了不同的濾波和放大效應，進而導致地表反應譜幅值及走勢的改變。

圖8 qSV波入射下場地地表水平加速度時程曲線Fig.8 Time-history curves of horizontal acceleration for the surface under qSV wave incidence

為研究場地飽和特性對場地地震響應的影響，給出了CNTEWGXE波(0.3g)作用在基巖露頭處時，不同工況下場地地表的加速度時程曲線及其反應譜等動力響應結果。土層參數(shù)除滲透率k外取表1中材料1參數(shù)，基巖參數(shù)與上述算例取值一致。其中k1=k3=1000 m2表示不考慮兩方向固液耦合作用(即粘滯性)，k1=k3=10-11m2表示考慮兩方向固液耦合作用。

圖9 qP1波入射下場地地表豎向加速度反應譜Fig.9 Response spectrum of vertical acceleration for the surface under qP1 wave incidence

圖11和圖12分別給出了qP1和qSV波以不同角度(θ=0°,30°,60°)入射時，層狀 TI飽和場地地表的豎向和水平加速度時程曲線。從圖中可以看出：對于qP1波，粘滯性對地表動力響應的影響很小，但在其垂直入射時的影響大于其斜入射時的影響；對于qSV波，粘滯性對地表動力響應影響顯著，表現(xiàn)為考慮固液耦合作用時，地表動力響應明顯減弱。這主要是因為波在飽和多孔介質中傳播時，固液耦合作用所表現(xiàn)出的飽和多孔介質粘滯性會造成剪切波速的衰減及削弱[8]。

圖10 qSV波入射下場地地表水平加速度反應譜Fig.10 Response spectrum of horizontal acceleration for the surface under qSV wave incidence

圖13和圖14進一步給出了相應的場地地表加速度反應譜。從圖中可以看出，飽和多孔介質的固液耦合作用對于入射下的層狀場地的動力響應具有削弱作用但未改變其峰值周期，且隨著入射角度的逐漸增加，這種削弱作用逐漸減小。

圖11 qP1波入射下場地地表豎向加速度時程曲線Fig.11 Time-history curves of vertical acceleration for the surface under qP1 wave incidence

圖12 qSV波入射下場地地表水平加速度時程曲線Fig.12 Time-history curves of horizontal acceleration for the surface under qSV wave incidence

圖13 qP1波入射下場地地表豎向加速度反應譜Fig.13 Response spectrum of vertical acceleration for the surface under qP1 wave incidence

4 結論

本文采用傳遞矩陣方法求解了層狀 TI飽和場地的地震動力響應，通過與已有文獻進行對比驗證了方法的正確性，進而進行了數(shù)值計算分析，研究了土體 TI參數(shù)及飽和特性對場地自由場響應的影響，得到以下主要結論：

(1)TI飽和場地與各向同性飽和場地地表的動力響應存在一定差異。總體上，TI飽和場地的加速度峰值小于各向同性飽和場地的加速度峰值，這一特性在qSV波入射下更為明顯。隨著Gh/Gv比值的增加，場地地表的反應譜峰值呈現(xiàn)減小的趨勢，且場地的TI性質對qSV波的濾波作用更為明顯。

(2)TI飽和場地的飽和特性對場地地表的動力響應存在重要影響且飽和特性對場地地表動力響應的影響主要體現(xiàn)在峰值方面。飽和多孔介質的固液耦合作用對波具有顯著的削弱作用，這一特點在qSV波垂直入射時更為明顯。

(3)地震波垂直入射與斜入射時場地地表的動力響應也有所不同。對于qP1波，其在垂直入射下的反應譜峰值大于其斜入射下的峰值；對于 qSV波，其在垂直入射和其斜入射下峰值變化并不顯著但，但峰值周期有所變化。

附錄A：

其中：

附錄B：

其中：

附錄C：

附錄D：

附錄E：

其中：

附錄F：

其中：