結(jié)構(gòu)面剛度對巖體彈性縱波傳播特性影響的數(shù)值模擬

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

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

0 引言

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

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

由于室內(nèi)試驗及解析計算等研究手段的局限性，研究人員嘗試?yán)秒x散元程序進(jìn)行結(jié)構(gòu)面對巖體彈性波傳播特性影響的研究[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ù)隨法向剛度的增大而增大.陳勇等[14]認(rèn)為一維應(yīng)力波垂直節(jié)理傳播時，隨著節(jié)理法向剛度的增大，反射波波幅減小、頻率降低，透射波的波幅增大、頻率降低.

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

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

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

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

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

(1)

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

1.1 數(shù)值模擬概況

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

vz=Asin(ωt)

(2)

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

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

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

1.2 可靠性論證

經(jīng)計算發(fā)現(xiàn)，同樣的程序每次計算的結(jié)果都有輕微差別，這可能是數(shù)值截斷誤差等導(dǎo)致的離散性.經(jīng)過大量的計算嘗試，發(fā)現(xiàn)取6次數(shù)值計算值的平均值作為最終計算結(jié)果是穩(wěn)定、可靠的.因此，本文中所有列出的數(shù)值模擬結(jié)果數(shù)據(jù)都是取6次數(shù)值計算結(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ù)值計算結(jié)果與理論解對比情況見圖2.由圖可知，數(shù)值模擬值與理論值相當(dāng)接近，兩者隨楊氏模量E變化的曲線也相當(dāng)吻合.同樣地，分別改變泊松比μ和密度ρ進(jìn)行數(shù)值模擬，結(jié)果也顯示彈性縱波波速理論值與模擬值較為接近，兩者的相對誤差皆小于2%.從而證實了程序的可行性和可靠性，并間接證明了上述數(shù)值模型尺寸、網(wǎng)格尺寸等程序設(shè)置的合理性，可應(yīng)用于后續(xù)結(jié)構(gòu)面巖體研究.

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

基于上述被論證為合理、可靠的3DEC數(shù)值模型及自編程序，為研究結(jié)構(gòu)面剛度對巖體彈性縱波傳播特性的影響，在連續(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)化蝕變特征、張開度、粗糙度及其吻合情況等有關(guān).因此，平直結(jié)構(gòu)面和結(jié)構(gòu)面剛度可等效模擬復(fù)雜天然結(jié)構(gòu)面.

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

(3)

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

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

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

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

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

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

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

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

不同級的切向剛度、不同kn/ks比值條件下，初至縱波振幅隨法向剛度kn的變化如圖4(a)所示，縱波波速隨法向剛度kn的變化如圖4(b)所示.給定激發(fā)縱波的振幅為100 mm·s-1，若無此結(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)切向剛度不變時，初至縱波振幅隨kn/ks比值的增加而增大，即初至縱波振幅隨法向剛度的增加而增大，但增長速率逐漸減小，最終初至縱波振幅趨近于激發(fā)縱波振幅； 當(dāng)切向剛度固定為3～12 GPa·m-1時，初至縱波振幅在kn/ks比值為1～3.5之間為陡增段，且切向剛度越小，增長速率越大，在kn/ks比值大于3.5后為緩增段，增長速率逐漸減?。?當(dāng)切向剛度固定為12～30 GPa·m-1時，初至縱波振幅在kn/ks比值為1～10之間皆為緩增段，增長速率逐漸減小.

由圖4(b)可知，當(dāng)切向剛度不變時，彈性縱波波速隨kn/ks比值的增加而增大，即彈性縱波波速隨結(jié)構(gòu)面法向剛度的增加而增大，但增長速率逐漸減小，最后趨于穩(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)，但無法判定其與結(jié)構(gòu)面切向剛度有關(guān).為此，在本節(jié)中，固定結(jié)構(gòu)面法向剛度為最低級(1.5 GPa·m-1)，將kn/ks比值由1逐漸增加至10，即逐漸減小結(jié)構(gòu)面切向剛度的大小； 再將結(jié)構(gòu)面法向剛度調(diào)至下一級，逐漸減小切向剛度； 最終將結(jié)構(gòu)面法向剛度調(diào)至最高級(30 GPa·m-1)，并同樣地逐漸減小切向剛度.從而在大跨度的法向剛度kn、法向剛度kn與切向剛度ks之比kn/ks范圍內(nèi)，研究切向剛度對巖體彈性縱波傳播特性的影響.

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

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

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

結(jié)構(gòu)面法向剛度kn與切向剛度ks之比為2.5時，是硬巖工程的常見工況，故本次保持kn/ks比值為2.5不變，同時改變結(jié)構(gòu)面法向剛度kn與切向剛度ks的大小，研究E/(S·kn)對巖體彈性縱波傳播特性的影響.初至縱波振幅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時，初至縱波振幅衰減速率極小，衰減幅度不超過2%，Ap/A0接近于1； ② 當(dāng)0.110時，初至縱波振幅衰減速率越來越小，但總體衰減量越來越大，最終Ap/A0減小至0.01.

用三階多項式對圖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時，縱波波速衰減速率極小，衰減幅度不超過3%，vp/Cp接近于1； ② 當(dāng)0.110時，衰減速率越來越小，但總體衰減量越來越大，最終vp/Cp減小至0.67.

用四階多項式對圖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é)語

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

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

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