低幅值應(yīng)力波在巖石三維節(jié)理面的折反射規(guī)律理論研究

劉傳正, 張建經(jīng), 崔 鵬

(1.中國(guó)科學(xué)院 山地災(zāi)害與地表過程重點(diǎn)實(shí)驗(yàn)室 成都山地災(zāi)害與環(huán)境研究所,成都 610041；2.中國(guó)科學(xué)院大學(xué) 工程科學(xué)學(xué)院，北京 100049；3.西南交通大學(xué) 土木工程學(xué)院，成都 610031)

巖體中的結(jié)構(gòu)面很大程度上影響了巖體的力學(xué)性質(zhì)。近些年來，由于地震事故頻發(fā)，并引發(fā)了大量的崩塌、滑坡災(zāi)害，地震作用下巖土體失穩(wěn)機(jī)理成為當(dāng)前十分迫切的研究任務(wù)之一。而巖石節(jié)理對(duì)應(yīng)力波傳播的影響規(guī)律對(duì)于研究地震巖質(zhì)邊坡動(dòng)力響應(yīng)與失穩(wěn)機(jī)理具有重要意義。

巖石節(jié)理力學(xué)性質(zhì)的研究始終貫穿于巖石力學(xué)中的研究中，是巖石工程常常所面對(duì)的不可忽視的客觀存在。巖石節(jié)理力學(xué)研究早期是對(duì)二維粗糙節(jié)理力學(xué)進(jìn)行研究，而對(duì)于三維粗糙節(jié)理的研究則相對(duì)較為有限，但是仍有大量的研究發(fā)現(xiàn)，如Belem等[1]通過對(duì)節(jié)理面幾何統(tǒng)計(jì)分析發(fā)現(xiàn)，對(duì)于自然形成節(jié)理，巖石節(jié)理面凸起的傾角和曲率以及節(jié)理的粗糙度與曲折度，都會(huì)表現(xiàn)出一定的各向異性性質(zhì)。孫輔庭等[2]運(yùn)用分形數(shù)學(xué)理論，提出了一種基于三維均方根抵抗角的節(jié)理面粗糙度分析方法，分形參數(shù)在方位角上表現(xiàn)出明顯的各向異性性質(zhì)。Jing等[3]通過巖石節(jié)理面剪切實(shí)驗(yàn)發(fā)現(xiàn)，在各方向上節(jié)理起伏角并不是均勻分布的，且節(jié)理的剪切變形剛度系數(shù)在各個(gè)方向上也不相等；進(jìn)一步研究發(fā)現(xiàn)巖石節(jié)理起伏角的各向分布與剪切剛度在方向上分布呈正相關(guān)，且具有中心對(duì)稱的兩個(gè)優(yōu)勢(shì)方向，可用橢圓描述。Jing等[4]在巖石節(jié)理表面形態(tài)各向異性研究的基礎(chǔ)上，給出了巖石節(jié)理三維本構(gòu)模型。類似的，王光綸等[5]給出了具有非線性彈性變形的巖石節(jié)理三維本構(gòu)方程。而Misra[6]則基于三維彈塑性接觸理論、節(jié)理面形態(tài)的幾何參數(shù)分布函數(shù)給出了積分形式的巖石節(jié)理三維本構(gòu)力學(xué)模型。以上研究均表明中皆從巖石三維節(jié)理面的形態(tài)作為主要出發(fā)點(diǎn)，有些給出三維節(jié)理本構(gòu)方程，然后都過去研究中可以發(fā)現(xiàn)三維節(jié)理變形本構(gòu)方程的剛度矩陣中交叉耦合項(xiàng)非零；由于節(jié)理變形方程中交叉耦合剛度系數(shù)是描述巖石節(jié)理力學(xué)性質(zhì)重要參數(shù)，如描述了節(jié)理的剪脹、減縮現(xiàn)象，而且節(jié)理本構(gòu)方程中交叉耦合剛度系數(shù)不僅僅描述這兩種現(xiàn)象，其對(duì)應(yīng)力波在節(jié)理上折反射規(guī)律的影響相關(guān)研究也十分缺乏。

