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

(1.中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點(diǎn)實(shí)驗(yàn)室(天津大學(xué))，天津 300350；2.濱海土木工程結(jié)構(gòu)與安全教育部重點(diǎn)實(shí)驗(yàn)室，天津300350；3.天津大學(xué)土木工程系，天津 300350)

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

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

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

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

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

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

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

1.1 TI飽和多孔介質(zhì)波動(dòng)方程及其基本解

根據(jù)文獻(xiàn)[19]，TI飽和多孔介質(zhì)平面內(nèi)應(yīng)力-應(yīng)變關(guān)系為：

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

直角坐標(biāo)系下[22]，以u(píng)-w形式表示的Biot動(dòng)力平衡方程為：

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

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

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

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

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

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

將式(5)代入式(1)中，可得到平面內(nèi)應(yīng)力及孔壓表達(dá)式如下：

1.2 層狀TI飽和半空間平面內(nèi)傳遞矩陣

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

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

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

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

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

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

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

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

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

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

1.3 邊界條件及地表動(dòng)力響應(yīng)

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

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

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

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

式中：當(dāng)邊界透水時(shí)，M1~M9具體表達(dá)式見附錄E；當(dāng)邊界不透水時(shí)；M1~M9具體表達(dá)式見附錄F。

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

2 方法驗(yàn)證

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

圖3 各向同性半空間中位移幅值與深度的關(guān)系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

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

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

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

3 算例和分析

算例模型為基巖上單一 TI飽和土層場地，研究了土體TI參數(shù)及飽和特性對(duì)場地自由場的影響，以歸一化峰值為0.3g的CNTEWGXE地震波(2007年 8月 16日四川珙縣地震長寧臺(tái)記錄，I0類基巖波)[28]作為輸入地震動(dòng)作用于露頭處，獲得了不同工況下場地地表的加速度時(shí)程曲線及其反應(yīng)譜等動(dòng)力響應(yīng)結(jié)果。文中選取了3種具有不同模量比的TI飽和介質(zhì)模型進(jìn)行計(jì)算，具體參數(shù)見表1，其中，材料1為各向同性飽和介質(zhì)。三種場地模型中，土層厚度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°)入射時(shí)，層狀 TI飽和場地地表的水平和豎向加速度時(shí)程曲線。從圖中可以看出：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

附錄A：

其中：

附錄B：

其中：

附錄C：

附錄D：

附錄E：

其中：

附錄F：

其中：