含不同半徑孔洞的顆粒體模型的力學行為數(shù)值模擬

伍小林，王學濱，潘一山

（遼寧工程技術大學 力學與工程學院，遼寧 阜新 123000）

0 引言

眾所周知，試樣的尺寸會影響其強度，甚至是峰后的應力-應變曲線。這種問題在固體力學和巖石力學中常稱之為尺寸效應或尺度律問題。這些問題的重要性怎么強調(diào)也不過分，這是因為在很多的物理和工程問題中，尺度律問題都占據(jù)其中心位置［1］。目前，關于不含孔洞的矩形試樣的尺寸效應問題研究得最為充分，尤其是單向壓縮條件下巖石、混凝土等準脆性材料的尺寸效應問題［1-6］。

尺寸效應也存在于巖石鉆孔的破裂之中［1］，這些現(xiàn)象的確認是通過室內(nèi)試驗得到的［7-8］。最近的現(xiàn)場觀測表明，巷道斷面的尺寸對巖爆的量級有重要的影響［9-10］。他們對某水電站探洞和交通洞開挖洞段的現(xiàn)場巖爆特征的調(diào)查統(tǒng)計表明，巖爆的基本地質條件基本相同，但是，在探洞（斷面為3m×3m）發(fā)生了II級巖爆（比I級巖爆弱），而交通洞（斷面為6m×6m）發(fā)生了I級巖爆，并由此推測，巖爆量級隨斷面開挖尺寸的增大而增強，巖爆具有尺寸效應。文獻［11］開展了不同尺寸的隧洞開挖之后，圍巖的剪切局部化過程的數(shù)值模擬研究。在計算中，采用了FLAC軟件的應變軟化本構模型，以“先加載，后挖洞”方式進行模擬。研究結果表明，對于小孔隧洞，圍巖中出現(xiàn)了4個對稱的V形坑，這與靜水壓力條件下有關的實驗結果和現(xiàn)場觀測結果較為吻合；對于大孔隧洞，圍巖中出現(xiàn)了多條狹長的剪切帶，圍巖發(fā)生了大面積的破壞。由此發(fā)現(xiàn)了隧洞的尺寸對圍巖的局部化有重要的影響。

本文將巖石視為非連續(xù)介質，即由顆粒和界面構成的顆粒體材料，采用FLAC對含不同半徑孔洞的顆粒體模型進行了數(shù)值模擬。在靜水壓力條件下，先將模型加載至平衡，然后再通過 FLAC的內(nèi)嵌編程FISH語言模擬開挖。研究了模型再次平衡后，巷道半徑不同時，圍巖中的剪切應變增量、最小主應力及最大主應力的分布規(guī)律，此外，還監(jiān)測了一些單元的一些力學量的分布規(guī)律，進一步確認了巷道的尺寸對圍巖變形破壞的重大影響。

1 顆粒的堆積算法

顆粒的堆積過程是在某一矩形區(qū)域內(nèi)由下至上進行的。為了易于實現(xiàn)顆粒堆積算法，需要事先在矩形區(qū)域下邊界上人為地放置一排顆粒。矩形區(qū)域的左邊界和右邊界均可視為半徑為無限大的顆粒。在本文中，表層顆粒是指眾多顆粒堆積成的堆積體的表面一層的顆粒。表層顆粒中的任一顆粒均未被其周圍的顆粒完全包圍，即顆粒的上方無任何顆粒而其下方至少有兩個顆粒與之相切。顆粒的堆積算法如下：

（1）記錄當前表層顆粒，即矩形區(qū)域下邊界上的一排顆粒。

（2）假設顆粒半徑服從均勻分布，隨機產(chǎn)生一個半徑為rn的新顆粒n。

（3）在表層顆粒中隨機抽取一個顆粒i。

（4）在顆粒i附近找到一個顆粒j（顆粒j可為左邊界、右邊界或表層顆粒），使得新顆粒n能放在顆粒i和j之上，并且與顆粒i和j都相切。

