黃真萍，曾煥接，曹洋兵，陳俊熙

(1.福州大學(xué)環(huán)境與資源學(xué)院； 地質(zhì)工程福建省高校工程研究中心，福建 福州 350108； 2.國(guó)土資源部丘陵山地地質(zhì)災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，福建 福州 350108)

0 引言

巖體中存在大量性質(zhì)不同、規(guī)模各異的結(jié)構(gòu)面，彈性波傳到這些結(jié)構(gòu)面時(shí)，將發(fā)生反射、透射現(xiàn)象，并由此導(dǎo)致能量、振幅、波速等發(fā)生變化.鑒于彈性波在巖體中的傳播特性問(wèn)題涉及巖體工程爆破、聲波檢測(cè)及微震監(jiān)測(cè)等領(lǐng)域，眾多研究人員對(duì)此開(kāi)展大量卓有成效的研究，取得重要進(jìn)展.

起初，研究人員將巖體看作連續(xù)、均勻、各向同性的介質(zhì)，采用等效連續(xù)介質(zhì)力學(xué)方法研究彈性波傳播問(wèn)題.考慮到這種模型過(guò)于簡(jiǎn)化，與實(shí)際情況差距甚大，因此位移不連續(xù)模型[1-4]被提出，此模型認(rèn)為彈性應(yīng)力波穿過(guò)節(jié)理面時(shí)，應(yīng)力場(chǎng)是連續(xù)的，而位移是不連續(xù)的，從而較好地模擬應(yīng)力波在節(jié)理巖體中的傳播特性及過(guò)程.除此之外，結(jié)構(gòu)面對(duì)于巖體彈性波傳播的振幅衰減、高頻濾波及信號(hào)延遲等效應(yīng)[4-8]也被實(shí)驗(yàn)揭示，研究人員基于試驗(yàn)結(jié)果開(kāi)展相關(guān)的理論建模及力學(xué)描述研究.

由于室內(nèi)試驗(yàn)及解析計(jì)算等研究手段的局限性，研究人員嘗試?yán)秒x散元程序進(jìn)行結(jié)構(gòu)面對(duì)巖體彈性波傳播特性影響的研究[9-15]，離散元模擬方法逐漸成為一種重要的研究方法.茹忠亮等[9]得出節(jié)理面越粗糙，法向剛度越小，彈性波衰減程度越大，彈性波的波速降低越快.盧文波[10]指出節(jié)理剛度越大，其透射系數(shù)也越大； 石崇等[11]認(rèn)為節(jié)理面法向(或切向)剛度越大，各點(diǎn)透射率越高； Myer等[1]發(fā)現(xiàn)節(jié)理剛度越大，透射波形越接近于入射波，剛度越小，透射波波幅減小、高頻濾波.王衛(wèi)華等[12-13]發(fā)現(xiàn)節(jié)理面粗糙度相同時(shí)，透射系數(shù)隨法向剛度的增大而增大.陳勇等[14]認(rèn)為一維應(yīng)力波垂直節(jié)理傳播時(shí)，隨著節(jié)理法向剛度的增大，反射波波幅減小、頻率降低，透射波的波幅增大、頻率降低.

綜上可知，目前的研究集中在巖體結(jié)構(gòu)面對(duì)應(yīng)力波的透射系數(shù)、反射系數(shù)的影響及波的振幅衰減、高頻濾波、信號(hào)延遲等方面，而關(guān)于巖體結(jié)構(gòu)面法向剛度、切向剛度如何定量影響巖體彈性縱波波速及初至縱波振幅等問(wèn)題則鮮有研究報(bào)道.本文基于3DEC離散元分析平臺(tái)，通過(guò)連續(xù)介質(zhì)模型檢驗(yàn)軟件及自編程序的準(zhǔn)確性，基于單裂隙巖體模型，定量研究結(jié)構(gòu)面法向剛度、切向剛度對(duì)巖體彈性縱波傳播特性的影響規(guī)律，從而為相關(guān)巖體工程動(dòng)力響應(yīng)、安全性評(píng)價(jià)及數(shù)值模擬參數(shù)設(shè)置等提供理論支撐.

1 3DEC軟件及自編程序的可靠性

