周期沖擊載荷下巷道頂板開裂機理數(shù)值模擬

王學濱，劉桐辛，田 鋒，錢帥帥

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

采礦、土木等工程中巖層周期破斷、斷層周期黏滑和反復爆破都能產(chǎn)生周期沖擊載荷。周期沖擊載荷是誘發(fā)巖爆、垮塌和片幫等災害的重要因素。在周期沖擊載荷作用下，巷道頂板和兩幫會在極短時間內(nèi)產(chǎn)生大量相互聯(lián)通的裂縫，進而將巷道圍巖切割成眾多塊體，它們會以猛烈方式脫離巷道圍巖，造成災害[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]基于能量守恒原理和破壞前后模型的應變能差估算的巖塊彈射初速度處于8.16～13.60 m/s?？陀^地講，巖爆過程較為復雜，受多種因素影響。在理論上想準確估算巖塊彈射初速度極為困難。而且，在基于連續(xù)方法(例如，有限元法)的數(shù)值模擬中，僅能通過應變能來估算巖塊彈射初速度，而不能模擬出巖塊從圍巖中脫離的過程，進而難以準確估算巖塊彈射初速度。

對于巷道的拉裂機理，目前已取得了一些重要進展。例如，李夕兵等[8]采用PFC2D對動載作用下深部巷道圍巖的動力響應進行了數(shù)值模擬，發(fā)現(xiàn)巷道頂板裂紋數(shù)量增加是由靜態(tài)的拉應力與應力波到達頂板后反射的拉應力波疊加所引起；ZHU等[9-10]采用Autodyn2D對爆破誘發(fā)的巷道圍巖破壞過程進行了數(shù)值模擬，發(fā)現(xiàn)從自由邊界反射的拉應力波引起了距離自由邊界一定距離的環(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)的動力響應都有一定的影響[13-18]。例如，陳國祥等[13]發(fā)現(xiàn)隨著頻率的減小，巷道兩幫圍巖的最大垂直和水平應力增加；宮鳳強等[16]在常規(guī)靜載和“預靜載+擾動”條件下對巖石斷裂特性的頻率影響進行了實驗，發(fā)現(xiàn)隨著頻率的增加，斷裂韌度呈線性減小的趨勢；SU等[17-18]在真三軸加載條件下對含孔洞巖樣徑向應力梯度的頻率影響進行了實驗，發(fā)現(xiàn)隨著頻率的增加，巖爆更易發(fā)生。目前，該方面研究主要集中在實驗方面和基于連續(xù)方法的數(shù)值模擬方面，而基于連續(xù)-非連續(xù)方法的研究還十分少見。

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

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

1.1 原始方法簡介

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

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

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

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

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

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

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

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

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

(1)

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

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

2 計算模型及方案

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

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

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

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

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

(2)

式中，Pmax為壓應力波幅值，取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)當N=20 000時，在模型的上端面施加豎直向下的周期沖擊載荷(半正弦壓應力波)，此后，繼續(xù)計算；當N=64 000時最后1個壓應力波波尾傳入模型，此后一段時間之內(nèi)計算仍在繼續(xù)。在最后1個壓應力波在圍巖中逐漸消失的過程中，由于圍巖應力場的略微調(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多個壓應力波沖擊下剪、拉裂縫的時空分布

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

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

圖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分別上升至最大值，此階段位于壓應力波傳入模型之后，vx上升的過程是監(jiān)測單元逐漸脫離圍巖的過程。例如，當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)測單元的位置和彈射初速度的不同導致vx發(fā)生變化的時刻不同，有的位于壓應力波傳出模型之前，而有的位于壓應力波傳出模型之后。

在vx有衰減趨勢的階段，2～6號監(jiān)測單元的右下角節(jié)點vx呈現(xiàn)總體衰減的趨勢。此階段存在vx反向現(xiàn)象。首先，當N=74 661～74 760時，2號監(jiān)測單元的右下角節(jié)點vx由6.8 m/s降至-5.16 m/s。這是因為該監(jiān)測單元與其他單元發(fā)生碰撞導致反向并產(chǎn)生能量損耗(圖2(c))；然后，該監(jiān)測單元不斷與其他單元發(fā)生碰撞，導致其右下角節(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)，在壓應力波傳入模型(N=20 000)之前，首先，各測點的σ3穩(wěn)定在-27 MPa，這對應于巷道開挖之前模型的靜力平衡狀態(tài)；隨后，受巷道開挖的影響，各測點的σ3先經(jīng)歷震蕩上升、后逐漸衰減至穩(wěn)定的變化過程。在壓應力波傳入模型之后，大部分測點的σ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個壓應力波，各壓應力波均會對該測點的σ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的增加而增加的趨勢。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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)鑒于單元畸變和局部自適應阻尼可能導致單元彈射出模型時速度失真，為了準確模擬彈射現(xiàn)象，在原始方法的基礎上對彈射單元進行了剛化處理，確保了脫離模型的單元在碰撞前后作剛體運動。結(jié)果表明，彈射單元任一節(jié)點的水平速度呈現(xiàn)上升—穩(wěn)定—衰減的變化趨勢，這與不進行剛化處理的結(jié)果相比更符合實際。

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

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

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