相比巖石節(jié)理力學(xué)性質(zhì)的研究，過去在研究巖石節(jié)理應(yīng)力波傳播規(guī)律中對(duì)節(jié)理本構(gòu)模型則做出了較大簡(jiǎn)化，其很大原因是由于問題求解的復(fù)雜性。Schoenberg[7]基于位移非連續(xù)模型，給出應(yīng)力波在線剛度非連續(xù)面上的折反射解析解，隨后很多學(xué)者對(duì)此方法進(jìn)行了研究和補(bǔ)充，如Myer等[8-15]通過理論、實(shí)驗(yàn)對(duì)線剛度DDM(Displacement Discontinue Method)模型進(jìn)行大量研究，證明了該理論的可靠性和實(shí)用性。Zhao等[16]給出了縱波垂直入射非線性變形節(jié)理面的解析解；俞縉等[17-18]對(duì)非線性節(jié)理面的應(yīng)力波傳播模型進(jìn)行了改進(jìn)。Daehnke等[19]給出考慮交叉耦合項(xiàng)的二維節(jié)理面DDM模型解析解。Misra等[20]在給出的三維節(jié)理本構(gòu)模型的基礎(chǔ)上，分析了應(yīng)力波的折反射行為，但是其理論框架仍是在傳統(tǒng)的認(rèn)識(shí)范圍內(nèi)。綜上述，以往研究沒有對(duì)三維粗糙節(jié)理的應(yīng)力波折反射行為進(jìn)行探討，沒有全面的分析三維節(jié)理本構(gòu)方程中的交叉耦合剛度系數(shù)項(xiàng)對(duì)三種簡(jiǎn)諧波傳播的影響。

本文基于巖石節(jié)理的細(xì)觀接觸模型，討論了巖石節(jié)理一般形式的三維應(yīng)力-應(yīng)變本構(gòu)方程，并以此為基礎(chǔ)分析了簡(jiǎn)諧波在節(jié)理面上的相互轉(zhuǎn)化機(jī)理，并據(jù)此推導(dǎo)了三維節(jié)理上應(yīng)力波折反射的理論解，最后重點(diǎn)討論了節(jié)理本構(gòu)方程中的交叉耦合項(xiàng)對(duì)應(yīng)力波折反射系數(shù)的影響。

1 巖石節(jié)理的形變剛度討論

1.1 三維粗糙節(jié)理接觸變形機(jī)理

(1)

圖1 節(jié)理凸起接觸示意圖

(2)

(3)

(4)

1.2 三維節(jié)理的本構(gòu)方程

由接觸面的本構(gòu)方程，即可以得出節(jié)理的整體本構(gòu)方程，但是需知節(jié)理上接觸點(diǎn)在空間的分布情況。如圖2所示，設(shè)節(jié)理單位面積內(nèi)的接觸點(diǎn)的分布密度為ρc(與節(jié)理上下壁面形態(tài)、應(yīng)力水平有關(guān))，而接觸點(diǎn)在高程上的分布密度函數(shù)為H(z)，則ρcH(z)dz則是指單位面積節(jié)理上的高程z～z+dz的所有接觸點(diǎn)數(shù)目。由此可以得出，此時(shí)單位面積節(jié)理上接觸點(diǎn)的數(shù)目為

圖2 節(jié)理表面高程分布

(5)

而對(duì)于節(jié)理上接觸點(diǎn)的分布，應(yīng)根據(jù)具體的情況選用適當(dāng)?shù)姆植己瘮?shù)來描述，例如：卡方分布χ2、伽馬分布Γ(α,β),F分布等。

圖3 接觸面與作用力

(6)

(7)

式中：Δz=[zmin,zmax]。將節(jié)理凸起的接觸變形方程式(3)代入式(7)，得

(8)

(9)

