王學(xué)濱，劉桐辛，田 鋒，錢(qián)帥帥

(1.遼寧工程技術(shù)大學(xué) 計(jì)算力學(xué)研究所，遼寧 阜新 123000； 2.遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院，遼寧 阜新 123000)

采礦、土木等工程中巖層周期破斷、斷層周期黏滑和反復(fù)爆破都能產(chǎn)生周期沖擊載荷。周期沖擊載荷是誘發(fā)巖爆、垮塌和片幫等災(zāi)害的重要因素。在周期沖擊載荷作用下，巷道頂板和兩幫會(huì)在極短時(shí)間內(nèi)產(chǎn)生大量相互聯(lián)通的裂縫，進(jìn)而將巷道圍巖切割成眾多塊體，它們會(huì)以猛烈方式脫離巷道圍巖，造成災(zāi)害[1-4]。

輕微巖爆發(fā)生時(shí)，可以聽(tīng)見(jiàn)清脆的噼啪、撕裂聲，偶有爆裂聲，巖塊彈射初速度小于1 m/s；而強(qiáng)烈?guī)r爆發(fā)生時(shí)，可以聽(tīng)見(jiàn)似炸藥爆破的爆裂聲，聲響強(qiáng)烈，巖塊彈射初速度處于5～10 m/s，極易造成重大危害。目前，關(guān)于巖塊彈射初速度已有不少來(lái)自實(shí)驗(yàn)室和現(xiàn)場(chǎng)的數(shù)據(jù)積累[5-6]。例如，王之東等[5]在單軸壓縮條件下對(duì)帶方孔的3種巖石試樣的能量釋放過(guò)程進(jìn)行了觀測(cè)，其中巖塊彈射初速度處于1.6～21.0 m/s；在加拿大和南非的巖爆現(xiàn)場(chǎng)[6]，巖塊彈射初速度處于7.65～12.60 m/s。對(duì)巖塊彈射初速度估算的方法主要包括理論方法和數(shù)值模擬研究方法。例如，在前者方面，秦劍鋒等[7]基于巖板的屈曲失穩(wěn)理論估算的巖塊彈射初速度在9 m/s；在后者方面，陳滔等[6]基于能量守恒原理和破壞前后模型的應(yīng)變能差估算的巖塊彈射初速度處于8.16～13.60 m/s?？陀^地講，巖爆過(guò)程較為復(fù)雜，受多種因素影響。在理論上想準(zhǔn)確估算巖塊彈射初速度極為困難。而且，在基于連續(xù)方法(例如，有限元法)的數(shù)值模擬中，僅能通過(guò)應(yīng)變能來(lái)估算巖塊彈射初速度，而不能模擬出巖塊從圍巖中脫離的過(guò)程，進(jìn)而難以準(zhǔn)確估算巖塊彈射初速度。

對(duì)于巷道的拉裂機(jī)理，目前已取得了一些重要進(jìn)展。例如，李夕兵等[8]采用PFC2D對(duì)動(dòng)載作用下深部巷道圍巖的動(dòng)力響應(yīng)進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)巷道頂板裂紋數(shù)量增加是由靜態(tài)的拉應(yīng)力與應(yīng)力波到達(dá)頂板后反射的拉應(yīng)力波疊加所引起；ZHU等[9-10]采用Autodyn2D對(duì)爆破誘發(fā)的巷道圍巖破壞過(guò)程進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)從自由邊界反射的拉應(yīng)力波引起了距離自由邊界一定距離的環(huán)向裂紋。

周期沖擊載荷3要素是幅值、頻率和持續(xù)時(shí)間。周期沖擊載荷3要素的影響研究一直受到重視[11]。前人已對(duì)周期沖擊載荷頻率的影響開(kāi)展了一定的研究。載荷頻率的取值范圍在不同文獻(xiàn)中不盡相同[12-13]。例如，閆長(zhǎng)斌等[12]采用FLAC在動(dòng)載作用下對(duì)地下巷道群的頻率影響進(jìn)行了數(shù)值模擬，頻率處于5～100 Hz；陳國(guó)祥等[13]采用FLAC對(duì)半正弦波作用下巷道圍巖的破壞過(guò)程進(jìn)行了數(shù)值模擬，頻率處于5～50 Hz?，F(xiàn)有的研究表明，載荷頻率對(duì)巖樣及巖石結(jié)構(gòu)的動(dòng)力響應(yīng)都有一定的影響[13-18]。例如，陳國(guó)祥等[13]發(fā)現(xiàn)隨著頻率的減小，巷道兩幫圍巖的最大垂直和水平應(yīng)力增加；宮鳳強(qiáng)等[16]在常規(guī)靜載和“預(yù)靜載+擾動(dòng)”條件下對(duì)巖石斷裂特性的頻率影響進(jìn)行了實(shí)驗(yàn)，發(fā)現(xiàn)隨著頻率的增加，斷裂韌度呈線性減小的趨勢(shì)；SU等[17-18]在真三軸加載條件下對(duì)含孔洞巖樣徑向應(yīng)力梯度的頻率影響進(jìn)行了實(shí)驗(yàn)，發(fā)現(xiàn)隨著頻率的增加，巖爆更易發(fā)生。目前，該方面研究主要集中在實(shí)驗(yàn)方面和基于連續(xù)方法的數(shù)值模擬方面，而基于連續(xù)-非連續(xù)方法的研究還十分少見(jiàn)。