（5）根據(jù)幾何關系求出新顆粒n的中心坐標，設其中心坐標為（xn，yn）。由顆粒i和顆粒n相切的關系，可得：

式中：ri——顆粒i的半徑；

xi，yi——顆粒 i的中心坐標。

如果顆粒j與新顆粒n相切，那么將可能出現(xiàn)三種情況：

a）若顆粒j為左邊界，則有

b）若顆粒j為右邊界，則有

式中a為矩形區(qū)域水平方向的長度。

c）若顆粒j為表層顆粒，則有

式中：rj——顆粒j的半徑；

xj，yj——為顆粒 j的中心坐標。

由式（1）與式（2）、（3）及（4）中的任一個可求出新顆粒n的中心坐標（xn，yn）的兩組解：其中一組 yn值較小的解對應于新顆粒n在顆粒i和j下方的情況，因而應當舍去；另一組解即為所求。

（6）判斷新顆粒n的中心橫坐標xn是否介于顆粒i的中心橫坐標xi與顆粒j的中心橫坐標xj（或左邊界、右邊界）之間，以確保新顆粒n在重力作用下能保持靜止。若xn不介于它們之間，則再回到第（3）步；否則，進行下一步。

（7）本文不允許兩個顆粒相互侵入，因此需要判斷任意兩個顆粒是否相互侵入。只要有一個顆粒與新顆粒相互侵入，則再回到第（3）步；否則，進行下一步。

（8）修改當前表層顆粒的記錄。

（9）重復（2）到（8）步，直至結束。

圖1給出了利用上述算法獲得的顆粒的堆積過程。

圖1 顆粒堆積過程Fig.1 The accumulating process of granulaes

2 計算模型、計算方案及本構關系

按上述算法在某一矩形區(qū)域內(nèi)生成顆粒堆積體，再從中切割出一個邊長為1m的正方形試樣，試樣的詳細切割過程見文獻［12］。試樣內(nèi)的顆粒數(shù)目為2795，顆粒的最小半徑為 0.00625m，最大半徑為0.0125m。為了模擬顆粒之間的接觸和滑動，在任意兩個相互接觸的顆粒之間規(guī)定了界面。在本文中，將顆粒和界面都離散為正方形單元構成的集合體。圖2中深灰色的單元為界面單元，而淺灰色的單元為顆粒單元。若包括孔隙在內(nèi)，則本文模型中有500×500個單元。模型中顆粒和界面單元總數(shù)分別為174937及29353。由此可以計算出孔隙度約為（500×500-174937-29353）/（500×500）=0.183。

采用了5個計算方案。方案1至5的模型的中心孔洞半徑R不同，其它參數(shù)均相同。方案1至5的R分別為0.075m、0.125m、0.175m、0.225m 及 0.275m。模型四周的水平和垂直方向所施加的壓應力均為3MPa。監(jiān)測了方案1和5中一行單元的剪切應變增量、最大主應力及最小主應力的分布規(guī)律。這一行單元是通過巷道中心的水平線上的單元，且位于巷道的左側。圖2（a）中的擬開挖的巷道與監(jiān)測單元的位置是以方案5為例，各方案中的巷道半徑和監(jiān)測單元的數(shù)量有所不同。

計算步驟如下：首先，在給定的本構關系、邊界條件及加載條件下，對未開挖的模型進行計算，阻尼由FLAC自行施加，直至達到靜力平衡狀態(tài)。在本文中，當最大失衡力小于0.05N時，停止計算。此時，認為模型已經(jīng)達到了靜力平衡狀態(tài)。然后，利用編寫的FISH函數(shù)［13］刪除孔洞內(nèi)的單元，即開挖巷道。由于開挖卸荷產(chǎn)生了失衡力，因此模型不再處于平衡狀態(tài)。最后，對開挖后的模型進行計算，阻尼由 FLAC自行施加，一直計算到4萬步為止。從失衡力-時步曲線可以發(fā)現(xiàn)，此時，各方案均達到了靜力平衡狀態(tài)。本文僅給出了此時的計算結果。

