王帥，盛謙，朱澤奇，周春梅

（1.長江科學院水利部巖土力學與工程重點實驗室，武漢 430010;2.中國科學院武漢巖土力學研究所巖土力學與工程國家重點實驗室，武漢 430071;3.武漢工程大學環(huán)境與城市建設(shè)學院，武漢 430073）

頻域分析方法在層狀節(jié)理巖體波動問題中的應(yīng)用

王帥1，2，盛謙2，朱澤奇2，周春梅3

（1.長江科學院水利部巖土力學與工程重點實驗室，武漢 430010;2.中國科學院武漢巖土力學研究所巖土力學與工程國家重點實驗室，武漢 430071;3.武漢工程大學環(huán)境與城市建設(shè)學院，武漢 430073）

考慮波的多次反射，采用頻域分析方法研究一維波在層狀節(jié)理巖體中的傳播問題。首先基于波動理論、線性位移不連續(xù)模型在頻域內(nèi)建立一維波沿多條平行節(jié)理傳播的總體傳遞矩陣，結(jié)合邊界條件解得層狀節(jié)理巖體中一維波傳播的頻域解，然后采用離散傅里葉變換和反變換將頻域解過渡到時域解，最后運用上述方法研究了SV波在不同間距、不同條數(shù)節(jié)理組處的傳播，并與離散元程序UDEC模擬結(jié)果進行對比。結(jié)果表明:隨節(jié)理間距增大，位移透射系數(shù)|TN|（透射波與入射波振幅比值）的變化趨勢呈現(xiàn)上升、下降、平穩(wěn)3個階段;在透射系數(shù)|TN|呈上升階段的節(jié)理區(qū)間內(nèi)，透射系數(shù)|TN|對節(jié)理條數(shù)依耐性很小;通過與UDEC計算結(jié)果對比表明本文求解方法是可行的。

多條節(jié)理;應(yīng)力波;離散傅里葉變換;頻域;線性位移不連續(xù)模型

1 研究背景

自然界中巖體通常是由成組節(jié)理切割形成。研究成組節(jié)理對應(yīng)力波在巖體中的衰減規(guī)律的影響，是了解巖體物理和動態(tài)行為必不可少的步驟，對地質(zhì)勘探、樁基測試以及爆破工程都有重要的意義。

目前應(yīng)力波在節(jié)理處傳播的解析方法主要采用線性位移不連續(xù)模型和等效連續(xù)模型。線性位移不連續(xù)模型，即指節(jié)理面位移不連續(xù)，應(yīng)力連續(xù)。Milie［1］采用該線性不連續(xù)模型，并對節(jié)理非線性本構(gòu)模型進行等效線性化處理，研究了單節(jié)理的滑動摩擦效應(yīng)。Schoenberg［2］得出任意角度入射到單節(jié)理的平面波的透反射系數(shù)。Pyrak-Nolte等［3-4］將其推廣到飽和充水節(jié)理中（開爾文模型和麥克斯韋模型），隨后在不考慮應(yīng)力波的多次反射條件下研究多層節(jié)理中應(yīng)力波傳播問題［5］。石崇［6］采用線性不連續(xù)模型研究了節(jié)理面上的波場及能量分解，建立了節(jié)理面隔振模型。Cai J G和Zhao J［7］，Zhao J和Zhao X B［8］將特征函數(shù)方法引入線性位移不連續(xù)模型，研究了一維彈性波在多條平行節(jié)理處的傳播。在此基礎(chǔ)上，趙堅［9］、Zhao Xiaobo［10］分別研究了節(jié)理法向、切向非線性對應(yīng)力波傳播的影響，發(fā)現(xiàn)節(jié)理的非線性變形行為會引起一種高諧波現(xiàn)象。該方法是在時域內(nèi)求解差分方程獲得層狀節(jié)理巖體動力響應(yīng)，計算時歩大小影響結(jié)果精確性。另外該方法很難考慮斜入射的情形。還有一種求解方法是等效連續(xù)介質(zhì)的方法，采用等效參數(shù)代表離散介質(zhì)整體力學特性［11-12］，該方法沒有考慮應(yīng)力波在節(jié)理之間的多次反射，不適用于節(jié)理尺寸、間距相對較大時的情況。最近，李建春［13-14］提出在等效黏彈性的模型中增設(shè)虛擬波源描述應(yīng)力波在平行節(jié)理之間的多次反射，但需要根據(jù)已有解析解確定黏彈性介質(zhì)模型中的未知參數(shù)。