鑒于單元畸變和局部自適應(yīng)阻尼可能導(dǎo)致單元彈射出模型時(shí)速度失真，為了準(zhǔn)確模擬彈射現(xiàn)象，在自主開(kāi)發(fā)的拉格朗日元與離散元耦合的連續(xù)-非連續(xù)方法[19]的基礎(chǔ)上，對(duì)彈射單元進(jìn)行了剛化處理，研究了周期沖擊載荷作用下巷道圍巖的變形-開(kāi)裂-運(yùn)動(dòng)過(guò)程，還初步分析了周期載荷頻率對(duì)單元彈射初速度和剪、拉裂縫區(qū)段數(shù)目的影響。

1 連續(xù)-非連續(xù)方法簡(jiǎn)介

1.1 原始方法簡(jiǎn)介

拉格朗日元和離散元耦合的連續(xù)-非連續(xù)方法主要包括4個(gè)計(jì)算模塊[19]：應(yīng)力應(yīng)變模塊、節(jié)點(diǎn)分離模塊、接觸力求解模塊和運(yùn)動(dòng)方程求解模塊。

在應(yīng)力應(yīng)變模塊中，采用了混合離散方法，通過(guò)節(jié)點(diǎn)的速度利用高斯定理求解單元的應(yīng)力和應(yīng)變，可以避免沙漏問(wèn)題。

在節(jié)點(diǎn)分離模塊中，通過(guò)引入虛擬裂縫模型，處理應(yīng)變軟化問(wèn)題。分別選取最大拉應(yīng)力準(zhǔn)則和莫爾-庫(kù)侖準(zhǔn)則作為拉裂和剪裂判據(jù)。當(dāng)節(jié)點(diǎn)分離后，虛擬裂縫產(chǎn)生，虛擬裂縫面之間存在法向及切向黏聚力。通過(guò)同時(shí)引入Ⅰ型和Ⅱ型斷裂能計(jì)算法向及切向黏聚力。當(dāng)法向或切向黏聚力降至0時(shí)，虛擬裂縫成為真實(shí)裂縫。

在接觸力求解模塊中，采用了基于空間劃分的接觸檢測(cè)方法和基于勢(shì)的接觸力求解方法，具有接觸檢測(cè)效率較高、無(wú)需對(duì)“角-角”接觸問(wèn)題進(jìn)行特殊處理的優(yōu)勢(shì)。

在運(yùn)動(dòng)方程求解模塊中，采用中心差分方法求解節(jié)點(diǎn)的速度和位移，具有計(jì)算效率高、計(jì)算精度較高的優(yōu)勢(shì)。

1.2 脫離模型的單元作剛體平動(dòng)的必要性及處理方法

網(wǎng)格法(有限元法、有限差分法等)容易出現(xiàn)網(wǎng)格或單元畸變問(wèn)題。本文的連續(xù)-非連續(xù)方法屬于網(wǎng)格法。當(dāng)某單元彈射出模型后，該單元的4個(gè)節(jié)點(diǎn)的速度一般并不相同，而且，各節(jié)點(diǎn)的運(yùn)動(dòng)是相互獨(dú)立的，所以，該單元的變形、運(yùn)動(dòng)規(guī)律將極為復(fù)雜，從而可能導(dǎo)致該單元發(fā)生畸變，這會(huì)使應(yīng)力、應(yīng)變和由應(yīng)力引起的彈性力等計(jì)算錯(cuò)誤，從而會(huì)進(jìn)一步加劇該單元畸變。當(dāng)某單元彈射出模型后，在與其他單元碰撞之前和之后，該單元應(yīng)僅受重力作用，而在水平方向不受力，這樣，該單元的各節(jié)點(diǎn)水平速度vx應(yīng)是常量。一旦該單元發(fā)生畸變，將不能保證這一點(diǎn)。另外，本文方法使用局部自適應(yīng)阻尼(與常見(jiàn)的黏性阻尼不同)，由于其方向取決于節(jié)點(diǎn)速度的方向，而大小取決于不平衡力(由應(yīng)力引起的彈性力和重力等構(gòu)成)的大小，這也可能導(dǎo)致vx不是常量。