基于3DEC軟件和FISH語(yǔ)言，自編可精確監(jiān)測(cè)模型邊界動(dòng)力響應(yīng)的離散元計(jì)算程序，研究裂隙巖體彈性縱波傳播特性.以等效連續(xù)介質(zhì)模型的解析解作為標(biāo)準(zhǔn)參考值，檢驗(yàn)3DEC軟件的適用性及自編程序的準(zhǔn)確性.

圖1 巖塊數(shù)值模型Fig.1 Numerical model of rock block

對(duì)于連續(xù)、均勻、各向同性的彈性材料，其彈性縱波傳播速度Cp的理論解如下：

(1)

式中：ρ為材料的密度；E為材料的楊氏模量；μ為材料的泊松比.

1.1 數(shù)值模擬概況

構(gòu)建1 m×1 m×10 m的長(zhǎng)方體模型如圖1所示，將模型的頂面及四個(gè)立面設(shè)為自由邊界條件(由于底面施加縱波，側(cè)面的自由邊界條件不會(huì)有反射和折射，能達(dá)到粘性邊界條件的目的)，在模型底面施加速度型正弦縱波如下式 .

vz=Asin(ωt)

(2)

式中：vz為底面z方向的速度(m·s-1)；A為振幅，本次取0.1 m·s-1；ω為角頻率，本次取200 π，即頻率f為100 Hz.

表1 巖塊物理力學(xué)參數(shù)

根據(jù)巖體動(dòng)力響應(yīng)特征，選擇瑞利阻尼模型，設(shè)置最小臨界阻尼比為0.05、最小臨界中心頻率為100 Hz，設(shè)置的巖塊物理力學(xué)參數(shù)如表1所示.在模擬彈性波傳播時(shí)，網(wǎng)格單元的尺寸對(duì)計(jì)算結(jié)果的準(zhǔn)確性有很大影響，網(wǎng)格尺寸過(guò)大容易導(dǎo)致波形失真及波峰、波谷的位置偏離真實(shí)值，進(jìn)而導(dǎo)致波速計(jì)算值誤差過(guò)大.錢七虎院士[15]建議，網(wǎng)格尺寸與波長(zhǎng)之比不應(yīng)大于1/10～1/8.根據(jù)以上設(shè)置及式(1)，該巖塊的縱波速度理論解為3 651 m·s-1，故輸入波的波長(zhǎng)為36.51 m，由此設(shè)置網(wǎng)格尺寸為0.5 m，可以較精確地模擬波在模型中的傳播.

1.2 可靠性論證