圖2 顆粒堆積體模型（a）及局部放大圖（b），圖2（a）中的巷道尺寸及監(jiān)測單元是以方案5為例Fig.2 The model of granular material（a）and a close-up figure（b），where the tunnel size and the position of monitored elements are for Scenario 5

在本文中，顆粒體被視為各向同性的彈性材料，界面在破壞之前也被視為各向同性的彈性材料，而界面在破壞之后被視為莫爾-庫侖材料。彈性模量取為26.52GPa，泊松比取為0.21，內(nèi)聚力取為2MPa，內(nèi)摩擦角取為30°，抗拉強度取為1MPa。

3 結果分析及討論

3.1 剪切應變增量的分布規(guī)律

圖3給出了方案1至5的剪切應變增量高值區(qū)分布的等值線圖。圖3中的黑色區(qū)域為剪切應變增量的高值區(qū)。為使剪切應變增量的分布顯示得更清晰，圖3（a-e）中均未顯示大于3×10-4的剪切應變增量。由圖3可以發(fā)現(xiàn)，剪切應變增量的高值區(qū)主要位于巷道表面附近，且呈圓環(huán)形分布。在方案1至5中，圓環(huán)的圈數(shù)大致分別為1、2、3、5及8。由此可以發(fā)現(xiàn)，圓環(huán)的圈數(shù)隨孔洞半徑的增大而呈先慢后快的增長趨勢。計算表明，方案1至5的剪切應變增量的最大 值 分 別 為 0.8466 ×10-3、1.3461 ×10-3、1.8376 ×10-3、3.0219 ×10-3及 5.4776 ×10-3。由此可以發(fā)現(xiàn)，剪切應變增量的最大值隨孔洞半徑的增大而先慢后快地增大。這些剪切應變增量的最大值均出現(xiàn)在巷道表面附近。圖4給出了方案1至5的剪切應變增量的最大值位置的局部放大圖。在圖4中，顏色越深表示剪切應變增量的值越高，可以發(fā)現(xiàn)，剪切應變增量的最大值位置的單元（箭頭所指單元）發(fā)生了較為明顯的剪切變形。為什么是巷道表面附近的這些位置的單元具有最高的剪切應變增量？或許是由于這些單元附近有較大的孔隙且它們周圍的顆粒較為稀疏所致，因此這些單元易發(fā)生剪切錯動而產(chǎn)生較高的剪切應變增量。

圖5給出了方案1和5中監(jiān)測單元的剪切應變增量的分布規(guī)律。監(jiān)測的單元被重新進行了編號，第一個單元為離模型左邊界最近的單元，圖5中的橫坐標表示單元編號（不包括孔隙在內(nèi)）。由圖5可以發(fā)現(xiàn)，剪切應變增量的分布曲線具有多個峰。方案1和5中巷道左側一定范圍內(nèi)的單元的剪切應變增量的值都較高。方案5的剪切應變增量的值及高值區(qū)的范圍都遠大于方案 1。此外，還可以發(fā)現(xiàn)，方案 1和5中峰的個數(shù)分別與圖3（a）和（e）中圓環(huán)的圈數(shù)大致相同。

圖3 不同方案的剪切應變增量分布Fig.3 Shear strain increment distribution in different scenarios

圖4 不同方案的剪切應變增量的最大值位置的局部放大圖Fig.4 C lose-up figures of the position with maxim um shear strain increment in different scenarios

圖5 方案1和5的監(jiān)測單元剪切應變增量分布Fig.5 Shear strain increment distributions of monitored elements in scenarios 1 and 5

3.2 最小主應力的分布規(guī)律