為此，筆者提出一種處理方法以確保彈射單元各節(jié)點(diǎn)速度的計(jì)算結(jié)果正常，即在與其他單元碰撞之前和之后，保持恒定。

對(duì)于彈射單元，首先，需找到應(yīng)力狀態(tài)較為接近于自由狀態(tài)的臨界時(shí)步數(shù)目?？紤]到該單元?jiǎng)偙粡椛涑鰰r(shí)常處于壓縮狀態(tài)，因而其剛進(jìn)入完全拉伸狀態(tài)時(shí)的應(yīng)力狀態(tài)較為接近于自由狀態(tài)。為此，將該單元最小主應(yīng)力σ1>0(在本文中，σ3≥σ1，σ3為最大主應(yīng)力)對(duì)應(yīng)的時(shí)步數(shù)目作為臨界時(shí)步數(shù)目。然后，將該單元4個(gè)節(jié)點(diǎn)的速度矢量取平均，獲得平均速度矢量，并賦給這些節(jié)點(diǎn)，并對(duì)該單元的應(yīng)力清零，即清除彈性力。各節(jié)點(diǎn)的新速度v0取為

(1)

式中，vi(i=1～4)為該單元各節(jié)點(diǎn)的原速度。

上述處理迫使彈射單元在碰撞之前和之后作剛體平動(dòng)。由于彈射單元的應(yīng)力已被清除，單元畸變不會(huì)進(jìn)一步發(fā)展。而且，由于各節(jié)點(diǎn)只受重力(或?yàn)椴黄胶饬Φ奈ㄒ怀煞?和局部自適應(yīng)阻尼力(僅位于垂直方向上，其大小取決于重力的大小，其方向取決于節(jié)點(diǎn)垂直方向速度vy)，各節(jié)點(diǎn)在水平方向上將作勻速直線運(yùn)動(dòng)。這樣，即可避免節(jié)點(diǎn)速度的計(jì)算過(guò)于失真，從而確保巖塊彈射現(xiàn)象的模擬結(jié)果盡量真實(shí)。

2 計(jì)算模型及方案

本文的計(jì)算模型(以下簡(jiǎn)稱(chēng)模型)建立依據(jù)某大巷[20]所處地質(zhì)條件，巖性中硬。

模型尺寸為40 m×40 m，被剖分成160×160個(gè)正方形單元。壓應(yīng)力波施加后的力學(xué)模型如圖1(a)所示。應(yīng)當(dāng)指出，模型的左、右側(cè)面均為透射邊界，即應(yīng)力波經(jīng)過(guò)時(shí)不會(huì)發(fā)生反射。

在巷道頂板，布置了5個(gè)測(cè)點(diǎn)，在巷道左幫，布置了6個(gè)監(jiān)測(cè)單元，具體如圖1(b)所示。計(jì)算在平面應(yīng)變、大變形條件下進(jìn)行。

圖1 力學(xué)模型、測(cè)點(diǎn)和監(jiān)測(cè)單元Fig.1 Mechanical model,monitored points and monitored elements

陳建君等[21]將沖擊載荷簡(jiǎn)化為半正弦波，采用的半正弦壓應(yīng)力波P(N)表達(dá)式為

(2)

式中，Pmax為壓應(yīng)力波幅值，取7.3 MPa；頻率f=ω/π；ω為圓頻率；N為時(shí)步數(shù)目。

計(jì)算可分為3個(gè)過(guò)程：

(1)對(duì)開(kāi)挖前的模型進(jìn)行計(jì)算，直到模型處于靜力平衡狀態(tài)(所用N=12 000)。

(2)開(kāi)挖尺寸為6 m×6 m的巷道(所用N=4 000)；此后，繼續(xù)計(jì)算，直到模型處于靜力平衡狀態(tài)(所用N=4 000)。

(3)當(dāng)N=20 000時(shí)，在模型的上端面施加豎直向下的周期沖擊載荷(半正弦壓應(yīng)力波)，此后，繼續(xù)計(jì)算；當(dāng)N=64 000時(shí)最后1個(gè)壓應(yīng)力波波尾傳入模型，此后一段時(shí)間之內(nèi)計(jì)算仍在繼續(xù)。在最后1個(gè)壓應(yīng)力波在圍巖中逐漸消失的過(guò)程中，由于圍巖應(yīng)力場(chǎng)的略微調(diào)整，一些裂縫仍有可能進(jìn)一步發(fā)展，而且，脫離圍巖的單元仍在運(yùn)動(dòng)。因此，多計(jì)算一段時(shí)間是必要的。

共采用4個(gè)計(jì)算方案。方案1～4的f分別為15，25，35及45 Hz。文獻(xiàn)[22]通過(guò)現(xiàn)場(chǎng)觀測(cè)發(fā)現(xiàn)，震動(dòng)波的頻率主要集中于50 Hz以內(nèi)。本文選取的f涵蓋了上述范圍。

