李 春 鵬,李 遠 遠,仝 霄 金
(濟南市勘察測繪研究院,山東 濟南 250101)
天然巖體往往存在各種軟弱結構面,這些結構面內部被黏土、巖屑、巖塊等填充。由于充填介質與兩側巖體物理力學性質相差較大,導致結構面容易產生失穩破壞,尤其在動載作用下,當由地震或工程爆破產生的應力波到達軟弱結構面時會發生透反射,這一過程中應力波引起的動態響應可能會引發結構面滑移、錯動、傾覆等,對工程以及周邊場地造成不利的影響。
國內外學者對應力波入射無充填節理的傳播規律開展了大量研究工作。Nolte[1]、Myer[2]、Schoenberg[3]等基于位移不連續理論,研究了應力波入射結構面的透反射特性,分析了入射角、節理剛度、節理數量與厚度等對應力波傳播規律的影響;Zhao等[4]通過開展非線性變形不連續模型,研究了P波垂直入射單一節理的傳播特性;柴少波等[5]基于時域遞歸法,分析了P波入射交叉節理的傳播規律;盧文波[6]對應力波作用于滑移結構面的特征進行研究,揭示了滑移結構面的高頻濾波作用;劉嘯等[7]基于應力波入射結構面的透反射關系,建立了結構面剪切滑移失穩力學模型。
對于軟弱結構面,由于充填介質的影響使得應力波傳播更加復雜,部分學者開展了相關研究。李夕兵[8]考慮軟弱結構面摩擦滑移邊界條件,得到了軟弱結構面對應力波傳播的影響特征;Li[9]、劉婷婷[10]等基于三單元模型,采用特征線法得到了P波在軟弱結構面處的傳播規律,并討論了充填介質、接觸剛度等的影響特征;范留明等[11]采用射線理論分析了P波入射軟弱夾層的傳播特性,從能量角度探討了波阻比等因素的影響;賈帥龍等[12]基于時域遞歸法研究了充填節理對應力波傳播規律的影響,分析了充填節理法向與切向力學性質的變化規律。
此外,一些學者基于室內試驗開展了軟弱結構面相關研究。劉漢香等[13]通過開展振動臺試驗研究了含軟弱夾層邊坡的地震響應特征;周飛等[14]基于模型試驗分析了含水平軟弱夾層邊坡的動力響應特性;Wu[15]等采用SHPB試驗探討了含水率對P波在充填結構面巖體中衰減規律的影響。隨著計算機技術的發展,在數值模擬方面也有一些研究。廖少波等[16]采用3DEC研究了結構面產狀對巖質邊坡動力響應的影響;胡建華等[17]采用LS-DYNA分析了充填結構面對于爆炸應力波傳播的影響特征;Chen等[18]基于顆粒離散元法,研究了地震波在含軟弱夾層巖質邊坡中的傳播特征,進行了邊坡穩定性分析。
綜上所述,對于軟弱結構面應力波傳播特征問題,國內外學者大多采用理論分析方法開展相關研究,明確了應力波穿過軟弱結構面的透反射特性,但大部分研究并未考慮上下巖體與軟弱夾層的接觸關系,僅僅考慮了單個接觸面邊界條件的影響。因此,本文基于軟弱結構面理論模型,采用應力波理論對P波入射軟弱結構面傳播特征進行分析,并采用FEM數值模擬方法驗證理論結果的正確性,分別討論P波入射角、頻率、結構面剪切參數、充填材料力學參數以及充填厚度等因素對應力波傳播的影響特征,對比研究軟弱結構面滑移前后應力波透反射特性,為含軟弱結構面地質構造場地穩定性控制提供理論參考。
天然巖體中存在大量軟弱結構面,填充材料不同于巖體,往往是容易受到剪切破壞的砂礫、黏土等。由地震或者工程爆破在巖體中產生的應力波主要為P波和SV波,應力波從一側入射結構面時,會分別在結構面兩側發生透射與反射。由于P波易使巖體產生拉伸破壞,在研究軟弱結構面動力特性時優先考慮P波入射情況。基于此過程,建立由厚度為h的軟弱夾層以及無限巖體組成的軟弱結構面理論模型,如圖1所示,各組分均假定為線彈性、均勻、各項同性介質。

