黏彈性節理巖體中應力波的傳播特性分析

汪書敏，王志亮，賈帥龍，王浩然，王昊辰

(合肥工業大學 土木與水利工程學院，合肥 230009)

節理裂隙作為一種普遍存在于巖體中的缺陷結構，對應力波(爆炸波、地震波等)在巖體中的傳播規律產生較大影響。應力波與節理之間的作用關系復雜，常常導致振幅下降、波形變化和能量損失[1-2]。巖體中的天然裂隙常常由飽和砂土、黏土以及風化巖石等具有黏彈性變形特性的介質填充，在自然條件下逐漸演化成黏彈性節理，故研究具有黏彈特性的節理對應力波傳播的影響機制在抗震工程、損傷測定和巖石力學等領域都具有重大意義。

國內外學者針對應力波在節理巖體中的傳播過程進行了大量研究。王志亮等[3]提出了一種三參數模型用于描述節理巖體的動態響應，基于位移不連續理論(DDM)推導了一維應力波的傳播方程并與實驗進行對比，發現當應力波持續時間較短時，理論結果與實驗結果吻合較好。Li等[4-7]根據波前動量守恒原理提出了一種新的時域遞歸方法，相繼導出了應力波通過線彈性節理、非線彈性節理以及考慮節理切向滑移特性時的傳播方程。而對于天然夾層裂隙中常常充填有飽和砂土、黏土等形成的節理，Fehler等[8-9]認為應考慮其黏彈特性，對此，也有不少學者進行了相關研究。Huang等[10]利用矩陣傳播法研究了應力波通過一組黏彈性節理的傳播規律，發現在Kelvin節理中存在一個特定黏度系數值使得能量耗散率最大。Yi等[11]采用等效介質法將節理充填介質等效為靜態的各向同性介質，發現該方法不能準確反映填充介質較疏時的振幅各向異性和波速隨頻率變化的規律。Zhu等[12]考慮填充介質質量產生的慣性應力，提出了位移與應力不連續理論(DSDM)，發現應力波通過填充節理時的能量耗散率由節理黏度系數和填充介質質量決定。Wu等[13]使用分離式霍普金森巖桿對充填節理進行的動態試驗結果表明，DSDM比DDM更適合描述厚度在10 mm以上的填充節理對波傳播的影響。此外，Li等[14-16]提出了一種薄層界面模型用于研究應力波入射充填節理時的傳播規律，結果表明，在剪切波入射該模型引起剪應力大于薄層界面的抗剪強度而產生滑移時該方法仍有效。

綜上，目前基于時域遞歸方法(TDRM)且考慮應力不連續條件的平行黏彈性節理中波傳播特性研究開展較少。因此，本研究先采用TDRM法，引入飽依丁-湯姆遜本構模型描述節理的黏彈特性，并考慮填充節理質量，將DDM擴展至DSDM，導出了應力波在一組平行黏彈性節理中傳播的完整解析方程。通過與頻域內的封閉解對比驗證該方法的正確性。最后,分析了節理模型參數、無量綱節理厚度和平行節理間距等參數對透反射系數和能量耗散率的影響，力求得出有參考價值的結論。

1 理論分析

1.1 問題描述

在離爆炸、地震等擾動源較遠的地層中，應力波相對于節理面可看作平面波。當應力波穿過一組等間距平行節理時，在節理處會發生透射和反射現象，分別形成透射P波與S波，反射P波與S波，如圖1所示，S為節理間距。圖2為應力波通過第X條節理的左右波形圖。其中，d為節理厚度，符號“-”和“+”分別代表節理的左側與右側，LP-、LS-、RP-和RS-分別為節理左側的左行P波、左行S波、右行P波和右行S波，LP+、LS+、RP+和RS+分別為節理右側的左行P波、左行S波、右行P波和右行S波。假設節理兩側的巖石均為線彈性介質，且密度與泊松比相同，P波入射角為α，S波入射角為β。應力波通過含黏彈特性的充填介質時，會產生一定的邊際彌散效應而損失部分到達節理右側的能量。本研究不考慮應力波入射充填節理時產生的散射效應，認為該透反射過程遵循Snell定律。故右側出射波角度與左側入射角相同，即透反射P波和S波與法線夾角均分別為α和β[12,15,17]。

圖1 應力波通過平行節理Fig.1 Stress wave traversing through parallel joint

圖2 第X條節理兩側的左行波和右行波Fig.2 Left-and right-running waves on both sides of the Xth joint

1.2 理論公式推導

(1)

(2)

根據Snell定律

(3)

式中：cp和cs分別為P波和S波在巖石中的波速。圖3(b)～3(h)分析過程同圖3(a)。

圖3 應力波波前與節理兩側應力Fig.3 Stress on wavefront and rock joint surfaces

結合上述結果及波前動量守恒原理，根據Li等[4]推導過程可得節理兩側的法向應力、切向應力、法向速度與切向速度：

(4)

(5)

(6)

(7)

式中：上標e分別代表“-”或“+”。下標r和l分別代表右行波和左行波。