3 計(jì)算結(jié)果及分析

3.1 周期沖擊載荷作用下巷道圍巖變形-開(kāi)裂-運(yùn)動(dòng)過(guò)程

3.1.1多個(gè)壓應(yīng)力波沖擊下剪、拉裂縫的時(shí)空分布

以方案1為例簡(jiǎn)單分析多個(gè)壓應(yīng)力波沖擊下剪、拉裂縫的時(shí)空分布。

圖2為方案1的剪裂縫與最大主應(yīng)力σ3的時(shí)空分布規(guī)律，黑色線段代表剪裂縫區(qū)段。圖3為方案1的拉裂縫與最大主應(yīng)力σ3的時(shí)空分布規(guī)律，黑色線段代表拉裂縫區(qū)段。由圖2，3可以發(fā)現(xiàn)，首先，第1個(gè)壓應(yīng)力波傳至巷道兩幫后，巷道兩幫產(chǎn)生剪裂縫，并逐漸發(fā)展形成V形坑，其內(nèi)產(chǎn)生少量拉裂縫(圖2(a)，(b)和圖3(a)，(b))。然后，后續(xù)壓應(yīng)力波陸續(xù)傳入模型，巷道兩幫既有V形坑外若干新的V形坑形成，其內(nèi)仍產(chǎn)生少量拉裂縫，與此同時(shí)，巷道頂板產(chǎn)生拉裂縫，并發(fā)展形成層裂結(jié)構(gòu)，巷道兩幫脫離圍巖的一些單元彈入巷道(圖2(b)，(c)和圖3(b)，(c))。文獻(xiàn)[1]采用真三軸試驗(yàn)機(jī)對(duì)含預(yù)制矩形巷道的立方體巖樣進(jìn)行壓縮實(shí)驗(yàn)，發(fā)現(xiàn)巷道兩幫出現(xiàn)巖塊彈射現(xiàn)象，巖塊明顯堆積于巷道底板；文獻(xiàn)[23]利用水泥類(lèi)膨脹膠凝材料與水反應(yīng)體積驟增的特點(diǎn)對(duì)巷道圍巖進(jìn)行單次沖擊實(shí)驗(yàn)，發(fā)現(xiàn)巷道頂板和兩幫均發(fā)生破壞，且兩幫存在V形坑。上述現(xiàn)象的條件盡管與本文有所不同，但上述現(xiàn)象與本文的結(jié)果(圖2(b)，(c)和圖3(b)，(c))基本一致。最終，最后1個(gè)壓應(yīng)力波傳出模型后，即在圍巖中消失后，剪、拉裂縫停止發(fā)展。本文方法在處理開(kāi)裂、接觸和摩擦等非線性問(wèn)題時(shí)不可避免存在微小誤差，隨著計(jì)算的進(jìn)行，這種誤差可能被放大，這會(huì)導(dǎo)致剪、拉裂縫的分布不具有嚴(yán)格的對(duì)稱(chēng)性(例如，圖2(c)，(d)和圖3(c)，(d))。這種現(xiàn)象是非線性數(shù)值模擬方法的共性。應(yīng)當(dāng)指出，2個(gè)單元之間的裂縫稱(chēng)之為1個(gè)裂縫區(qū)段，裂縫區(qū)段的形狀為四邊形。若干裂縫區(qū)段連在一起構(gòu)成裂縫?？紤]到單元脫離圍巖后裂縫將變得很大，圖2～3僅顯示了各邊長(zhǎng)度均不大于1個(gè)單元邊長(zhǎng)的裂縫區(qū)段。

3.1.2巷道左幫監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx演化

圖4為方案1的1～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx-N曲線。由圖4可以發(fā)現(xiàn)，在壓應(yīng)力波傳入模型之后，隨著N的增加，2～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx(對(duì)于任一監(jiān)測(cè)單元，彈入巷道后其4個(gè)節(jié)點(diǎn)vx均相同)呈現(xiàn)上升—穩(wěn)定—衰減的變化過(guò)程，而1號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx呈現(xiàn)波動(dòng)—穩(wěn)定的變化過(guò)程。應(yīng)當(dāng)指出，當(dāng)N=20 000時(shí)第1個(gè)壓應(yīng)力波開(kāi)始傳入模型，當(dāng)N=64 000時(shí)最后1個(gè)壓應(yīng)力波波尾傳入模型。

圖4 方案1的1～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx-N曲線Fig.4 Evolution of horizontal velocities of lower right corner nodes of monitored elements 1-6 with time steps of Scheme 1