圖1 軟弱結構面理論模型Fig.1 Theoretical model of weak filled structural plane
設有一均勻平面簡諧P波從一側巖體(Ⅰ-巖體)入射軟弱結構面,入射角為α。在穿過結構面時會產生反射波和透射波,如圖1所示,其中在z<0區域內(I-巖體)存在入射P波、反射P波和反射SV波,振幅分別為Aip1、Arp1、Brs1;在0
在Ⅰ-巖體中應力波的位移勢函數可表示為
(1)
在Ⅱ-軟弱夾層中應力波的位移勢函數可表示為
(2)
在Ⅲ-巖體中應力波的位移勢函數可表示為
(3)
式中:φ(n)、ψ(n)分別為第n層介質P波、SV波位移勢函數(n=1,2,3);ω為應力波圓頻率;i為虛數單位;t為傳播時間;kwx、kwz分別為應力波在x和z方向的波矢量,w表示波型,沿傳播路徑的波矢量kw表示為
(4)
式中:cp、cs為介質P波波速與S波波速,對于各向同性彈性介質,可通過公式(5)計算得到
(5)
式中:λ為介質拉梅常數,μ為介質剪切模量,代入理論模型中的介質參數可得到相應波速。
根據Snell定理,各應力波在x方向的波矢量相等,即
kip1x=kip1sinα=krp1x=krs1x=ktp2x=kts2x
=krp2x=krs2x=ktp3x=kts3x=kx
(6)
對于同一介質中相同類型的應力波,波數是一致的,即
(7)
確定了位移勢函數形式后,需要結合邊界條件求解各應力波大小。由于軟弱結構面填充介質與兩側巖體性質相差較大,在接觸面上容易發生剪切滑移,這種剪切滑移行為一般采用摩爾-庫倫強度準則進行描述。基于此特點,可以將應力波入射軟弱結構面的動力作用過程分為結構面未滑移、結構面滑移兩個階段。對于結構面未滑移階段,軟弱夾層與巖體接觸面相互粘結,接觸面兩側介質位移與應力均滿足連續條件,即
(8)
(9)
對于結構面滑移階段,考慮兩側接觸面均發生滑移的情況,假定結構面未發生張裂,則接觸面兩側切向應力和切向位移不連續,且切向應力符合摩爾-庫倫準則,設結構面黏聚力為c、內摩擦角為θ,則有
(10)
(11)
上式中θe為等效內摩擦角,滿足以下關系[7]:
σzztanθe=σzztanθ+c
(12)
采用Helmholtz分解定理,可將基巖位移場u分解為矢量場φ與標量場ψ:
u=?φ+?×ψ
(13)
式(13)展開為x與z方向,得到
(14)
理論模型中各組分均假定為各向同性彈性介質,其應力應變關系應滿足線彈性胡克定律,應力由位移勢函數形式表示為
(15)
將式(1)~(3)代入式(13)、(14)中得到應力場、位移場,聯立式(8)~(11),結合式(4)~(7)條件得到振幅方程式。對于結構面未滑移的,有
(16)
對于結構面滑移的,有
(17)
矩陣中各分量見附錄。應力波的透反射關系一般采用透反射系數表示,即透反射波振幅與入射波振幅比值。定義該比值為Z,各應力波透反射系數可表示為
(18)
應力波理論解析的正確性可通過相似研究和數值方法進行驗證,在本文中采用有限元數值模擬進行對比驗證。ANSYS/LS-DYNA是一款顯式動力有限元軟件,適用于進行動力分析,因此用于本文模型驗證是可靠的。圖1中兩側巖體均為半無限巖體,結構面也具有延展性,這種無限邊界在數值模型中可通過設置透射邊界來實現,因此模型大小只要滿足計算要求即可。根據圖1理論模型概況,建立整體尺寸為4 m×4 m 的二維有限元模型,其中軟弱夾層厚度為5 cm,傾斜角度為10°,兩側巖體平均高度為197.5 cm。模型采用4節點Solid162實體單元,同時應用Lagrange網格進行劃分。為了更為準確地反映應力波的傳播,單元網格尺寸控制在應力波波長的1/6~1/12,本文選取網格尺寸為波長的1/10,考慮到軟弱夾層厚度較小,在劃分網格時適當減小尺寸,數值模型如圖2所示。