本研究引入飽依丁-湯姆遜本構模型對巖體節理的黏彈特性進行描述，如圖4所示。由一個Maxwell體與彈簧元件并聯而成，本構方程為

圖4 飽依丁-湯姆遜本構模型[18]Fig.4 Poyting-Thomson constitutive model[18]

(8)

節理左右界面的連接如圖2所示，考慮節理填充介質質量的應力不連續條件[12,19]：

(9)

(10)

將式(10)作為連接節理左右兩側法向和切向的位移不連續條件可以得

(11)

將式(4)～(7)代入式(9)和(11)，并整理為矩陣形式：

(12)

(13)

式中的相關系數矩陣如下：

假定節理間的巖石為完全彈性，則應力波由第X-1條節理右側出射到達第X條節理左側或由第X+1條節理左側反射至第X條節理右側的過程中，只相差一個時間差，故時移函數方程表達式為

(14)

(15)

式中：下標X表示第X條節理,S表示相鄰節理間的間距。在計算過程中，時間步長為Δt，則兩側速度時移步數差為np=S/(cosα·cp·Δt)或ns=S/(cosβ·cs·Δt)。

式(12)～(15)即為應力波在一組平行節理中傳播的完整方程。通過編制MATLAB程序，輸入相應已知條件和初始條件并采用滿足精度要求的時間步長，對上述方程進行時域遞歸計算即可得解。

2 理論驗證

為方便研究節理對應力波在巖體中傳播特征的影響，定義反射和透射系數分別為

(16)

(17)

式中：下標a和b分別表示P波或S波，下標N表示平行節理條數。

引入兼具Maxwell體與Kelvin體黏彈特性的飽依丁-湯姆遜黏彈本構模型作為位移不連續的變形條件。當不考慮節理的黏性特征和應力不連續條件時，可將黏性系數η及慣性質量m均設為0，便退化為線彈性節理，容易發現退化后的方程與Li等[4]一致。當彈簧元件k1趨于正無窮時，原模型演變為Kelvin模型，適合描述飽和土填充節理的動態響應規律[21-22]。而將彈簧元件k2設為0時，原模型便退化為適合描述橫波穿過黏土填充的薄層節理情形的Maxwell模型[23]。為驗證上述所推導方程的正確性，將計算結果與Zhu等[12]所得的結果進行對比。假定在完整巖石中P波和S波的傳播速度分別為6 131和3 830 m/s，節理填充介質中cplate=1 000 m/s，這些參數均可通過巖石和節理填充介質的楊氏模量、密度和泊松比求得[12,19]。根據Snell定律可得入射P波和入射S波的入射臨界角分別為90°與38.65°，分析中僅討論入射角小于臨界角的情況。定義兩個法向特定剛度系數分別為Kn1=kn1/(zpω)與Kn2=kn2/(zpω)，特定黏度系數Hn=ηn/zp，兩個切向特定剛度系數分別為Kτ1=kτ1/(zsω)與Kτ2=kτ2/(zsω)，特定黏度系數為Hτ=ητ/zs。ω為入射波的圓頻率，即ω=2πf，zp和zs分別為P波和S波的波阻抗。不同入射角下的透反射系數計算結果如圖5、6所示，可以看出，Kelvin節理和Maxwell節理隨入射角的變化趨勢大致相同，臨界角之前透射系數變化較為平穩，臨界角附近急劇變化，這與Zhu等[12]所得頻域內封閉解具有良好的一致性。

3 參數討論與分析

節理模型參數往往通過霍普金森壓桿(SHPB)對充填介質進行動態沖擊實驗的方式獲取[12]。但由于SHPB實驗中S波難以獲得，切向參數的確定存在困難。現假設節理法向與切向同性，即取Kτ1=Kn1，Kτ2=Kn2，Hτ=Hn，除待研究參數外，其余參數同上。設入射波為半個周期的正弦P波或S波，如式(18)所示，Is=Ip。

(18)

式中f為正弦波的頻率，取1 kHz。

定義法向入射波的能量耗散率為[12]

eloss=1-(R2+T2)

(19)

3.1 節理模型參數的影響分析

為探究模型參數對透反射系數與能量耗散率的影響，以P波通過單條節理為例，特定剛度系數Kn2分別取為0.2、0.4、0.8和1.6。透反射系數及能量耗散率隨1/Kn1和1/Hn的變化曲線分別如圖7、8所示。可以看出，當其他參數為定值時，Kn2對節理模型的影響較為明顯，透射系數隨著Kn2的增加而增大，而反射系數與能量耗散率均減小。

由圖7可知，當Kn2較小時，隨著1/Kn1的增加，反射系數逐漸增加，而透射系數初期有一小段上升，隨后緩慢下降；當Kn2較大時，反射和透射系數的初期上升趨勢趨弱，且后續隨1/Kn1的增大基本持平。能量耗散率隨著1/Kn1的增大而逐漸下降，并趨于一定值e0。由于Kn2的增加，模型中Maxwell體的影響逐漸弱化，節理的變形特征由黏彈性逐漸向線彈性演變，因而能量耗散率相應下降。故e0主要由Kn2決定，且隨Kn2的增大而減小。