在vx有上升趨勢(shì)的階段，2～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx分別上升至最大值，此階段位于壓應(yīng)力波傳入模型之后，vx上升的過(guò)程是監(jiān)測(cè)單元逐漸脫離圍巖的過(guò)程。例如，當(dāng)N=28 120～31 990時(shí)，2號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx由0.396 m/s增至峰值6.8 m/s，此時(shí)，2號(hào)監(jiān)測(cè)單元脫離巷道圍巖(圖2(b))。期間，vx的最大值就是彈射初速度，分別為6.8，5.8，9.7，9.5和10.3 m/s，上述彈射初速度的平均值，即平均彈射初速度，為8.43 m/s。這位于強(qiáng)烈?guī)r爆[24]發(fā)生時(shí)平均巖塊彈射初速度范圍(5.0～10.0 m/s)之內(nèi)。

在vx穩(wěn)定的階段，2～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx保持不變。期間，這些監(jiān)測(cè)單元未與其他單元發(fā)生碰撞，即在水平方向上不受力，因而沒(méi)有產(chǎn)生能量損耗。由于巷道的空間有限，這些監(jiān)測(cè)單元一定會(huì)與其他單元發(fā)生碰撞，因此這些監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx無(wú)法總保持不變。這些監(jiān)測(cè)單元的位置和彈射初速度的不同導(dǎo)致vx發(fā)生變化的時(shí)刻不同，有的位于壓應(yīng)力波傳出模型之前，而有的位于壓應(yīng)力波傳出模型之后。

在vx有衰減趨勢(shì)的階段，2～6號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx呈現(xiàn)總體衰減的趨勢(shì)。此階段存在vx反向現(xiàn)象。首先，當(dāng)N=74 661～74 760時(shí)，2號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx由6.8 m/s降至-5.16 m/s。這是因?yàn)樵摫O(jiān)測(cè)單元與其他單元發(fā)生碰撞導(dǎo)致反向并產(chǎn)生能量損耗(圖2(c))；然后，該監(jiān)測(cè)單元不斷與其他單元發(fā)生碰撞，導(dǎo)致其右下角節(jié)點(diǎn)vx發(fā)生多次反向現(xiàn)象，并呈衰減趨勢(shì)(圖2(d))。

圖5為剛化處理前后2號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx-N曲線。為了呈現(xiàn)對(duì)彈射單元進(jìn)行剛化處理前后的差異，對(duì)方案1重新進(jìn)行了計(jì)算。由此可以發(fā)現(xiàn)，剛化處理前2號(hào)監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx雜亂無(wú)章，這并不是該單元與其他單元碰撞的結(jié)果，而剛化處理后2號(hào)監(jiān)測(cè)單元在未與其他單元碰撞時(shí)，可以保持恒定的vx，這較為符合該單元在水平方向上不受力的事實(shí)，這在一定程度上說(shuō)明了剛化處理的正確性。

圖5 剛化處理前后2號(hào)監(jiān)測(cè)單元的 右下角節(jié)點(diǎn)vx-N曲線Fig.5 Evolution of horizontal velocities of lower right corner nodes of monitored elements 2 with time steps before and after rigid treatments

3.1.3巷道頂板測(cè)點(diǎn)σ3的演化規(guī)律及開(kāi)裂機(jī)理

以方案1為例，簡(jiǎn)單分析巷道頂板各測(cè)點(diǎn)的σ3隨N的演化規(guī)律。圖6為方案1的1～5號(hào)測(cè)點(diǎn)的σ3-N曲線。由圖6可以發(fā)現(xiàn)，在壓應(yīng)力波傳入模型(N=20 000)之前，首先，各測(cè)點(diǎn)的σ3穩(wěn)定在-27 MPa，這對(duì)應(yīng)于巷道開(kāi)挖之前模型的靜力平衡狀態(tài)；隨后，受巷道開(kāi)挖的影響，各測(cè)點(diǎn)的σ3先經(jīng)歷震蕩上升、后逐漸衰減至穩(wěn)定的變化過(guò)程。在壓應(yīng)力波傳入模型之后，大部分測(cè)點(diǎn)的σ3呈現(xiàn)近似正弦波動(dòng)上升—衰減—穩(wěn)定的變化過(guò)程。

圖6 方案1的1～5號(hào)測(cè)點(diǎn)的σ3-N曲線Fig.6 Evolution of σ3 of monitored points 1-5 with time steps of Scheme 1

下面以4號(hào)測(cè)點(diǎn)為例，詳細(xì)分析近似正弦波動(dòng)上升階段σ3的演化規(guī)律。在此階段中，在σ3-N曲線上，可隱約觀察到13次波動(dòng)，這是因?yàn)榉桨?中施加了13個(gè)壓應(yīng)力波，各壓應(yīng)力波均會(huì)對(duì)該測(cè)點(diǎn)的σ3產(chǎn)生影響。圖7為方案1的4號(hào)測(cè)點(diǎn)的σ3-N曲線。由圖7可以發(fā)現(xiàn)：