這里剛度系數(shù)Kij為此時(shí)節(jié)理變形的剛度系數(shù)。由式(9)可以看出節(jié)理的變形剛度系數(shù)Kij與節(jié)理的形態(tài)特征、應(yīng)力狀態(tài)密切關(guān)。由于在巖石節(jié)理的變形過程中，節(jié)理上下兩壁之間的接觸面的分布、形態(tài)等是動(dòng)態(tài)變化的，因此變形剛度矩陣Kij是動(dòng)態(tài)變化的。三維節(jié)理的變形剛度矩陣為3×3的矩陣，非對(duì)角線上的剛度系數(shù)一般非零-由式(9)可知，在節(jié)理為平直節(jié)理時(shí)，非對(duì)角線上的剛度系數(shù)為零；或者對(duì)于傾角在傾向上為對(duì)稱分布，且當(dāng)節(jié)理當(dāng)前應(yīng)力狀態(tài)為0時(shí)，或僅法向應(yīng)力存在時(shí)節(jié)理變形矩陣中的非對(duì)角線上的剛度系數(shù)為零；對(duì)于其他情形，需要根據(jù)節(jié)理面的形態(tài)進(jìn)行具體的分析。然而，這里并不給直接給出節(jié)理本構(gòu)方程的具體形式，且對(duì)于動(dòng)力問題，應(yīng)力波在三維非線性本構(gòu)節(jié)理上折反射求解異常復(fù)雜，不利于理論解析。然而對(duì)于低幅值應(yīng)力波入射時(shí)，應(yīng)力波動(dòng)幅值相對(duì)節(jié)理的強(qiáng)度較小時(shí)，可以認(rèn)為此時(shí)節(jié)理的變形為線剛度模型，因此這里對(duì)巖石節(jié)理本構(gòu)方程為三維線剛度模型時(shí)的應(yīng)力波折反射規(guī)律進(jìn)行分析。

2 三維節(jié)理的應(yīng)力波折反射理論分析

令節(jié)理的平均面為XOZ坐標(biāo)平面，法向?yàn)閅軸；設(shè)應(yīng)力波入射線所在平面為XOY平面，則P波、SV波的偏振方向在XOY平面內(nèi)，而SH波的偏振方向?yàn)閆軸方向。節(jié)理上下兩壁巖石塊體A，B為線彈性介質(zhì)。由前面的討論可知，當(dāng)?shù)头祽?yīng)力波入射時(shí)，節(jié)理本構(gòu)方程可以采用的三維線剛度本構(gòu)方程dFi=Kijduj來表示。由線性問題的可疊加性可知，此時(shí)在應(yīng)力波入射與折、反射過程中，節(jié)理本構(gòu)方程可以改寫為

(10)

式中:[δσy,δτyx,δτyz]T為節(jié)理上的應(yīng)力變化分量，Pa；[δuy,δux,δuz]T為節(jié)理上下壁面的相對(duì)位移變化分量，m；kij，i,j=x,y,z為節(jié)理的變形剛度系數(shù)，Pa/m。且在本文中將處于矩陣K對(duì)角線上的kxx,kyy,kzz稱為主剛度系數(shù)，而位于對(duì)角線之外剛度系數(shù)稱為交叉耦合剛度系數(shù)。根據(jù)熱力學(xué)第二定律可以證明，節(jié)理剛度系數(shù)矩陣Kij必為實(shí)對(duì)稱陣，即kij=kji。由波動(dòng)理論可知，沿同一方向傳播的應(yīng)力波，P波、SV波與SH波之間的偏振方向是相互正交的，如圖4所示當(dāng)P波垂直入射時(shí)，節(jié)理面會(huì)產(chǎn)生擾動(dòng)δuy,δσy；SV波垂直入射時(shí)，會(huì)產(chǎn)生擾動(dòng)δux,δτxy；而P波、SV波傾斜入射時(shí)，都會(huì)在節(jié)理面產(chǎn)生擾動(dòng)δuy,δσy和δux,δτxy。而SH波無論是垂直入射還是傾斜入射僅會(huì)產(chǎn)生擾動(dòng)δuz,δτyz。垂直方向的擾動(dòng)δuy,δσy會(huì)因?yàn)楣?jié)理剛度系數(shù)kyy產(chǎn)生垂直方向的變形和應(yīng)力反饋δuy,δσy；再由波動(dòng)原理可知，一個(gè)平面上無數(shù)質(zhì)點(diǎn)作為振源輻射出垂直擾動(dòng)會(huì)形成平面波陣面，從而形成反射和透射的P波。類似的原理，垂直擾動(dòng)δuy,δσy會(huì)因節(jié)理剛度系數(shù)kxy反射和透射SV波，因剛度系數(shù)kyz反射和透射SH波。同理，剪切擾動(dòng)δux,δτxy和δuz,δτyz也能夠分別產(chǎn)生P波、SV波與SH波。