本文擬考慮波在節(jié)理面的多次反射與透射，借鑒經(jīng)典的水平土層動力響應(yīng)問題中頻域分析方法求解應(yīng)力波在層狀節(jié)理巖體傳播解答。首先根據(jù)波動理論在頻域內(nèi)建立層狀節(jié)理巖體中每個巖層頂部和底部位移、應(yīng)力傳遞矩陣，基于節(jié)理線性不連續(xù)模型建立每個節(jié)理面上下表面的位移、應(yīng)力傳遞矩陣;進而得到應(yīng)力波穿過多條平行節(jié)理傳播的總體傳遞矩陣，聯(lián)合邊界條件解得層狀節(jié)理巖體中應(yīng)力波傳播的頻域解;隨后通過離散傅里葉變換和反變換將頻域解過渡到時域解。最后運用上述方法研究了SV波在不同間距的節(jié)理組處的傳播，與離散元程序UDEC模擬結(jié)果及現(xiàn)有研究成果對比，表明本文頻域-時域求解方法是可行的。

2 層狀節(jié)理巖體中一維波傳播頻域解

利用頻域分析方法討論分層介質(zhì)內(nèi)一維波動，在每層介質(zhì)中因上下界面的反射與透射出現(xiàn)的多個上行和下行諧波都可合成為一個上行諧波和下行諧波［15］。由文獻［2］知，線彈性節(jié)理對應(yīng)力波的頻率沒有影響，只是改變應(yīng)力波波幅、相位，因此上述頻域分析方法同樣適用于節(jié)理巖體中。

層狀節(jié)理巖體模型如圖1所示。巖體內(nèi)含n組節(jié)理，巖層共有n-1個水平層構(gòu)成，層序編號由下到上為2，3，…，n;上部介質(zhì)為第n+1層，下部介質(zhì)為第1層。采用局部坐標系并將z軸的坐標原點設(shè)置在各層底面，正方向垂直向上。設(shè)有一個SV波I自n+1層垂直入射到巖層頂面，此時，在第n+1層中有一反射波R，在第1層中有透射波T。

設(shè)在n+1層巖層入射橫波位移（省略時間因子exp（iωt），下同）為

圖1 一維波通過多條平行節(jié)理示意圖Fig.1Schematic of one-dimensional wave through parallel joints

式中:μ為巖層剪切模量;k為剪切波的圓波數(shù)，k= ω/c，ω為入射波頻率為剪切波的波速; In和Rn分別為第n層巖層透射波和反射波波幅系數(shù)，n=2，…，n+1;T1為第一層巖層透射波波幅系數(shù)。

2.1 層間傳遞矩陣

設(shè)第n層巖層厚度為h，取式（2）、式（3）中z= h，計算第n層巖層頂部的位移分量和應(yīng)力分量，記為

另一方面，取式（2）、式（3）位移和應(yīng)力分量在z=0處的值，可得第n層巖層底部的位移和應(yīng)力值

式中h為巖層厚度。

在式（9）中建立了第n層巖層頂部和底部位移、應(yīng)力分量之間的關(guān)系。

2.2 節(jié)理面?zhèn)鬟f矩陣

第n根節(jié)理上表面位移、應(yīng)力，即第n+1層巖層底部位移、應(yīng)力分量，記為，節(jié)理下表面位移、應(yīng)力，即第n層巖層頂部位移、應(yīng)力分量，記為，根據(jù)節(jié)理線性不連續(xù)模型和本構(gòu)關(guān)系有:

在式（12）中建立了節(jié)理上下面位移分量和應(yīng)力分量之間的關(guān)系。

2.3 總體傳遞矩陣及邊值問題解答

由式（9）、式（12）得到層狀節(jié)理巖體的位移分量和應(yīng)力分量的總體傳遞矩陣

式中Y（AY）n-1為位移分量和應(yīng)力分量的總體傳遞矩陣為第n+1層巖層底部位移、應(yīng)力分量為第1層巖層頂部的位移和應(yīng)力分量。