圖7 方案1的4號(hào)測(cè)點(diǎn)的σ3-N曲線Fig.7 Evolution of σ3 of the monitored point 4 with time steps of Scheme 1

(1)σ3首先呈現(xiàn)有規(guī)律的波動(dòng)，然后波動(dòng)幅度突然增大，并伴隨著劇烈震蕩。σ3-N曲線的波峰和波谷均有隨著N的增加而增加的趨勢(shì)。

當(dāng)N=20 460～23 760時(shí)，σ3出現(xiàn)第1次波動(dòng)。這表明，第1個(gè)壓應(yīng)力波傳至并逐漸傳過(guò)4號(hào)測(cè)點(diǎn)。

當(dāng)N=20 460～22 130時(shí)，σ3由-15.8 MPa減小至-30.94 MPa。這表明第1個(gè)壓應(yīng)力波的峰前部分逐漸傳至該測(cè)點(diǎn)。應(yīng)當(dāng)指出，當(dāng)N=20 800～20 970時(shí)，σ3有小幅度震蕩。隨后，隨著N的增加，σ3繼續(xù)下降，σ3-N曲線斜率的絕對(duì)值明顯小于N=20 800之前的。也就是說(shuō)，σ3小幅度震蕩前后σ3-N曲線的斜率有所不同，前面曲線更陡，而后面更平緩。σ3小幅震蕩的原因是第1個(gè)壓應(yīng)力波的前緣經(jīng)由巷道頂板表面反射的拉應(yīng)力波傳至該測(cè)點(diǎn)。此后，第1個(gè)壓應(yīng)力波的后續(xù)部分繼續(xù)傳過(guò)該測(cè)點(diǎn)，并與反射的拉應(yīng)力波發(fā)生疊加，致使σ3-N曲線的斜率改變。

當(dāng)N=22 131～22 480時(shí)，σ3-N曲線出現(xiàn)第1個(gè)波谷，這是因?yàn)榈?個(gè)壓應(yīng)力波的峰值部分剛傳至該測(cè)點(diǎn)，與此同時(shí)，該測(cè)點(diǎn)仍位于第1個(gè)壓應(yīng)力波反射的拉應(yīng)力波中。

當(dāng)N=22 481～23 760時(shí)，第1個(gè)壓應(yīng)力波的峰后部分逐漸傳至該測(cè)點(diǎn)，σ3由-30.04 MPa增大至-11 MPa。當(dāng)N=23 760時(shí)，第1個(gè)壓應(yīng)力波已完全傳過(guò)該測(cè)點(diǎn)，與此同時(shí)，上述拉應(yīng)力波尚未完全傳出該測(cè)點(diǎn)，σ3-N曲線出現(xiàn)第1個(gè)波峰。第2個(gè)壓應(yīng)力波立即傳至該測(cè)點(diǎn)，導(dǎo)致σ3下降。應(yīng)當(dāng)指出，與4號(hào)測(cè)點(diǎn)剛受第1個(gè)壓應(yīng)力波影響時(shí)(N=20 460)的σ3相比，第1個(gè)壓應(yīng)力波剛傳出該測(cè)點(diǎn)時(shí)(N=23 760)的σ3更大，這是因?yàn)榈?個(gè)壓應(yīng)力波反射的拉應(yīng)力波的作用。

當(dāng)N=20 460～46 889時(shí)，σ3規(guī)律性波動(dòng)，期間，波峰和波谷隨著N的增加均有增加趨勢(shì)。波峰變化的原因同前所述。波谷變化的原因?yàn)椋悍瓷涞亩鄠€(gè)拉應(yīng)力波會(huì)對(duì)傳至該測(cè)點(diǎn)的壓應(yīng)力波有一定的衰減作用，致使壓應(yīng)力波的作用不再?gòu)?qiáng)烈。

當(dāng)N=46 890～47 040時(shí)，與此前幾次有規(guī)律的波動(dòng)相比，σ3的波動(dòng)幅度突然增大，震蕩加劇，這與該測(cè)點(diǎn)下方節(jié)點(diǎn)的分離(介質(zhì)的拉裂)有關(guān)。該測(cè)點(diǎn)下方節(jié)點(diǎn)的分離將使σ3向圍巖的深部轉(zhuǎn)移，從而將提升該節(jié)點(diǎn)的σ3。此外，該測(cè)點(diǎn)下方節(jié)點(diǎn)的分離將對(duì)應(yīng)力波的已有傳播產(chǎn)生影響，因而該測(cè)點(diǎn)σ3的震蕩將加劇。

(2)節(jié)點(diǎn)發(fā)生分離后，σ3的震蕩幅度有衰減趨勢(shì)。