圖4 應(yīng)力波在三維節(jié)理上轉(zhuǎn)化機(jī)理示意圖

Fig.4 Mutual transformation mechanism of stress wave on 3D rock joint

因此對(duì)三維巖石節(jié)理，當(dāng)交叉耦合剛度系數(shù)不為零時(shí)，如圖5所示，P波、SV波或SH波會(huì)因剛度系數(shù)kxy,kyz,kxy而能夠相互轉(zhuǎn)化。

這里采用DDM(位移不連續(xù)方法)來研究應(yīng)力波折反射問題，且節(jié)理的變形本構(gòu)方程式(10)與應(yīng)力連續(xù)條件作為求解的邊界條件，入射波、反射波與透射波的波動(dòng)方程與符號(hào)約定如表1所示。

節(jié)理的變形[δuy,δux,δuz]T可以表示為

(11)

式中：[δuAy,δuAx,δuAz]T為巖石塊體A上節(jié)理的位移

(a) P波入射

(b) SV波入射

(c) SH波入射

變化量；[δuBy,δuBx,δuBz]T為巖石塊體B上節(jié)理的位移變化量。而由應(yīng)力波的傳播方向和質(zhì)點(diǎn)的偏振方向，可以得出巖石塊體A與塊體B的位移分量為

(12)

(13)

再由幾何變形方程，可以得出塊體A與塊體B中的應(yīng)變分量εAij與εBij。使用廣義胡克定律，即可以得出塊體A、塊體B中的應(yīng)力分量σAij與σBij。由DDM模型的假設(shè)：節(jié)理面上的應(yīng)力連續(xù)條件σAij|y=0=σBij|y=0和變形方程式(10)，并令應(yīng)力波的反射系數(shù)與透射系數(shù)為未知量，可以得出求解方程

X·E·C=ein·B

(14)

式中:X為求解方程的系數(shù)矩陣，與節(jié)理變形本構(gòu)方程和入射波的頻率、角度有關(guān)，具體形式為

(15)

式中：λA,μA為巖石塊體A的拉梅常數(shù);λB,μB為巖石塊體B的拉梅常數(shù)。E為指數(shù)矩陣

(16)

式中：C為反射應(yīng)力波與透射應(yīng)力波的幅值系數(shù)，具體形式與入射應(yīng)力波的類型有關(guān)，當(dāng)入射波為P波時(shí)，C為

(17)

式中：簡(jiǎn)寫為C=[RP→P,RP→SV,RP→SH,TP→P,TP→SV,TP→SH]T，類似的當(dāng)入射波為SV波時(shí)，折反射系數(shù)可以寫為C=[RSV→P,RSV→SV,RSV→SH,TSV→P,TSV→SV,TSV→SH]T；當(dāng)入射波SH波時(shí)，折反射系數(shù)可以寫為C=[RSH→P,RSH→SV,RSH→SH,TSH→P,TSH→SV,TSH→SH]T。

表1 波動(dòng)方程與符號(hào)約定

ein·B為式(14)的非齊次項(xiàng)，具體形式與入射波類型有關(guān)，當(dāng)入射波為P波時(shí)，ein·B可以寫為

einBP=eiκAd(xsin α1-CAdt)×

(18)

當(dāng)入射波為SV波時(shí)