將式（15）、式（16）代入式（14）解得反射波波幅Rn+1、透射射波波幅T1。

定義透、反射系數(shù)分別為透射波、反射波與入射波波幅比值，可得:

3 層狀節(jié)理巖體的一維波傳播時域解

采用離散傅里葉變換和反變換將頻域解過渡到時域解，具體的工作包括以下2個方面:①利用離散傅里葉變換將暫態(tài)輸入地震波轉(zhuǎn)換為穩(wěn)態(tài)輸入波;②在頻域內(nèi)求解計算層狀節(jié)理巖體穩(wěn)態(tài)反應(yīng)，再用快速傅里葉反變換得到暫態(tài)反應(yīng)。

3.1 離散傅里葉變換和輸入波動的分解

為了將暫態(tài)輸入波u'I（t）轉(zhuǎn)換為穩(wěn)態(tài)輸入波，可以采用離散傅里葉變換。將暫態(tài)時間函數(shù)u'I（t）用周期為T的時間函數(shù)uI（t）替換（如圖2所示）［15］，即

式中周期T不僅大于u'I（t）的持續(xù)時間，而且大于層狀節(jié)理巖體暫態(tài)反應(yīng)的持續(xù)時間。

圖2 入射波的周期化Fig.2Cycle of the incident wave

將單個周期內(nèi)的uI（t）以時間歩距Δt離散化得

uIj展開成如下傅里葉級數(shù)形式，

uIj的角頻率為ωm（m=0，1，…，J-1），傅里葉譜CIm為

式中ωm=m·2πΔf，Δf=1/T。

uIj與CIm構(gòu)成一一對應(yīng)的離散傅里葉變換對，可用快速傅里葉算法FFT計算。

3.2 層狀節(jié)理巖體的暫態(tài)波動反應(yīng)

由于離散傅里葉變換計算時間長度大于輸入波uI（t）及系統(tǒng)反應(yīng)量的持續(xù)時間，因此系統(tǒng)反應(yīng)量計算亦可用傅里葉變換對式（22）、式（23）計算，則第1層巖層頂部位移u1（t）可表示為

式中:u1j=u1（jΔt）;C1m為u1（t）的離散傅里葉譜。

利用上述離散傅里葉變換對，暫態(tài)反應(yīng)歸結(jié)為計算離散傅里葉譜C1m。由式（18）知第1層巖層中頂面上透射波的離散傅里葉譜C1m為

4 參數(shù)研究

本文主要研究節(jié)理間距h、節(jié)理條數(shù)N對應(yīng)力波在層狀節(jié)理巖體傳播的影響。節(jié)理剛度Ks為2 GPa/m，巖石密度ρ為2 650 kg/m3，彈性模量E為26.5 GPa，泊松比v為0.25，SV波在巖石中傳播速度c為2 000 m/s。假設(shè)入射SV波為含半個循環(huán)的正旋波，頻率50 Hz，振幅為1 mm，波長λ為40 m。將暫態(tài)輸入波轉(zhuǎn)換成周期T=0.512 s的穩(wěn)態(tài)波，后對單個周期內(nèi)穩(wěn)態(tài)波進行離散傅里葉變換，歩距Δt=0.512/1 024=5×10-4s，離散數(shù)據(jù)點1 024個，得到穩(wěn)態(tài)波動一個周期內(nèi)時程及其傅里葉譜如圖3所示。

當節(jié)理條數(shù)N=2和4時，對于不同的節(jié)理間距h，基于本文解析方法和離散元程序UDEC得到的透射波如圖4所示。由圖看出基于2種方法得到的結(jié)果相同，從而驗證了本文方法的正確性。同時可知當節(jié)理間距h為波長λ時，第一次到達透射波和第二次到達透射波有明顯的時間間隔，第一次到達透射波是初始入射應(yīng)力波的直接傳播的結(jié)果，第二次到達的透射波來自波在節(jié)理之間的多次透射和反射;當節(jié)理間距h為波長λ的1/10時，第一次到達透射波則疊加了在節(jié)理面上多次反射與透射的波，波幅發(fā)生一定改變。