圖2 FEM數值模型(尺寸單位:m)Fig.2 FEM numerical model
由于理論模型中各組分均為線彈性,因此在數值模型中巖體與軟弱夾層均設置為線彈性材料(*Mat_Elastic),具體參數如表1所列。

表1 數值模型材料參數Tab.1 Material parameters of numerical model
在1.2節中計算了粘結邊界條件和滑移邊界條件兩種動力形式解,因此在數值模型中需要考慮兩種不同的邊界條件。如圖3所示,粘結邊界條件通過共節點方式實現,該方法要求接觸面兩側共用同一節點,滿足應力與位移連續,而滑移邊界條件通過添加接觸方式實現,在LS-DYNA中提供了摩爾-庫倫接觸方式。本文模型采用二維自動面-面接觸,設置等效摩擦角θe為10°。

圖3 粘結接觸(左)與滑移接觸(右)示意Fig.3 Schematic diagram of bonding contact (left) and sliding contact (right)
考慮簡單垂直入射情況,將輸入荷載加載于Ⅱ-軟弱結構面下邊界處。假定輸入位移是振幅為0.01 cm、頻率為200π的正弦波,只考慮一個周期,即0 ≤T≤0.01 s,表達式為
(19)
由數值模型可知軟弱結構面與初始波陣面的距離與角度相關,假定初始波陣面距離結構面L,應力波經過t0=L/crp入射結構面,此時應力波位移函數為
(20)
為了便于對比數值模擬與理論解析結果,定義同一水平位置下Ⅰ-巖體處接觸面位移峰值u0與Ⅱ-巖體處接觸面位移峰值ui之比為D,可表示為
(21)
將式(1)代入式(14)中得到Ⅰ-巖體處位移
(22)
將式(3)代入式(14)中得到Ⅲ-巖體處位移表達式為
(23)
聯立式(20)~(23),結合透反射系數可得到水平位移比Dx與垂直位移比Dz。選取數值模擬中相同水平位置下位移比D,得到粘結界面與滑移界面情況數值模擬與理論解析位移比如圖4所示。對于兩種接觸條件下,有限元數值模擬結果均與理論結果相近,誤差在允許范圍,表明應力波理論解析是準確可靠的,可用于后續分析。此外,對比滑移界面與粘結界面可知,結構面在粘結狀態下水平位移比約為垂直位移比的1.2倍,而結構面在滑移狀態下水平位移比約為垂直位移比的3.4倍,表明滑移失穩下上盤巖體水平位移過大,容易沿結構面向下運動造成滑坡等災害,因此需進一步明確P波入射軟弱結構面的傳播特征。

圖4 不同接觸條件下理論與數值模擬位移比DFig.4 Theoretical and numerical simulation D under different contact conditions
以2.1節中材料參數為基本物理參數,展開對平面P波入射軟弱結構面傳播特征研究。
在地震或工程爆破過程中應力波從震源發出,波陣面到達結構面時往往存在多個入射角。根據第1節應力波理論分析可知,不同入射角下P波透反射系數不同,因此在0°~90°范圍內均勻選取10個入射角,得到粘結邊界與滑移邊界條件透反射系數隨入射角變化關系如圖5所示。從圖中可以看到,當不同入射角P波入射結構面時,結構面在粘結狀態和滑移狀態下隨入射角呈現相似的變化趨勢,但從整體上看滑移狀態下透反射系數偏小,這是由于結構面滑動摩擦產生損耗導致的。當P波入射角為60°時,反射SV波達到峰值,但其余透反射波均為0,表明此時P波入射下只在Ⅰ-巖體中產生反射SV波,這種反射波應當為不均勻波。