由圖8可知，當Kn2較小時，反射系數初期小幅下降，且下降的趨勢隨Kn2的增加而逐漸趨緩，但總體隨1/Hn增加而上升。同時，透射系數隨1/Hn的增加而逐漸下降，且下降幅度隨Kn2增加而逐漸減小。能量耗散率eloss隨1/Hn的增加逐漸上升，而當1/Hn達2～3之后則基本持平。

3.2 無量綱節理厚度的影響分析

應力波在節理巖體傳播的過程中，節理厚度同樣會對應力波的傳播特性產生影響。考慮到透反射系數的波頻依賴性，定義無量綱節理厚度h為ωd/(2πcplate)。圖9為IS和IP法向入射單條節理時透反射系數及能量耗散率隨h的變化曲線。可以看出，當S波法向入射時，透反射系數及能量耗散率的大小隨lgh的變化趨勢與P波法向入射時基本相同。當h<10-2時，無量綱節理厚度對透反射系數的影響不明顯。隨著h的增大，透射系數出現明顯的下降，并逐漸接近于0。而反射系數隨著無量綱節理厚度的增加呈現先降低后增加的現象。此外，隨著厚度增加，S波入射時透反射系數的變化趨勢要先于P波，而能量耗散率隨無量綱節理厚度的增加表現出先上升后下降的趨勢。

圖9 透反射系數與能量耗散率隨無量綱節理厚度的變化Fig.9 Variation of transmission and reflection and energy dissipation rate with non-dimensional joint thickness

3.3 應力波入射角的影響分析

取入射S波表達式同P波，圖10、11分別為不同入射角下P波和S波在節理間距為1 m的一組平行等間距節理中傳播的透射系數TPP、TPS、TSS和TSP。可以看出，當P波和S波以任意角度入射時，隨著節理條數的增加，TPP和TSS均出現了不同程度的下降。在入射P波與S波分別接近臨界角αc=90°與βc=arcsin(cs/cp)=38.6°時，TPP和TSS分別出現急劇下降或上升的現象。隨著節理條數N的增加，透反射系數在臨界角處的“尖點”顯著下降，而TPS和TSP隨著節理條數的增加而逐漸增大。

圖10 一組平行節理對不同入射角度P波透射系數的影響Fig.10 Influence of a set of parallel joints on the transmission coefficients of incident P-wave at different angles

圖11 一組平行節理對不同入射角度S波透射系數的影響Fig.11 Influence of a set of parallel joints on the transmission coefficients of incident S-wave at different angles

3.4 無量綱節理間距的影響分析

定義無量綱節理間距ξ為節理間距S與入射應力波波長的比值。本節采用半個周期的正弦入射P波，如式(18)所示。圖12為透射系數與無量綱節理間距ξ的關系，易見透射系數隨ξ的變化可以分為上升區間(ξ≤ξcri)、下降區間(ξcri<ξ<ξthr)以及不變區間(ξthr≤ξ)。在節理條數較多時，TPP隨無量綱節理間距的增加呈線性下降，而黏彈節理則是分為兩段——陡峭段和平緩段下降。圖13為臨界值ξcri和ξthr隨節理數N的變化曲線，可以看出,隨著節理數的增加，臨界值ξcri緩慢下降，且基本保持線性關系；對于臨界值ξthr，當N<5時，增長較緩；當N>6時，增長較快。

圖12 透射系數與無量綱節理間距ξ的關系Fig.12 Transmission coefficients versus non-dimensional joint spacing (ξ)

圖13 ξcri和ξthr隨節理數N的變化Fig.13 Variation of ξcri and ξthr with joint number N

4 結 論

1)時域遞歸方法可用于推導任意角度的入射應力波通過黏彈性節理的傳播方程。引入的飽依丁-湯姆遜模型退化為Kelvin模型和Maxwell模型的計算結果不僅驗證了所用模型及推導過程的正確性，同時印證了其兼具兩種模型的特點，為描述節理復雜的黏彈特性提供了新的模型選擇。

2)參數化研究表明，剛度系數、黏度系數、無量綱節理厚度等都影響透反射系數及能量耗散率。其中，特定剛度系數Kn2與透射系數呈正相關，與反射系數及能量耗散率呈負相關。當無量綱節理厚度小于10-2時，其變化對透反射系數及能量耗散率的影響不明顯。

3)透反射系數隨黏彈性節理數及節理間距的變化而變化。對于入射P波(S波)，透射P波(S波)隨節理數增加而逐漸下降，透射S波(P波)隨節理數的增加逐漸上升。當節理條數較少時，節理條數的變化對透射系數影響較大；節理條數較多時，其變化對透射系數的影響則較小。此外，ξcri隨節理數的增長基本保持線性關系，但ξthr隨節理數的增加呈現先緩后陡的增長趨勢。