黏彈性節(jié)理巖體中應(yīng)力波的傳播特性分析

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

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

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

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

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

1 理論分析

1.1 問題描述

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

圖1 應(yīng)力波通過平行節(jié)理Fig.1 Stress wave traversing through parallel joint

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

1.2 理論公式推導(dǎo)

(1)

(2)

根據(jù)Snell定律

(3)

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

圖3 應(yīng)力波波前與節(jié)理兩側(cè)應(yīng)力Fig.3 Stress on wavefront and rock joint surfaces

結(jié)合上述結(jié)果及波前動量守恒原理，根據(jù)Li等[4]推導(dǎo)過程可得節(jié)理兩側(cè)的法向應(yīng)力、切向應(yīng)力、法向速度與切向速度：

(4)

(5)

(6)

(7)

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

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

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

(8)

節(jié)理左右界面的連接如圖2所示，考慮節(jié)理填充介質(zhì)質(zhì)量的應(yīng)力不連續(xù)條件[12,19]：

(9)

(10)

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

(11)

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

(12)

(13)

式中的相關(guān)系數(shù)矩陣如下：

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

(14)

(15)

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

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

2 理論驗證

為方便研究節(jié)理對應(yīng)力波在巖體中傳播特征的影響，定義反射和透射系數(shù)分別為

(16)

(17)

式中：下標(biāo)a和b分別表示P波或S波，下標(biāo)N表示平行節(jié)理條數(shù)。

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

3 參數(shù)討論與分析

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

(18)

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

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

eloss=1-(R2+T2)

(19)

3.1 節(jié)理模型參數(shù)的影響分析

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

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

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

3.2 無量綱節(jié)理厚度的影響分析

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

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

3.3 應(yīng)力波入射角的影響分析

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

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

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

3.4 無量綱節(jié)理間距的影響分析

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

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

圖13 ξcri和ξthr隨節(jié)理數(shù)N的變化Fig.13 Variation of ξcri and ξthr with joint number N

4 結(jié) 論

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

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

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