圖5 不同入射角透反射系數Fig.5 Transmittance and reflection coefficient at different incident angles
對于處于粘結狀態的結構面,在入射角為0~60°之間,Ⅰ-巖體中反射P波一直處于較小值且呈減小趨勢,而反射SV波呈現快速增大趨勢,表明反射能量基本轉化為SV波,結構面處于剪切狀態;Ⅱ-軟弱夾層中,透射P波和反射P波呈現持續減小的趨勢,透射SV波和反射SV波呈現先增大后減小的趨勢,在50°~60°之間均表現為迅速衰減,表明在該角度范圍內軟弱夾層剪切作用逐漸顯著;Ⅲ-巖體中,透射SV波系數一直處于較小值,而透射P波在0°~50°之間一直保持為1,表明軟弱充填結構面并未對P波起到削弱作用。在入射角為60°~90°之間,Ⅰ-巖體中反射P波逐漸增大,在入射角為80°時開始迅速增長至1,反射SV波則逐漸減小至0,這與入射角為0~60°時完全相反;在Ⅱ-軟弱夾層與Ⅲ-巖體中,各透反射波均表現為先增大后減小的趨勢,并均在70°~80°之間達到峰值。因此,處于粘結狀態的軟弱結構面在P波入射角較大時受到擾動更小,尤其在入射角為50°~70°之間時。
對于處于滑移狀態下的結構面,相比于粘結狀態,Ⅰ-巖體中反射P波與Ⅲ-巖體中透射P波以及透射SV波表現較為不同。其中,Ⅰ-巖體中的反射P波隨入射角增大呈現波動趨勢,其反射系數遠大于粘結狀態,表明在滑移狀態下反射波能量顯著增長,而Ⅱ-軟弱夾層中透射SV波普遍較小,也證明結構面滑移會影響透反射能量分配。Ⅲ-巖體中透射P波在入射角為0°~50°之間并未保持恒定,而是呈現衰減趨勢,同樣在入射角為70°~80°之間透射系數仍未達到1,在這一傳播過程中P波能量部分提供給結構面滑動做功。此外,Ⅲ-巖體中透射SV波隨入射角變化呈現波動趨勢,相比粘結狀態具有較大的透射系數,表明滑移狀態下結構面剪切作用增強。綜合兩種狀態來看,當軟弱充填結構面從粘結狀態進入滑移狀態時,結構面整體響應會降低,但剪切作用增強。
當軟弱結構面離震源較近時,入射P波頻率較大,當軟弱結構面離震源較遠時,入射P波頻率較小,根據以往研究可知軟弱介質往往具有高頻濾波作用,因此需要討論軟弱結構面對頻率的敏感程度。選取入射波頻率(f)分別為10,100,1 000,10 000 Hz,其余參數保持不變。根據圖5結果,在Ⅲ-巖體中的透射P波在不同入射角下均保持較大值,因此著重分析應力波穿過結構面后透射P波的透射系數。不同入射波頻率下|Ztp3|隨入射角的變化趨勢如圖6所示。