當(dāng)N=53 190時(shí)，4號(hào)測(cè)點(diǎn)的σ3達(dá)到σt(5 MPa)，頂板在此處將發(fā)生拉裂。下面，闡明頂板的拉裂機(jī)理。若σ3的波峰不發(fā)生改變，則頂板不可能被拉裂。σ3的波峰隨著N的增加而增加使頂板拉裂成為可能。當(dāng)σ3極小時(shí)，波峰的σ3為負(fù)，代表壓縮；當(dāng)σ3極大時(shí)，波峰的σ3為正，代表拉伸。前一個(gè)壓應(yīng)力波剛完全傳過(guò)該測(cè)點(diǎn)之時(shí)，恰是下一個(gè)壓應(yīng)力波剛傳至該測(cè)點(diǎn)之時(shí)。下一個(gè)壓應(yīng)力波剛傳入該測(cè)點(diǎn)時(shí)的σ3總比前一個(gè)壓應(yīng)力波的大，這說(shuō)明該測(cè)點(diǎn)的最大σ3伴隨著應(yīng)力波的不斷傳入和傳出而越來(lái)越高，在不斷累積，這種累積只能源于反射的多個(gè)拉應(yīng)力波的作用，即反射的拉應(yīng)力波每通過(guò)一次該測(cè)點(diǎn)將其σ3提升一次，直至達(dá)到σt。

此后，σ3的震蕩規(guī)律變得復(fù)雜。與此同時(shí)，震蕩幅度有衰減的趨勢(shì)，直至穩(wěn)定在零值附近。水平展布的裂縫會(huì)阻斷應(yīng)力波的傳播，并引起應(yīng)力波的反射，從而引起該測(cè)點(diǎn)σ3的劇烈震蕩。經(jīng)由巷道頂板表面反射的拉應(yīng)力波產(chǎn)生的拉伸應(yīng)力應(yīng)在接近垂直方向上。所以，頂板將產(chǎn)生水平展布的拉伸裂縫。與此同時(shí)，若干傳入模型的壓應(yīng)力波相當(dāng)于增加了模型所受的垂直方向應(yīng)力。由此，頂板將產(chǎn)生垂直方向展布的拉伸裂縫。

此外，由圖6還可以發(fā)現(xiàn)，從整體上看，在近似正弦波動(dòng)上升階段，隨著N的增加，大部分測(cè)點(diǎn)的σ3波動(dòng)上升。越靠近巷道頂板表面，σ3的波動(dòng)幅度越小，其中，1號(hào)測(cè)點(diǎn)的最小，5號(hào)測(cè)點(diǎn)的最大。下面對(duì)其原因進(jìn)行解釋。當(dāng)?shù)?個(gè)壓應(yīng)力波的前緣傳至1號(hào)測(cè)點(diǎn)不久，在巷道頂板表面發(fā)生反射；隨后，由于1號(hào)測(cè)點(diǎn)離巷道頂板表面很近，反射的拉應(yīng)力波立即傳至1號(hào)測(cè)點(diǎn)，由于拉應(yīng)力波σ3的絕對(duì)值與此時(shí)傳至1號(hào)測(cè)點(diǎn)的壓應(yīng)力波的σ3相差應(yīng)不大，二者疊加將造成壓應(yīng)力波的能量極大衰減，即對(duì)1號(hào)測(cè)點(diǎn)的影響大幅降低。當(dāng)壓應(yīng)力波的前緣經(jīng)由巷道頂板表面反射的拉應(yīng)力波傳至5號(hào)監(jiān)測(cè)節(jié)點(diǎn)時(shí)，由于該拉應(yīng)力波在有阻尼的介質(zhì)中已傳播了一定距離，其能量將有所衰減，傳至該測(cè)點(diǎn)的壓應(yīng)力波的σ3的絕對(duì)值將大于該拉應(yīng)力波的，即該拉應(yīng)力波對(duì)5號(hào)測(cè)點(diǎn)的影響將不如對(duì)1號(hào)測(cè)點(diǎn)的大。

3.1.4頻率的影響

(1)對(duì)剪、拉裂縫區(qū)段數(shù)目的影響。

圖8，9分別為方案1～4的剪裂縫區(qū)段數(shù)目Ns和拉裂縫區(qū)段數(shù)目Nt隨N的演變規(guī)律，統(tǒng)計(jì)的Ns和Nt包括圖2～3中顯示與否的裂縫區(qū)段。

由圖8可以發(fā)現(xiàn)，各方案的Ns均呈現(xiàn)恒為0—近似階梯增長(zhǎng)—恒定的變化趨勢(shì)。只標(biāo)注了方案1的Ns-N曲線和Nt-N曲線。由圖9可見(jiàn)，各方案的Nt均呈現(xiàn)恒為0—近似階梯增長(zhǎng)—恒定的變化趨勢(shì)。