(19)

當(dāng)入射波為SH波時(shí)

(20)

由穩(wěn)態(tài)應(yīng)力波場(chǎng)可知，應(yīng)力波在非連續(xù)面上傳播時(shí)應(yīng)該滿足Snell定律。對(duì)于三維節(jié)理面而言，Snell定律這里可以表示為

ω=κAdCAd=κAsCAs=κBdCBd=κBsCBs

(21)

因此，由以上分析和推導(dǎo)可知，對(duì)于三維粗糙節(jié)理面，由于剛度系數(shù)矩陣中交叉耦合項(xiàng)存在，簡(jiǎn)諧P波、SV波與SH波同時(shí)相互作用，之間能夠發(fā)生能量相互轉(zhuǎn)換，且SV波與SH波具有相同的反射角β2=γ2與透射角β3=γ3。將式(21)代入式(14)，即可以求出節(jié)理上的折反射系數(shù)C。

3 三維節(jié)理應(yīng)力波折反射規(guī)律分析

三維巖石節(jié)理對(duì)應(yīng)力波傳播的影響，與以往研究相比最大的不同點(diǎn)在于：三種類型的簡(jiǎn)諧波在節(jié)理上可以相互轉(zhuǎn)化，且本構(gòu)方程中剛度矩陣中交叉耦合項(xiàng)起到重要的作用?；谏瞎?jié)給出的應(yīng)力波在三維粗糙節(jié)理折反射的理論解，這里針對(duì)交叉耦合剛度系數(shù)對(duì)節(jié)理折反射應(yīng)力波的幅值系數(shù)影響進(jìn)行分析。

3.1 入射角對(duì)應(yīng)力波折反射系數(shù)的影響

圖6為P波入射時(shí)，折反射系數(shù)RP→P,TP→P,RP→SV,TP→SV,RP→SH和TP→SH隨入射角α1的變化規(guī)律。應(yīng)力波的幅值系數(shù)隨入射角的變化曲線與過去DDM(無交叉耦合剛度項(xiàng))所得的幅值系數(shù)類似，但是交叉耦合剛度系數(shù)的存在使得曲線有所變化，隨著節(jié)理中交叉耦合項(xiàng)的增大，與入射波同型的反射波的反射系數(shù)RP→P增大，且與入射波同型的透射波的幅值系數(shù)TP→P減小。若無交叉耦合剛度系數(shù)存在，則不會(huì)生成轉(zhuǎn)換波SH波；對(duì)于轉(zhuǎn)換波SV波，交叉耦合剛度系數(shù)的存在一定程度上改變了SV波折反射系數(shù)RP→SV,TP→SV的變化規(guī)律。但是總體上來說，隨著交叉耦合剛度系數(shù)的增加，轉(zhuǎn)換波的幅值系數(shù)是增大的。類似的規(guī)律在SV波入射時(shí)同樣成立，見圖7所示。SH波入射時(shí)，若無交叉耦合剛度系數(shù)存在，則轉(zhuǎn)換波P波與SV波則不會(huì)產(chǎn)生，如圖8所示；隨著交叉耦合剛度系數(shù)的增大，同型波反射系數(shù)RSH→SH增大，同型波透射系數(shù)TSH→SH減小，且轉(zhuǎn)換波折反射系數(shù)增大。

3.2 節(jié)理剛度系數(shù)對(duì)應(yīng)力折反射系數(shù)的影響

如圖9所示，P波垂直入射僅會(huì)導(dǎo)致節(jié)理面上的質(zhì)點(diǎn)在y軸方向運(yùn)動(dòng)，因此剛度系數(shù)kzx對(duì)應(yīng)力波傳播沒有影響，見圖9(c)。隨著剛度系數(shù)kxy的增大，與入射波同型波的反射系數(shù)RP→P和轉(zhuǎn)換SV波RP→SV，TP→SV的增大，與入射波同型波的透射系數(shù)TP→P減小，此時(shí)沒有SH波產(chǎn)生，見圖9(a)。類似的，當(dāng)變量為剛度系數(shù)kyz時(shí)，同型波的反射系數(shù)RP→P和轉(zhuǎn)換SH波RP→SH，TP→SH隨著kyz的增大而增加，同型入射波TP→P減小，且沒有SV波形成，見圖9(b)。相似的規(guī)律在SV波與SH波垂直入射時(shí)同樣存在，見圖10與圖11，這里不再贅述。