圖6 不同頻率下|Ztp3|隨入射角變化趨勢Fig.6 Variation trend of |Ztp3| with incident angle at different frequencies
從圖6中可以看到,不同頻率入射P波作用下,粘結狀態和滑移狀態軟弱充填結構面表現為相同的變化趨勢。隨著頻率的增大,同一入射角下|Ztp3|呈現減小趨勢,其中當入射波頻率處于中低頻時,同一入射角下|Ztp3|未有太大衰減,當入射波頻率處于高頻時,同一入射角下|Ztp3|產生較大衰減,這充分體現了軟弱充填介質的高頻濾波特性。因此,當軟弱充填結構面離震源或爆源較遠時也可能發生失穩,當軟弱充填結構面離震源或爆源較近時將產生直接破壞。
在結構面的剪切滑移影響下,兩側巖體在水平向發生相對滑動,不同的剪切參數必然導致相對滑移的位移不同。從3.1節分析可知,結構面為滑移和粘結狀態下,透反射系數在入射角較大和較小時有差別,說明在不同接觸條件下應力波耗散程度不同,為了明確結構面剪切參數對應力波傳播的影響,計算等效摩擦角為0°,5°,10°,15°,20°時的透反射系數,如圖7所示。從圖7(a)中可知,同一P波入射角下,隨著等效摩擦角的增大,P波反射系數逐漸減小,而透射系數逐漸增大,但 Ⅱ-軟弱夾層中的P波透反射系數對等效摩擦角不敏感。從圖7(b)中可知,等效摩擦角為0°時S波透反射關系不同,與其他等效摩擦角最大的區別在于當P波入射角為0°時,SV波透反射系數均為0,從整體來看,隨著等效摩擦角增大,巖體中SV波系數逐漸減小,而軟弱夾層中SV波系數逐漸增大,這是軟弱夾層易發生剪切破壞的原因。綜合P波和S波透反射系數可知,隨著等效摩擦角增大,穿過軟弱結構面透射P波能量增大,透射S波能量減小,表明剪切變形越難發生。
充填材料的力學性質會影響應力波的透反射系數,特別是應力波穿過界面兩側介質的波阻抗,對應力波的傳播衰減影響較大。在結構面兩側巖體參數不變的情況下,考慮將充填材料波阻抗作為變量研究應力波的透反射關系。由于波阻抗存在密度和波速兩個變量,因此固定密度保持不變,只改變波速。為了更好描述材料力學參數的影響,定義充填材料波阻抗與巖體波阻抗之比為γ,可知在第3節材料參數下γ為0.291,在0.05~0.5之間均勻選取10個值,得到入射角為30°時不同γ影響下結構面的透反射系數的變化規律如圖8所示。

圖8 不同波阻抗比γ下透反射系數隨入射角變化趨勢Fig.8 Variation trend of transmittance and reflection coefficient with incident angle for different γ
從圖8(a)中可知,結構面透反射系數在粘結狀態下隨著波阻抗比γ逐漸增大,Ⅰ-巖體中反射波逐漸減小并趨于0;Ⅱ-軟弱夾層中透射P波逐漸增大,在γ>0.1后增長趨勢較緩且恒定,反射P波先增大后減小,在γ>0.1 后減小趨勢較緩且恒定,而反射SV波與透射SV波均隨γ增長并在γ>0.2后保持恒定;Ⅲ-巖體中透射P波在0<γ<0.25時迅速增大,而后維持在1.0保持不變,透射SV波表現為先增大后減小,在γ>0.1時逐漸減小并趨于0。由此可見,當P波穿過結構面引起的振動受γ控制,當γ>0.25后軟弱充填結構面無法起到削弱應力波的作用。
從圖8(b)中可知,結構面在滑移狀態下表現出和粘結狀態相似的變化趨勢,不同之處在于穿過結構面后的透射P波對γ更不敏感,即當γ>0.2后,透射系數保持在1.0附近,并且稍大于1.0,表明在滑移狀態下高波阻抗比將對入射P波起到放大作用。
應力波在含軟弱充填結構面處傳播時,充填介質厚度是不可忽視的。不同充填介質厚度使得應力波傳播路徑不同,由此導致的透反射特性也不同。為了充分考慮充填介質的厚度,在10~100 mm范圍內均勻選取10個充填厚度,分析粘結界面與滑移界面的透反射特性。入射波參數、介質參數與接觸參數見第3.0節,其中固定入射角為30°,得到不同充填厚度下透反射系數如圖9所示。