圖3 單個周期內(nèi)穩(wěn)態(tài)輸入波動時程及其傅里葉譜Fig.3Steady-state incident wave within a single cycle and the corresponding Fourier amplitude spectrum

圖4 解析方法與離散元程序UDEC得到的透射波時程Fig.4Transmitted waves obtained from the analytic solution and the discrete element program UDEC

圖5為入射波穿過N條節(jié)理的透射系數(shù)|TN| （透射波位移振幅與入射波位移振幅的比值）與正交化節(jié)理間距ζ（h/λ）關(guān)系曲線。可知:

（1）隨節(jié)理間距增大，|TN|的變化趨勢呈現(xiàn)上升、下降、平穩(wěn)3個階段。

（2）在|TN|呈現(xiàn)上升階段的節(jié)理區(qū)間內(nèi)，透射系數(shù)|TN|對節(jié)理條數(shù)依耐性很小，其他區(qū)間隨節(jié)理條數(shù)增大減小。

圖5 透射系數(shù)|TN|與正交化節(jié)理間距ζ的關(guān)系曲線Fig.5Curves of transmission coefficient|TN|vs. orthogonalized joint spacing ζ

5 結(jié)論

本文考慮波在節(jié)理面的多次反射與透射，借鑒經(jīng)典的水平土層動力響應(yīng)問題中頻域分析方法建立應(yīng)力波在層狀節(jié)理巖體傳播解答。即首先基于波動理論、節(jié)理線性位移不連續(xù)模型，建立一維波穿過多條平行節(jié)理傳播的總體傳遞矩陣，結(jié)合邊界條件獲得一維波傳播的頻域解，通過離散傅里葉變換和反變換將頻域解過渡到時域解;隨后研究了SV波在不同間距、不同條數(shù)節(jié)理組處的傳播規(guī)律。結(jié)果表明:隨節(jié)理間距增大，位移透射系數(shù)|TN|（透射波與入射波振幅比值）的變化趨勢呈現(xiàn)上升、下降、平穩(wěn)3個變化階段，在透射系數(shù)|TN|呈上升階段的節(jié)理區(qū)間內(nèi)，透射系數(shù)|TN|對節(jié)理條數(shù)依耐性很小。

與離散元程序UDEC模擬結(jié)果對比表明，本文頻域-時域求解方法是可行的，而且相對于特征值方法，本文方法能得出透反射系數(shù)解析公式，便于參數(shù)研究，同時可進一步推廣到應(yīng)力波傾斜入射的情形，相對于增設(shè)虛擬波源的等效方法，本文方法計算中輸入?yún)?shù)不依賴于應(yīng)力波屬性。

［1］MILLER R K.Effects of Boundary Friction on Propagation of Elastic-Waves［J］.Bulletin of the Seismological Society of America，1978，68（4）:987-998.

［2］SCHOENBERG M.Elastic Wave Behavior Across Linear Slip Interfaces［J］.Journal of the Acoustical Society of America，1980，68（5）:1516-1521.

［3］PYRAK L J.Seismic Visibility of Fractures［D］.California，US:University of California，Berkeley，1988.

［4］PYRAK-NOLTE L J，MYER L R，COOK N G W.Transmission of Seismic Waves Across Single Natural Fractures［J］.Journal of Geophysical Research，1990，95（B6）: 8617-8638.

［5］PYRAK-NOLTE L J，MYER L R，COOK N G W.Anisotropy in Seismic Velocities and Amplitudes From Multi-ple Parallel Fractures［J］.Journal of Geophysical Research，1990，95（B7）:11345-11358.

［6］石崇，徐衛(wèi)亞，周家文，等.節(jié)理面透射模型及其隔振性能研究［J］.巖土力學，2009，30（3）:729-734. （SHI Chong，XU Wei-ya，ZHOU Jia-wen，et al.Transmission Model of Joint Interface and Its Performance of Vibration Isolation［J］.Rock and Soil Mechanics，2009，30（3）:729-734.（in Chinese））

［7］CAI J G，ZHAO J.Effects of Multiple Parallel Fractures on Apparent Attenuation of Stress Waves in Rock Masses［J］.International Journal of Rock Mechanics and Mining Sciences，2000，37（4）:661-682.