經(jīng)計(jì)算發(fā)現(xiàn)，同樣的程序每次計(jì)算的結(jié)果都有輕微差別，這可能是數(shù)值截?cái)嗾`差等導(dǎo)致的離散性.經(jīng)過(guò)大量的計(jì)算嘗試，發(fā)現(xiàn)取6次數(shù)值計(jì)算值的平均值作為最終計(jì)算結(jié)果是穩(wěn)定、可靠的.因此，本文中所有列出的數(shù)值模擬結(jié)果數(shù)據(jù)都是取6次數(shù)值計(jì)算結(jié)果的平均值.

圖2 縱波波速理論值與模擬值隨楊氏模量E的變化 Fig.2 Theoretic calculations and numerical simulation results of P-wave velocity as a functionof the Young's modulus E

按照表1的參數(shù)設(shè)置，保持泊松比μ、密度ρ不變，依次改變楊氏模量E，得出的數(shù)值計(jì)算結(jié)果與理論解對(duì)比情況見(jiàn)圖2.由圖可知，數(shù)值模擬值與理論值相當(dāng)接近，兩者隨楊氏模量E變化的曲線也相當(dāng)吻合.同樣地，分別改變泊松比μ和密度ρ進(jìn)行數(shù)值模擬，結(jié)果也顯示彈性縱波波速理論值與模擬值較為接近，兩者的相對(duì)誤差皆小于2%.從而證實(shí)了程序的可行性和可靠性，并間接證明了上述數(shù)值模型尺寸、網(wǎng)格尺寸等程序設(shè)置的合理性，可應(yīng)用于后續(xù)結(jié)構(gòu)面巖體研究.

2 結(jié)構(gòu)面剛度影響的數(shù)值模擬

基于上述被論證為合理、可靠的3DEC數(shù)值模型及自編程序，為研究結(jié)構(gòu)面剛度對(duì)巖體彈性縱波傳播特性的影響，在連續(xù)介質(zhì)模型中部施加一條水平光滑的結(jié)構(gòu)面(為接觸面單元)，如圖3所示.巖塊的物理力學(xué)參數(shù)、邊界條件、阻尼模型及參數(shù)、網(wǎng)格尺寸等設(shè)置情況與上節(jié)相同.

圖3 單結(jié)構(gòu)面巖體數(shù)值模型 Fig.3 Numerical model of rock mass

2.1 研究方案

結(jié)構(gòu)面剛度是反映結(jié)構(gòu)面幾何構(gòu)成的函數(shù)，與結(jié)構(gòu)面的風(fēng)化蝕變特征、張開(kāi)度、粗糙度及其吻合情況等有關(guān).因此，平直結(jié)構(gòu)面和結(jié)構(gòu)面剛度可等效模擬復(fù)雜天然結(jié)構(gòu)面.

為研究巖體結(jié)構(gòu)面法向剛度、切向剛度對(duì)彈性波傳播特性的影響，需要首先確定法向剛度、切向剛度的合理取值范圍.需要說(shuō)明的是，3DEC軟件中結(jié)構(gòu)面剛度的基本單位為Pa·m-1，并非N·m-1.由結(jié)構(gòu)面的泊松效應(yīng)可知，結(jié)構(gòu)面法向剛度與切向剛度的比值應(yīng)為[1，10].一般來(lái)說(shuō)，對(duì)于離散元計(jì)算，法向剛度和切向剛度應(yīng)小于與該節(jié)理鄰接塊體的等效剛度的10倍，公式如下：

(3)

式中：K為塊體的體積模量；G為塊體的剪切模量； Δzmin為與該節(jié)理鄰接塊體的最小棱邊長(zhǎng)度.此要求一般僅為提高數(shù)值計(jì)算效率，對(duì)本文的簡(jiǎn)單模型而言，結(jié)構(gòu)面剛度超出此范圍不會(huì)顯著增加計(jì)算復(fù)雜度，反而能核實(shí)相關(guān)規(guī)定的正確性.

根據(jù)數(shù)值模擬基本情況，確定出的結(jié)構(gòu)面法向剛度總體取值范圍為7.5～7.5×106MPa·m-1，切向剛度總體取值范圍為3.0～3.0×106MPa·m-1，并通過(guò)以下基本思路研究結(jié)構(gòu)面剛度對(duì)巖體彈性縱波傳播特性的影響：

1) 保持結(jié)構(gòu)面切向剛度ks不變，將kn/ks由1逐漸增加至10.

2) 保持結(jié)構(gòu)面法向剛度kn不變，將kn/ks由1逐漸增加至10.

3) 保持結(jié)構(gòu)面法向剛度與切向剛度之比為2.5不變，同時(shí)改變法向剛度與切向剛度的大小.

2.2 數(shù)值模擬結(jié)果分析

2.2.1 結(jié)構(gòu)面法向剛度的影響分析

固定結(jié)構(gòu)面切向剛度為最低級(jí)(3 GPa·m-1)，將kn/ks由1逐漸增加至10，即逐漸增大結(jié)構(gòu)面法向剛度的大?。?再將結(jié)構(gòu)面切向剛度調(diào)至下一級(jí)，逐漸增加法向剛度； 最終將結(jié)構(gòu)面切向剛度調(diào)至最高級(jí)(30 GPa·m-1)，并同樣地逐漸增加法向剛度.從而在大跨度的切向剛度ks、法向剛度kn與切向剛度ks之比kn/ks范圍內(nèi)，研究法向剛度對(duì)彈性縱波傳播特性的影響.

不同級(jí)的切向剛度、不同kn/ks比值條件下，初至縱波振幅隨法向剛度kn的變化如圖4(a)所示，縱波波速隨法向剛度kn的變化如圖4(b)所示.給定激發(fā)縱波的振幅為100 mm·s-1，若無(wú)此結(jié)構(gòu)面，巖塊的彈性縱波波速理論值為3 651 m·s-1.

圖4 不同切向剛度下初至縱波振幅和縱波波速隨kn/ks的變化Fig.4 Amplitude of first break and P-wave velocity as a function of kn/ks at variousdiscontinuity tangential stiffness

由圖4(a)可知，當(dāng)切向剛度不變時(shí)，初至縱波振幅隨kn/ks比值的增加而增大，即初至縱波振幅隨法向剛度的增加而增大，但增長(zhǎng)速率逐漸減小，最終初至縱波振幅趨近于激發(fā)縱波振幅； 當(dāng)切向剛度固定為3～12 GPa·m-1時(shí)，初至縱波振幅在kn/ks比值為1～3.5之間為陡增段，且切向剛度越小，增長(zhǎng)速率越大，在kn/ks比值大于3.5后為緩增段，增長(zhǎng)速率逐漸減?。?當(dāng)切向剛度固定為12～30 GPa·m-1時(shí)，初至縱波振幅在kn/ks比值為1～10之間皆為緩增段，增長(zhǎng)速率逐漸減小.

由圖4(b)可知，當(dāng)切向剛度不變時(shí)，彈性縱波波速隨kn/ks比值的增加而增大，即彈性縱波波速隨結(jié)構(gòu)面法向剛度的增加而增大，但增長(zhǎng)速率逐漸減小，最后趨于穩(wěn)定，最終穩(wěn)定值介于3 566～3 670 m·s-1之間，與連續(xù)介質(zhì)彈性縱波波速理論值基本相當(dāng).

2.2.2 結(jié)構(gòu)面切向剛度的影響分析

由上節(jié)可知，初至縱波振幅和縱波波速與結(jié)構(gòu)面法向剛度密切相關(guān)，但無(wú)法判定其與結(jié)構(gòu)面切向剛度有關(guān).為此，在本節(jié)中，固定結(jié)構(gòu)面法向剛度為最低級(jí)(1.5 GPa·m-1)，將kn/ks比值由1逐漸增加至10，即逐漸減小結(jié)構(gòu)面切向剛度的大??； 再將結(jié)構(gòu)面法向剛度調(diào)至下一級(jí)，逐漸減小切向剛度； 最終將結(jié)構(gòu)面法向剛度調(diào)至最高級(jí)(30 GPa·m-1)，并同樣地逐漸減小切向剛度.從而在大跨度的法向剛度kn、法向剛度kn與切向剛度ks之比kn/ks范圍內(nèi)，研究切向剛度對(duì)巖體彈性縱波傳播特性的影響.

不同級(jí)的法向剛度條件下，初至縱波振幅隨kn/ks比值的變化如圖5(a)所示，縱波波速隨kn/ks比值的變化如圖5(b)所示.

由圖5可知，當(dāng)法向剛度不變時(shí)，初至縱波振幅和縱波波速隨kn/ks的增加都幾乎不變(曲線上局部微小的波動(dòng)是由于數(shù)值誤差引起).由此表明，結(jié)構(gòu)面切向剛度對(duì)巖體初至縱波振幅和縱波波速幾乎沒(méi)有影響.

圖5 不同法向剛度下初至縱波振幅和縱波波速隨kn/ks的變化Fig.5 Amplitude of first break and P-wave velocity as a function of kn/ks at variousdiscontinuity normal stiffness

2.2.3 考慮結(jié)構(gòu)面間距和巖塊楊氏模量的結(jié)構(gòu)面剛度的影響分析

由上兩節(jié)可知，結(jié)構(gòu)面法向剛度對(duì)巖體縱波傳播特性有重要作用，而切向剛度對(duì)縱波傳播幾乎沒(méi)有影響.但是，上述研究以巖塊縱波波速為基準(zhǔn)值，僅孤立地研究法向剛度絕對(duì)值的影響規(guī)律，忽視了結(jié)構(gòu)面法向剛度與巖塊楊氏模量之間的緊密關(guān)聯(lián)，因而難以將巖體彈性縱波傳播特性進(jìn)行整體考量.并且，巖體中的結(jié)構(gòu)面數(shù)量龐大，僅一條結(jié)構(gòu)面無(wú)法反映工程巖體實(shí)際情況.為此，引入比值E/(S·kn)，用以衡量巖塊楊氏模量E(本文設(shè)置為30 GPa)、結(jié)構(gòu)面間距S(對(duì)于本文的單結(jié)構(gòu)面巖體模型，S為巖體模型高度，即10 m)和結(jié)構(gòu)面法向剛度kn的組合關(guān)系對(duì)巖體彈性縱波傳播特性的綜合影響，從而使得研究結(jié)論具有更好的普適性.

結(jié)構(gòu)面法向剛度kn與切向剛度ks之比為2.5時(shí)，是硬巖工程的常見(jiàn)工況，故本次保持kn/ks比值為2.5不變，同時(shí)改變結(jié)構(gòu)面法向剛度kn與切向剛度ks的大小，研究E/(S·kn)對(duì)巖體彈性縱波傳播特性的影響.初至縱波振幅Ap/激發(fā)縱波振幅A0隨E/(S·kn)的變化如圖6(a)所示，縱波波速vp/連續(xù)介質(zhì)縱波波速理論值Cp隨E/(S·kn)的變化如圖6(b)所示.

圖6 Ap/A0和vp/Cp隨E/(S·kn)的變化Fig.6 Ap/A0 and vp/Cp as a function of E/(S·kn)

1) 由圖6(a)可知，初至縱波振幅Ap/激發(fā)縱波振幅A0隨E/(S·kn)的增大而減小，具體的減小規(guī)律為： ① 當(dāng)E/(S·kn)≤0.1時(shí)，初至縱波振幅衰減速率極小，衰減幅度不超過(guò)2%，Ap/A0接近于1； ② 當(dāng)0.110時(shí)，初至縱波振幅衰減速率越來(lái)越小，但總體衰減量越來(lái)越大，最終Ap/A0減小至0.01.

用三階多項(xiàng)式對(duì)圖6(a)的曲線進(jìn)行擬合，獲得以下擬合公式：

y1=-0.000 011 2x3+0.005x2-0.205x+0.994

(4)

式中：y1為初至縱波振幅Ap/激發(fā)縱波振幅A0比值；x為E/(S·kn).

2) 由圖6(b)可知，彈性縱波波速vp/連續(xù)介質(zhì)縱波波速理論值Cp隨E/(S·kn)的增大而減小，具體減小規(guī)律為： ① 當(dāng)E/(S·kn)≤0.1時(shí)，縱波波速衰減速率極小，衰減幅度不超過(guò)3%，vp/Cp接近于1； ② 當(dāng)0.110時(shí)，衰減速率越來(lái)越小，但總體衰減量越來(lái)越大，最終vp/Cp減小至0.67.

用四階多項(xiàng)式對(duì)圖6(b)的曲線進(jìn)行擬合，獲得以下擬合公式：

y2=0.000 003 16x4-0.001 4x3+0.058x2-0.28x+1.004 2

(5)

式中：y2為彈性縱波波速vp/連續(xù)介質(zhì)縱波波速理論值Cp比值；x為E/(S·kn).

3 結(jié)語(yǔ)

1) 結(jié)構(gòu)面法向剛度對(duì)巖體彈性縱波傳播特性有重要作用，而結(jié)構(gòu)面切向剛度對(duì)巖體彈性縱波傳播特性幾乎沒(méi)有影響.

2) 依次固定結(jié)構(gòu)面切向剛度為不同級(jí)大小(3 GPa·m-1到30 GPa·m-1)，在每級(jí)取值下，逐漸增大結(jié)構(gòu)面法向剛度，當(dāng)結(jié)構(gòu)面切向剛度不變時(shí)，初至縱波振幅和縱波波速皆隨法向剛度的增加而增大，但增長(zhǎng)速率逐漸減小，最終初至縱波振幅趨近于激發(fā)縱波振幅，縱波波速趨近于連續(xù)介質(zhì)縱波波速理論值.

3) 考慮巖塊楊氏模量E、結(jié)構(gòu)面間距S與結(jié)構(gòu)面法向剛度kn的綜合影響，引入比值E/(S·kn)，在保持kn/ks為2.5，同時(shí)改變結(jié)構(gòu)面法向剛度kn與切向剛度ks時(shí)，初至縱波振幅和縱波波速都隨E/(S·kn)的增加而減小，衰減速率經(jīng)歷極慢-變快-變慢-極慢四個(gè)階段，以E/(S·kn)=0.1、1、10為界，最終Ap/A0減小至0.01，vp/Cp減小至0.67.