圖6給出了各方案的最小主應力的高值區(qū)分布的等值線圖。為了清晰地顯示最小主應力的分布規(guī)律，圖6（a-e）中只顯示了-10MPa至0之間的最小主應力。FLAC中的最小主應力一般為負值，代表很高的壓應力。在圖6中，顏色越深表示壓應力的值越大。圖6與圖3中剪切應變增量高值區(qū)的分布規(guī)律較為類似。計算表明，方案1至5的最小主應力的最小值（即最大的 壓 應 力）分 別 為 -15.798MPa、-21.916MPa、-18.787MPa、-28.950MPa及 -28.443MPa。由此可以發(fā)現(xiàn)，各方案的最大的壓應力約為模型四周所施加的壓應力（3MPa）的5～10倍。

圖7給出了方案1和5中監(jiān)測單元的最小主應力的分布規(guī)律。由圖7可以發(fā)現(xiàn)，最小主應力的分布曲線也出現(xiàn)了多個峰。方案1和5中巷道左側一定范圍內(nèi)的單元的最小主應力的值（壓應力）都較大，而方案5的最小主應力的值及高值區(qū)的范圍都遠大于方案1。此外，圖7中具有較大的壓應力（絕對值大的負值）的位置大部分都與圖5中剪切應變增量的峰的位置相同。

3.3 最大主應力的分布規(guī)律

圖6 不同方案的最小主應力分布Fig.6 Minim um principle stress distribution in different scenarios

圖7 方案1和5的監(jiān)測單元的最小主應力的分布Fig.7 Minim um principle stress distributions of monitored elem en ts in scenarios 1 and 5

圖8給出了不同方案的最大主應力的高值區(qū)分布的等值線圖。在圖8中，為了更清晰地顯示最大主應力的分布規(guī)律，只顯示了0至3MPa之間的最大主基本上等于或大于模型四周所施加的壓應力（3MPa）。應力。圖8中的黑色區(qū)域為最大主應力的高值區(qū)。FLAC中的最大主應力通常為正值，表示較高的拉應力。由圖8（a-e）可以發(fā)現(xiàn)，最大主應力呈輻射狀分布，且隨著孔洞半徑的增大，最大主應力的高值區(qū)由巷道表面附近逐漸延伸到圍巖內(nèi)部。方案1至5中最大主應力的高值區(qū)延伸的長度分別約為15、28、42、60及79個單元。由此可以發(fā)現(xiàn)，最大主應力的高值區(qū)延伸的長度也隨孔洞半徑的增大而呈先慢后快的增長趨勢。由圖8及圖3可以發(fā)現(xiàn)，各方案中最大主應力的高值區(qū)延伸的長度與剪切應變增量的高值區(qū)分布的范圍大致相同。計算表明，方案1至5的最大主應力的最大值（最大的拉應力）分別為2.841MPa、 3.233MPa、 3.541MPa、 4.536MPa 及4.339MPa。由此可以發(fā)現(xiàn)，各方案的最大的拉應力

圖8 不同方案的最大主應力分布Fig.8 Maxim um principle stress distribution in different scenarios

圖9給出了方案1和5中監(jiān)測單元的最大主應力的分布規(guī)律。由圖9可以發(fā)現(xiàn)，最大主應力的分布曲線也具有多個峰。方案1和5中巷道左側一定范圍內(nèi)的單元的最大主應力的值（拉應力）都較大，而方案5的最大主應力的值及高值區(qū)的范圍都遠大于方案1。兩個方案中高值區(qū)的范圍分別與圖8（a）和（e）中拉應力的延伸的長度大致相同。為了揭示剪切應變增量與差應力（最大主應力與最小主應力之差）之間的關系，圖10同時給出了方案5中監(jiān)測單元的剪切應變增量與差應力的分布規(guī)律?？梢园l(fā)現(xiàn)，二者的峰值的位置基本相同。因而，具有較高的差應力的位置，也具有較高的剪切應變增量。