［8］ZHAO J，ZHAO X B，CAI J G.A Further Study of PWave Attenuation across Parallel Fractures with Linear Deformational Behaviour［J］.International Journal of Rock Mechanics＆Mining Sciences，2006，43:776-788.

［9］趙堅，蔡軍剛，趙曉豹，等.彈性縱波在具有非線性法向變形本構(gòu)關(guān)系的節(jié)理處的傳播特征［J］.巖石力學與工程學報，2003，22（1）:9-17.（ZHAO Jian，CAI Jun-gang，ZHAO Xiao-bao，et al.Transmission of Elastic P-Waves across Single Fracture with Nonlinear Normal Deformation Behavior［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（1）:9-17. （in Chinese））

［10］ZHAO J，ZHAO X B，HEFNY A M，et al.Normal Transmission of S-Wave across Parallel Fractures with Coulomb Slip Behavior［J］.Journal of Engineering Mechanics，2006，132（6）:641-650.

［11］SCHOENBERG M，MUIR F.A Calculus for Finely Layered Anisotropic Media［J］.Geophysics，1989，54（5）: 581-589.

［12］HUDSON J A.Wave Speeds and Attenuation of Elastic Waves in Material Containing Cracks［J］.Geophysical Journal of the Royal Astronomical Society，1981，64（1）: 133-150.

［13］LI J C，MA G W，ZHAO J.An Equivalent Viscoelastic Model for Rock Mass with Parallel Joints［J］.Journal of Geophysical Research:Solid Earth（1978-2012），2010，115（B3）:B3305.

［14］李建春，李海波，MA G W，等.節(jié)理巖體的一維動態(tài)等效連續(xù)介質(zhì)模型的研究［J］.巖石力學與工程學報，2010，29（2）:4063-4067.（LI Jian-chun，LI Hai-bo，MA G W，et al.An Equivalent 1D Dynamic Continuum Model for Rock Mass［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（2）:4063-4067. （in Chinese））

［15］廖振鵬.工程波動理論導論［M］.北京:科學出版社，2002.（LIAO Zhen-peng.Introduction to Theory of Engineering Waves［M］.Beijing:Science Press，2002.（in Chinese ））

（編輯:王慰）

Frequency Domain Analysis Method for Solving the Stress Wave Propagation in Layered Jointed Rock Mass

WANG Shuai1，2，SHENG Qian2，ZHU Ze-qi2，ZHOU Chun-mei3
（1.Key Laboratory of Geotechnical Mechanics and Engineering of MWR，Yangtze River Scientific Research Institute，Wuhan430010，China;2.State Key Laboratory of Geomechanics and Geotechnical Engineering，Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，Wuhan430071，China; 3.School of Environment and Civil Construction，Wuhan Institute of Technology，Wuhan430073，China）

In consideration of multiple wave reflections，frequency domain analysis method was adopted to solve one-dimensional wave propagation in layered jointed rock.Firstly，the overall transfer matrix of one-dimensional wave propagation along multiple parallel joints was constructed by using wave theory and linear displacement discontinuity model.Solution for one-dimensional wave propagation in frequency domain can be obtained with boundary conditions.Subsequently，discrete Fourier transform and inverse discrete Fourier transform were used to transform the solution from frequency domain to time domain.Finally，SV wave propagation through joints of different spaces and numbers was studied by the above methods，and the simulation result was compared with that from discrete element program UDEC.The results showed that the transmission coefficient|TN|（ratio of transmitted wave amplitude to incident wave amplitude）rose first，then declined，and finally became stable with the increasing of joint spacing.During the rising stage of，it has less dependence on the number of joints.Through comparison with UDEC results，the method in this paper was found to be feasible.

multiple joints;stress wave;discrete Fourier transform;frequency domain;discontinuity model of linear displacement

TV321

A

1001-5485（2013）04-0062-05

10.3969/j.issn.1001-5485.2013.04.014

2013，30（04）:62-66

2012-02-13;

2012-04-09

國家自然科學基金重大研究計劃資助項目（90715042）;國家自然科學基金項目（51009130）;國家科技支撐計劃課題子題（2006BAB04A06）

王帥（1986-），男，安徽六安人，工程師，博士，主要從事大型地下洞室群的地震穩(wěn)定研究工作，（電話）15102766897（電子信箱） wsswjtu@qq.com。