圖8 方案1～4的Ns-N曲線Fig.8 Evolution of the number of shear crack segments with time steps of Schemes 1-4

圖9 方案1～4的Nt-N曲線Fig.9 Evolution of the number of tensile crack segments with time steps of Schemes 1-4

在恒定為0階段，第1個(gè)壓應(yīng)力波的前緣尚未傳至巷道頂板表面。

在近似階梯增長(zhǎng)階段，Nt時(shí)增時(shí)穩(wěn)。以方案1(f=15 Hz)為例。當(dāng)N=23 200～27 520時(shí)，Nt-N曲線存在若干個(gè)增長(zhǎng)階梯，期間，巷道兩幫V形坑內(nèi)產(chǎn)生拉裂縫。當(dāng)N=27 521～46 300時(shí)，Nt穩(wěn)定在280附近，期間，V形坑內(nèi)拉裂縫幾乎停止發(fā)展，同時(shí)，巷道頂板尚未產(chǎn)生拉裂縫。在N=46 301之后，Nt快速增長(zhǎng)，期間，反射的多個(gè)拉應(yīng)力波的累積作用不斷提升頂板某一位置最大主應(yīng)力的波峰，巷道頂板產(chǎn)生拉裂縫。

在恒定階段，方案1～4的Nt分別最終穩(wěn)定在2 681，1 110，741和546。期間，模型中不再有壓應(yīng)力波傳入，拉裂縫停止發(fā)展。由此可見(jiàn)，隨著f的增加，Nt的最終穩(wěn)定值減小。

綜上所述，方案1～4的Ns和Nt的最終穩(wěn)定值均隨著f的增加而減小。已有研究表明，應(yīng)力波的頻率越高，在介質(zhì)中傳播衰減程度越大[25]。在本文模型中，采用局部自適應(yīng)阻尼，應(yīng)力波在其中傳播，能量自然也會(huì)衰減。上述Ns和Nt的最終穩(wěn)定值依賴于f的數(shù)值結(jié)果可以在一定程度上得到解釋。

(2)對(duì)平均彈射初速度的影響。

由圖4，10可以計(jì)算得到，方案1～4的平均彈射初速度分別為8.43，8.56，10.14和11.07 m/s。由此可見(jiàn)，隨著f的增加，平均彈射初速度增加。應(yīng)當(dāng)指出，方案1～4傳入模型的壓應(yīng)力波分別為13，22，30和33個(gè)。傳入的壓應(yīng)力波數(shù)量越多，則巷道圍巖對(duì)即將彈射單元的推動(dòng)作用越頻繁，這會(huì)導(dǎo)致平均彈射初速度越大。

圖10 方案2～4的監(jiān)測(cè)單元的右下角節(jié)點(diǎn)vx-N曲線Fig.10 Evolution of horizontal velocities of lower right corner nodes of monitored elements with time steps of Schemes 2-4

4 結(jié) 論

(1)鑒于單元畸變和局部自適應(yīng)阻尼可能導(dǎo)致單元彈射出模型時(shí)速度失真，為了準(zhǔn)確模擬彈射現(xiàn)象，在原始方法的基礎(chǔ)上對(duì)彈射單元進(jìn)行了剛化處理，確保了脫離模型的單元在碰撞前后作剛體運(yùn)動(dòng)。結(jié)果表明，彈射單元任一節(jié)點(diǎn)的水平速度呈現(xiàn)上升—穩(wěn)定—衰減的變化趨勢(shì)，這與不進(jìn)行剛化處理的結(jié)果相比更符合實(shí)際。

(2)在壓應(yīng)力波傳入模型之后，巷道頂板左、右對(duì)稱(chēng)線上大部分測(cè)點(diǎn)的最大主應(yīng)力呈現(xiàn)近似正弦波動(dòng)上升—衰減—穩(wěn)定的變化過(guò)程。在某一節(jié)點(diǎn)分離前，最大主應(yīng)力首先呈現(xiàn)有規(guī)律的波動(dòng)；然后，波動(dòng)幅度突然增大，并伴隨著劇烈震蕩。

(3)當(dāng)應(yīng)力波傳入模型之后，頂板某一位置的最大主應(yīng)力的波峰和波谷隨著時(shí)間的增加而有增加的趨勢(shì)，反射的多個(gè)拉應(yīng)力波的累積作用不斷提升該位置最大主應(yīng)力的波峰，致使拉裂，此即為巷道頂板拉裂機(jī)理。

(4)周期載荷頻率對(duì)單元彈射初速度和剪、拉裂縫區(qū)段數(shù)目有一定影響。隨著頻率的增加，脫離巷道圍巖的單元的平均彈射初速度增加，剪、拉裂縫區(qū)段數(shù)目減小。