圖9 方案1和5的監(jiān)測單元的最大主應力的分布Fig.9 Maxim um principle stress distributios of monitored elements in scenarios 1 and 5

圖10 方案5的監(jiān)測單元的差應力與剪切應變增量的分布Fig.10 The distributions of the shear strain increm ent and the differential stress of monitored elem ents in scenario 5

3.4 簡短的討論

對于 Berea砂巖，其孔隙度在 17％ ～25％［14］，本文模型的孔隙度為18.3％，與其下限更接近。通常，砂巖的平均顆粒半徑越大，則孔隙度越大。當砂巖的孔隙度為25％時，其平均顆粒半徑可以達到0.25mm［14］，而 在 本 文 中，顆 粒 的 平 均 半 徑 為0.009375m，將近1cm。因此，這是顆粒實際尺寸的40倍以上。不僅本文作者，其他作者在研究尺度稍大的問題中，顆粒尺寸的選取都普遍高于顆粒的真實尺寸。例如，文獻［15］中顆粒半徑為 0.2m；文獻［16］中顆粒半徑為0.01m，這與本文接近；文獻［17］中顆粒半徑約為1.5cm。和上述三篇文獻相比，本文中顆粒的尺寸是最小的。限于本文作者所使用的微機的硬件水平，目前還無法在1m×1m的模型中，考慮顆粒的真實尺寸。盡管如此，并不妨礙對洞室尺寸效應問題的理解，這是本文工作的重點。

洞室尺寸效應問題的研究，在大規(guī)模巖石工程不斷向深部發(fā)展，以及需要越來越大的地下洞室被開挖的今天，顯得愈加重要和緊迫。在高地應力條件下開挖大斷面的地下洞室時，如果僅在設計階段基于均勻、連續(xù)介質力學方法進行分析和計算，可能會低估巖石工程建設及運行階段發(fā)生巖爆等地質災害的風險。因此，應該尋找更加行之有效的研究手段和技術，巖石的細觀結構不應該被忽視。

4 結論

（1）剪切應變增量與最小主應力的高值區(qū)（壓應力區(qū)）都呈圓環(huán)形，主要分布于巷道表面附近，且隨著孔洞半徑的增大，圓環(huán)的圈數(shù)及剪切應變增量的最大值都呈先慢后快的增長趨勢。最大主應力的高值區(qū)（拉應力區(qū)）呈輻射狀，其延伸范圍隨著孔洞半徑的增大，也呈先慢后快的增長趨勢。模型中最大的拉應力接近于模型邊界上所施加的壓應力，而最大的壓應力約為所施加的壓應力的5～10倍。

（2）巷道左側圍巖中監(jiān)測單元的剪切應變增量、最小主應力及最大主應力的分布曲線都具有多個峰，這表明它們的分布是高度不均勻的。在巷道表面附近，如果有些位置的單元周圍的孔隙較大且顆粒較為稀疏，則這些單元可能更易發(fā)生剪切變形而產(chǎn)生較高的剪切應變增量。具有較高的差應力的位置與具有較高的剪切應變增量的位置具有很好的一致性。

［1］Ba?ant Z P，Chen E P.結構破壞的尺度律［J］.力學進展，1999，29（3）：383-433.Ba?ant Z P，Chen E P.Scaling of structural failure［J］.Advances in Mechanics，1999，29（3）：383-433.

［2］Hudson J A，Crouch S L，F(xiàn)airhurst C.Soft，stiff and servo-controlled testing machine：a review with reference to rock failure［J］.Eng Geol，1972，6（3）：155-189.

［3］Jansen D C，Shah S P.Effect of length on compressive strain softening of concrete［J］.J Eng Mech，ASCE，1997，123（1）：25-35.

［4］Van Mier J G M，Shah S P，Arnaud M，et al.Strain-softening of concrete in uniaxial compression-report of the round robin test carried out by RILEM TC 148-SSC［J］.Mat Struct，1997，30（198）：195-209.