4 結(jié) 論

本文在討論三維巖石節(jié)理本構(gòu)方程的基礎(chǔ)上，通過合理的簡(jiǎn)化，分析了平面諧波在節(jié)理面上的作用機(jī)理，并給出了平面諧波在三維節(jié)理上的折反射理論解。并通過參數(shù)分析初步探討了P波、SV波和SH波入射時(shí)，反射波與透射波的幅值系數(shù)隨入射角以及交叉耦合剛度系數(shù)的變化規(guī)律。本研究得出以下有意義結(jié)論：

(1) 簡(jiǎn)諧波P波，SV波或SH波在三維節(jié)理上可以相互轉(zhuǎn)化，由于節(jié)理的交叉耦合變形剛度系數(shù)的存在，使得P波、SV波與SH波在節(jié)理面上可以發(fā)生波動(dòng)能量的相互傳遞和轉(zhuǎn)化。

(2) 三維粗糙節(jié)理本構(gòu)方程中的交叉耦合剛度系數(shù)越大，則越多的入射波能量轉(zhuǎn)化為轉(zhuǎn)換波能量；且同型反射波幅值隨交叉耦合剛度系數(shù)增大而增大，同型透射波幅值隨之減小。

(a) RP→P

(b) TP→P

(c) RP→SV

(d) TP→SV

(e) RP→SH,TP→SH

(a) RSV→SV

(b) TSV→SV

(c) RSV→P

(d) TSV→P

(e) RSV→SH,TSV→SH

(a) RSH→SH

(b) TSH→SH

(c) RSH→P

(d) TSH→P

(e) RSH→SV

(f) TSH→SV

(a) 變量為kxy,且此時(shí)kzy=kzx=0

(b) 變量為kzy,且此時(shí)kxy=kzx=0

(c) 變量為kzx,且此時(shí)kzy=kxy=0

(a) 變量為kxy,且此時(shí)kzy=kzx=0

(b) 變量為kzy,且此時(shí)kxy=kzx=0

(c) 變量為kzx,且此時(shí)kzy=kxy=0

(a) 變量為kxy,且此時(shí)kzy=kzx=0

(b) 變量為kzy,且此時(shí)kxy=kzx=0

(c) 變量為kzx,且此時(shí)kzy=kxy=0

(3) 相對(duì)于以往采用僅存在主剛度系數(shù)節(jié)理描述應(yīng)力波折反射規(guī)律的研究，節(jié)理中交叉耦合剛度系數(shù)的存在相當(dāng)于是對(duì)原模型的擴(kuò)展，交叉耦合剛度系數(shù)控制了應(yīng)力波能量流向。

本文對(duì)三維巖石節(jié)理應(yīng)力波傳播規(guī)律初步探討和分析，通過變形剛度系數(shù)考慮了由于復(fù)雜節(jié)理表面形態(tài)對(duì)應(yīng)力波折反射規(guī)律的影響。本研究有利于加深對(duì)巖體動(dòng)力學(xué)的理解，對(duì)于巖土地震工程和工程物探具有一定的積極意義。然而本文由于采用線剛度模型描述低幅值應(yīng)力波入射時(shí)的力-變形行為，忽略由于較大變形節(jié)理面嵌合度變化等現(xiàn)象帶來的非線性力學(xué)變化規(guī)律，因此不適用于高幅值應(yīng)力波入射的情形。因此下一步工作應(yīng)該結(jié)合數(shù)值計(jì)算與實(shí)驗(yàn)工作對(duì)高應(yīng)力幅值應(yīng)力波入射三維節(jié)理的傳播行為進(jìn)行研究。