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

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

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

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

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

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

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

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

1.1 原始方法簡介

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

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

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

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

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

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

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

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

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

(1)

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

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

2 計算模型及方案

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

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

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

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

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

(2)

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

計算可分為3個過程：

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

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

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

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

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

3.1 周期沖擊載荷作用下巷道圍巖變形-開裂-運動過程

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

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

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

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

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

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

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

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

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

圖5 剛化處理前后2號監(jiān)測單元的 右下角節(jié)點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巷道頂板測點σ3的演化規(guī)律及開裂機理

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

圖6 方案1的1～5號測點的σ3-N曲線Fig.6 Evolution of σ3 of monitored points 1-5 with time steps of Scheme 1

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

圖7 方案1的4號測點的σ3-N曲線Fig.7 Evolution of σ3 of the monitored point 4 with time steps of Scheme 1

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

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

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

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

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

當(dāng)N=20 460～46 889時，σ3規(guī)律性波動，期間，波峰和波谷隨著N的增加均有增加趨勢。波峰變化的原因同前所述。波谷變化的原因為：反射的多個拉應(yīng)力波會對傳至該測點的壓應(yīng)力波有一定的衰減作用，致使壓應(yīng)力波的作用不再強烈。

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

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

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

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

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

3.1.4頻率的影響

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

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

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

圖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個壓應(yīng)力波的前緣尚未傳至巷道頂板表面。

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

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

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

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

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

圖10 方案2～4的監(jiān)測單元的右下角節(jié)點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)致單元彈射出模型時速度失真，為了準確模擬彈射現(xiàn)象，在原始方法的基礎(chǔ)上對彈射單元進行了剛化處理，確保了脫離模型的單元在碰撞前后作剛體運動。結(jié)果表明，彈射單元任一節(jié)點的水平速度呈現(xiàn)上升—穩(wěn)定—衰減的變化趨勢，這與不進行剛化處理的結(jié)果相比更符合實際。

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

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

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