［5］王學濱，潘一山，宋維源.巖石試件尺寸效應的塑性剪切應變梯度模型［J］.巖土工程學報，2001，23（6）：711-714.WANG Xuebin， PAN Yishan， SONG Weiyuan.The model of plastic shear strain gradient on size effect in uniaxial compression of rock specimen［J］.Chinese Journal of Geotechnical Engineering，2001，23（6）：711-713.

［6］WANG X B.Prediction of height effect， plastic deformation and fracture energy for high-strength concrete by linear shear softening constitutive relation based on energy conservation method［J］.Magazine of Concrete Research，2007，59（5）：341-350.

［7］Carter B C.Size and stress gradient effects on fracture around cavities［J］.Rock Mech and Rock Eng，1992，25（3）：167-186.

［8］Carter B C，Lajtai E Z，Yuan Y.Tensile fracture from circular cavities loaded in compression［J］.Int J Fracture，1992，57：211-236.

［9］孟陸波，趙建壯.巖爆的洞室效應［J］.中國地質災害與防治學報，2008，19（2）：59-62.MENG Lubo，ZHAO Jianzhuang.The chamber effects of rock burst［J］.The Chinese Journal of Geological Hazard and Control，2008，19（2）：59-62.

［10］孟陸波，李天斌，周春宏，等.某水電站洞室?guī)r爆的尺寸效應［J］.水力發(fā)電，2008，34（1）：29-31.MENG Lubo，LITianbin，ZHOU Chunhong，et al.Rock burst′s dimensional effect on the tunnels of a hydropower station［J］.Water Power，2008，34（1）：29-31.

［11］王學濱，潘一山，陶帥.不同尺寸的圓形隧洞剪切應變局部化過程模擬［J］.中國地質災害與防治學報，2009，20（4）：101-108.WANG Xuebin，PAN Yishan，TAO Shuai.Numerical simulation of shear strain localization processes along circular tunnels with different diameters［J］.The Chinese Journal of Geological Hazard and Control，2009，20（4）：101-108.

［12］伍小林，王學濱，潘一山.單向壓縮條件下圓形顆粒體含孔洞試樣的力學模擬［J］.水利水運工程學報，2010，（3）：40-44.WU Xiaolin，WANG Xuebin，PAN Yishan.Numerical simulation of mechanical behavior of a specimen composed of circular granular material［J］.Hydro-Science and Engineering，2010，（3）：40-44.

［13］王學濱，王瑋，潘一山.孔隙壓力條件下圓形巷道圍巖的應變局部化數(shù)值模擬［J］.煤炭學報，2010，35（5）：723-728.WANG Xuebin，WANG Wei，PAN Yishan.Numerical simulation of the strain localization of the surrounding rock of a circular tunnel at different pore pressures［J］.Journal of China Coal Society，2010，35（5）：723-728.

［14］Haimson B C.Borehole breakouts in berea sandstone reveal a new fracture mechanism［J］.Pure and Applied Geophysics，2003，160：813-831.

［15］王濤，呂慶，李楊，等.顆粒離散元方法中接觸模型的開發(fā)［J］.巖石力學與工程學報，2009，28（增2）：4040-4045.WANG Tao，LU Qing，LI Yang，et al.Development of contact model in particle discrete element method［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（Supp.2）：4040-4045.

［16］周健，鄧益兵，賈敏才，等.基于顆粒單元接觸的二維離散 -連續(xù)耦合分析方法［J］.巖土工程學報，2010，32（10）：1479-1484.ZHOU Jian，DENG Yibing，JIA Mincai，et al.Coupling method of two-dimensional discontinuum -continuum based on contact between particle and element［J］.Chinese Journal of Geotechnical Engineering，2010，32（10）：1479-1484.

［17］Hazzard J F， Young R P.Moment tensors and m icromechanical models［J］.Tectonophysics， 2002，356：181-197.