圖9 不同充填厚度下透反射系數隨入射角變化趨勢Fig.9 Variation trend of transmittance and reflection coefficient with incident angle for different filling thicknesses
從圖中可以看到粘結狀態和滑移狀態下結構面透反射系數隨著充填厚度的增加呈現波動變化,其中Ⅰ-巖體中反射P波和反射SV波波動趨勢與其余波正好相反。重點關注穿過結構面后的透射P波。從圖9中可以看到,當充填厚度在0.04,0.08 m時,透射系數都達最大值,對于粘結狀態下透射系數峰值相對一致,而滑移狀態下透射系數峰值相對降低。可見隨著充填厚度增加,透反射系數呈現周期變化,每隔0.04 m透射系數達到峰值,滑移狀態下透射系數呈減小趨勢。
本文通過應力波理論求解得到平面P波入射軟弱結構面時粘結狀態與滑移狀態下的透反射關系,并采用LS-DYNA模擬驗證了理論的正確性,最后分析了入射角、頻率、結構面剪切參數、充填材料性質與厚度等對P波傳播的影響特征,得到結論如下:
(1) 不同入射角下粘結狀態與滑移狀態下應力波傳播特征相近,但滑移狀態下應力波透反射系數稍低,這是由于滑移過程中存在能量損耗,此外入射角為60°時只存在反射S波,在該角度下只存在剪切破壞。
(2) 軟弱充填結構面具有高頻濾波特性,表明即使在低頻下,軟弱充填結構面也容易發生破壞,而高頻下產生的破壞更多是在入射波處巖體,因此對于工程爆破軟弱結構面起到削弱作用,導致結構面上盤容易產生大塊巖體,對于地震作用下需要時刻關注結構面兩側的穩定性。
(3) 等效摩擦角可控制軟弱結構面的剪切滑移,當等效摩擦角增大時,穿過軟弱結構面透射P波能量會隨之增大,而透射SV波能量減小,表明剪切變形越難發生。
(4) P波穿過軟弱結構面振幅受充填介質與巖體波阻抗比控制,粘結狀態下波阻抗比大于0.25時,透射系數基本保持恒定,滑移狀態下波阻抗比大于0.2時基本保持恒定。
(5) 隨著軟弱結構面厚度增大,應力波幅值呈現周期波動式變化,周期為0.04 m,該現象可能與入射P波波長相關,需要進一步研究。
(6) 巖土工程中應根據震源與軟弱結構面分布確定潛在滑移位置,在抗震方面需注重軟弱結構面的加固,在工程爆破方面需充分利用軟弱結構面選擇合理的爆破位置達到破巖、減振等效果。
附錄:
矩陣[a]中分量aij表達式為
a17=a18=0,a21=-2krp1zkxμr,
a31=-krp1z,a32=-a34=-a36=-kx,
a33=-ktp2z,a35=krp2z,a37=a38=0,
a41=-a43=-a45=-kx,a42=krs1z,a44=kts2z,
a46=-krs2z,a47=a48=0,a51=a52=0,
a54=-2eih(ktp2z+ktp3z+kts3z)kts2zkxμs,

a56=2eih(krs2z+ktp2z+ktp3z+kts2z+kts3z)krs2zkxμs,
a58=2eih(ktp2z+ktp3z+kts2z)kts3zkxμr,a61=a62=0,
a63=-2eih(ktp3z+kts2z+kts3z)ktp2zkxμs,
a65=2eih(krp2z+ktp2z+ktp3z+kts2z+kts3z)krp2zkxμs,
a67=2eih(ktp2z+kts2z+kts3z)ktp3zkxμr,
a71=a72=0,a73=-eih(ktp3z+kts2z+kts3z)ktp2z,
a74=eih(ktp2z+ktp3z+kts3z)kx,
a75=eih(krp2z+ktp2z+ktp3z+kts2z+kts3z)krp2z,
a76=eih(krs2z+ktp2z+ktp3z+kts2z+kts3z)kx,
a77=eih(ktp2z+kts2z+kts3z)ktp3z,
a78=-eih(ktp2z+ktp3z+kts2z)kx,a81=a82=0,
a83=-eih(ktp3z+kts2z+kts3z)kx,
a84=-eih(ktp2z+ktp3z+kts3z)kts2z,
a85=-eih(krp2z+ktp2z+ktp3z+kts2z+kts3z)kx,
a86=eih(krs2z+ktp2z+ktp3z+kts2z+kts3z)krs2z,
a87=eih(ktp2z+kts2z+kts3z)kx,a88=eih(ktp2z+ktp3z+kts2z)kts3z
矩陣[b]中分量bij表達式為
b31=-kip1z,b41=kx,b51=b61=b71=b81=0
矩陣[c]中分量cij表達式為
c87=eihkts3z{-2ktp3zkxμr
c43=c44=c45=c46=c47=c48=c81
=c82=c83=c84=c85=c86=0,
其余項與矩陣[a]對應項一致。
矩陣[d]中分量dij表達式為
其余項與矩陣[b]對應項一致。