圓形巷道圍巖層裂或板裂化的等效連續(xù)介質(zhì)模型 及側(cè)壓系數(shù)的影響

王學(xué)濱，伍小林，潘一山

（遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院，遼寧 阜新 123000）

1 引 言

巷道圍巖的層裂現(xiàn)象是指巷道周圍的巖石在壓應(yīng)力的作用下，巖石被劈裂成薄板的一種奇特現(xiàn) 象[1－4]，這種現(xiàn)象在礦柱中[5]也十分常見。層裂現(xiàn)象與巖爆、沖擊地壓也密切相關(guān)。在巖爆和沖擊地壓的機制分析方面，有些理論模型認為巷道附近的巖石中包含著等間距的裂紋或薄板（薄板位于兩裂紋之間）[1－4]，據(jù)此可對模型進行某種簡化和處理。在巖爆后的巷道圍巖中，經(jīng)常能發(fā)現(xiàn)各種形狀及尺寸的凹坑，常稱這些凹坑為巖爆坑，其中以V 形或三角形的巖爆坑最為常見。文獻[1]認為，V 形巖爆坑形成的主要過程包括4 個階段：損傷階段、擴容階段、劈裂成板的剝離階段及穩(wěn)定階段。在第3 階段，斷裂過程區(qū)的發(fā)展導(dǎo)致了薄板的形成，它們通過剪切、劈裂和屈曲產(chǎn)生，其厚度在1～5 cm 之間。文獻[1]還認為，這些薄板之間的凝聚力已損失高達100%，如果不對巷道施加及時的支護，這些像一疊卡片似的薄板易于在重力作用下從開挖邊界剝落下來，從而形成呈V 形的巖爆坑。文獻[2]認為，緊挨著V 形坑的尖端的巖體中也包含著若干薄板，但凝聚力的損失并不完全，因而它們?nèi)匀徽吃谝黄稹?/p>

在巖石破壞過程的數(shù)值模擬方面，通常都采用非均質(zhì)模型，目前，已經(jīng)取得了一些重要的研究進展[6－7]，在這些模型中，巖石均被處理成連續(xù)介質(zhì)，只不過單元的力學(xué)參數(shù)是隨機分布的。但巖石不僅是一種非均勻介質(zhì)，也是一種非連續(xù)介質(zhì)，巖石內(nèi)部存在各種形式的孔洞和微裂紋等微觀缺陷，僅將巖石視為非均勻介質(zhì)，并不能全面地反映巖石的物理本質(zhì)。以砂巖為例，砂巖屬于典型的顆粒體材料，它們與散體不同，是由圓形或橢圓形的砂粒通過膠結(jié)物的膠結(jié)而形成的巖石[8]。這類巖石通常都具有較高的孔隙度，以砂巖為例，孔隙度可達30%以上，一般為10%～20%。

本文將巖石視為顆粒堆積體，其中包含眾多的顆粒和一定的孔隙。如果顆粒之間發(fā)生接觸，則在顆粒之間設(shè)置界面。另外，將顆粒和界面均視為連續(xù)介質(zhì)，以區(qū)別于離散元方法。本文提出的方法在FLAC（連續(xù)介質(zhì)快速拉格朗日分析）中得以實現(xiàn)。本文的特點在于采用處理連續(xù)介質(zhì)的方法處理包含一定孔隙的等效連續(xù)介質(zhì)。

2 方法介紹

本文方法包括6 個步驟：（1）生成顆粒堆積體；（2）切割模型；（3）確定界面的尺寸和位置；（4）剖分單元，修正界面單元；（5）賦予單元本構(gòu)關(guān)系；（6）加載、計算、結(jié)果保存。

2.1 生成顆粒堆積體

在本文中，巖石被視為眾多圓形顆粒在重力作用下堆積而成的顆粒堆積體。在某一范圍A 內(nèi)隨機生成顆粒堆積體模型。顆粒為圓形，其堆積過程從下到上進行。為了模擬單個顆粒在重力作用下的自然堆積過程，需要保證任一顆粒在重力作用下能靜力平衡，為此任一顆粒至少要與其下方或兩側(cè)的另外兩個顆粒接觸，此外，擬堆積的任一顆粒的中心應(yīng)位于已堆積的兩個顆粒的中心之間。在顆粒堆積過程中，需要自動記錄和修改當(dāng)前顆粒堆積體的最上一層（稱之為表層）的顆粒，例如，當(dāng)任一新的顆粒被堆積在顆粒堆積體上時，需要對當(dāng)前表層各個顆粒的記錄進行修改（刪除舊記錄和增加新記錄）。在堆積一個新的顆粒時，將優(yōu)先堆積在顆粒堆積體表層的最低凹處附近，目的是使顆粒堆積體的表層盡可能平坦。

2.2 切割模型

在A 區(qū)域內(nèi)切割一個正方形區(qū)域B。選擇中心位于正方形區(qū)域B 之內(nèi)的顆粒的集合體作為試樣或計算模型。在切割試樣時，不允許一個顆粒被切成兩半。盡管有些顆粒的中心位于區(qū)域B 之內(nèi)，但如果它們與B 內(nèi)的其他顆粒僅有一個接觸點，則需要刪除它們。在本文中，區(qū)域B 的尺寸為1 m×1 m，區(qū)域B 內(nèi)的顆粒數(shù)目為2 795 個。顆粒的最大半徑、最小半徑分別為0.012 5 m、0.006 25 m。假設(shè)顆粒的半徑服從均勻分布，對顆粒半徑的統(tǒng)計結(jié)果見圖1，橫坐標(biāo)為在顆粒半徑的范圍內(nèi)將其等分成20 個區(qū)間，縱坐標(biāo)為落在各半徑區(qū)間內(nèi)的顆粒數(shù)目占顆粒總數(shù)的百分比。

圖1 顆粒半徑分布圖 Fig.1 Distribution of granular radius

2.3 確定界面的尺寸和位置

如果顆粒堆積體中的兩個顆粒發(fā)生接觸或者二者之間的距離非常近，則需要在二者之間規(guī)定界面。設(shè)置界面的目的是模擬顆粒之間的滑動和接觸。界面尺寸的確定原則如圖2(a)所示。在圖2(a)中，圓O1和O2代表2 個圓形顆粒，其半徑分別為r1和r2，假定r1＜r2。點O1、O2之間的距離設(shè)為d。假定這2個顆粒之間存在一個界面，用矩形ABCD 表示，界面的中心點為O，界面的長和寬分別為l 和b，界面的長度方向是指垂直于O1O2的方向，而界面的寬度方向是指平行于O1O2的方向。界面的位置（由O點確定）及尺寸（由l 和b 確定）的確定需要根據(jù)已知的參數(shù)，例如點O1的坐標(biāo)（x1，y1）、點O2的坐標(biāo)（x2，y2）、r1、r2、d 及將來剖分的正方形單元的尺寸（邊長）R。

在本文中，當(dāng)2 個顆粒之間的距離d0小于一個正方形單元的尺寸R 時，假定這2 個顆粒發(fā)生接觸，此時需要在二者之間設(shè)置界面，即

圖2 兩個顆粒及它們之間的界面 Fig.2 Two granular and an interface between them

由圖2(a)可以發(fā)現(xiàn)，界面的寬度b 應(yīng)大于R。b的值的選取需要憑借經(jīng)驗，可以先給b 一個稍大一些的值，例如取R 的2 倍，在將顆粒和界面剖分成單元之后，如果必要可再對顆粒和界面單元進行適當(dāng)?shù)男拚?，以保證界面不致太厚，也不致太薄（致使無法將兩個顆粒的單元完全隔開）。

假定界面的2 個角點A、B 位于較大的顆粒的圓周上。設(shè)界面的中心點O 的坐標(biāo)為（x，y）。由于O 點在線段O1O2上，因此，只需確定x 即可。這里假定，x 與O1點的橫坐標(biāo)x1、O2點的橫坐標(biāo)x2有下列簡單的關(guān)系：

在此，可以對式（2）的依據(jù)進行簡略的討論。由式（2）可以發(fā)現(xiàn)，當(dāng)r1→0 時，x→x1。也就是說，當(dāng)小顆粒很小時，界面ABCD 的中心點O 向小顆粒的中心趨近。這么做的目的是為了確保界面能將2個相差較大的顆粒聯(lián)接起來（見圖2（b））。否則如果r1→0 時，讓x→x2（不再遵循式（2）），則完全可能出現(xiàn)圖2(c)的情形。此時，界面沒有起到聯(lián)結(jié)2 個顆粒的作用，這樣，這2 個被判為相互接觸的顆粒由于界面位置設(shè)置不當(dāng)將不發(fā)生任何作用。

假設(shè)線段AB 與線段O1O2相交于E 點，線段EO2的長度為a，線段OO2的長度為c，則c 為

在ΔAEO2中，根據(jù)幾何關(guān)系，可得

這樣，利用式（3）、（4）可以得到界面長度l的表達式。界面的位置和尺寸確定之后，需要將界面和顆粒都剖分成正方形單元，其余的部分是孔隙。

2.4 剖分單元，修正界面單元

為了采用FLAC 模擬等效連續(xù)介質(zhì)，顆粒O1、O2和界面ABCD 都將被離散為正方形單元。因此，實際上，圖2 中的矩形ABCD 和顆粒O1、O2僅是離散后的界面和顆粒的輪廓。位于矩形ABCD 內(nèi)的正方形單元的集合構(gòu)成了真實的界面，這些單元稱之為界面單元，而位于圓內(nèi)的單元稱之為顆粒單元。

如果界面單元不能將2 個顆粒的單元完全隔開，則2 個顆粒的一些單元可能共用一些節(jié)點，由于在計算中節(jié)點無法分離，這對于模擬2 個顆粒沿界面的相對運動是不利的，因為這將使顆粒無法做相對運動或相對運動受到阻礙。因此，適當(dāng)修正是必要的。如果界面單元顯得過厚，則要將那些多余的界面變成顆粒單元或孔隙；如果界面顯得過薄，則需要適當(dāng)將顆粒單元或孔隙當(dāng)成界面單元。

將顆粒和界面剖分成單元，并經(jīng)過適當(dāng)修正之后，經(jīng)過統(tǒng)計可以得到下列數(shù)據(jù)：計算模型水平及垂直方向上的最大單元數(shù)約為500，顆粒單元的總數(shù)為174 937，界面單元的總數(shù)為29 353。據(jù)此，計算模型的孔隙度大約為0.182 8（（500×500－ 174 937－29 353）/（500×500））。對于Berea 砂巖，其孔隙度在17%～25%之間[9]，本文的孔隙度與其下限更接近。

2.5 賦予單元本構(gòu)關(guān)系

顆粒單元被視為各向同性的彈性材料，彈性模量取為26.52 GPa，泊松比取為0.21，而界面單元在破壞之后被視為摩爾-庫侖材料，界面單元的凝聚力取為2 MPa，內(nèi)摩擦角取為30°。在界面單元破壞之前，其本構(gòu)關(guān)系與顆粒單元一致。

2.6 加載、計算、結(jié)果保存

加載步驟如下：首先，給定本構(gòu)關(guān)系、邊界條件及模型四周的圍壓，對未開挖的顆粒堆積體模型進行計算，阻尼由FLAC 自行施加，直到達到靜力平衡狀態(tài)。在本文中，若最大失衡力小于0.05 N 時則停止計算，認為模型已經(jīng)達到了靜力平衡狀態(tài)。由于模型尺寸較小，因此，在計算中沒有考慮重力。

然后，為了模擬開挖，將巷道內(nèi)部的單元刪除。以模型的中心做1 個直徑為模型尺寸1/3 的圓。將那些位于圓內(nèi)的顆粒和界面單元刪除，以模擬巷道的一次性開挖。為了使巷道的表面看起來更圓一些，這里允許將顆粒和界面一分為二，僅刪除那些位于圓內(nèi)的單元。

圖3 顆粒堆積體模型 Fig.3 Model of circle granular material

開挖后的模型見圖3，從中可以看到兩種單元和孔隙。巷道開挖之后，模型的四周仍然受到壓應(yīng)力作用。由于開挖卸荷，因而在巷道附近的節(jié)點上將存在失衡力，模型將不再處于靜力平衡狀態(tài)。

最后，對開挖后的模型重新進行計算，阻尼由FLAC 自行施加。對于每個方案，從開挖開始之時開始計算，一直計算到4×104步。從失衡力-時步的演變規(guī)律來看，各個方案均達到了靜力平衡狀態(tài)。

本文采用了4 個計算方案（見表1）。方案1～4的模型四周的壓應(yīng)力不同，而其余參數(shù)全相同。

表1 計算方案 Table 1 Computational schemes

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

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

圖4 給出了不同方案的剪切應(yīng)變增量高值區(qū)分布的等值線圖，是從巷道開挖之后再計算4×104步的結(jié)果。在圖中，顏色越深代表單元的剪切應(yīng)變增量越大。本文中隨后的其他圖片也是這一時步時的結(jié)果，不再贅述。此時，開挖后的模型已經(jīng)重新達到平衡。由圖還可以發(fā)現(xiàn)，當(dāng)水平方向的應(yīng)力小于垂直方向的應(yīng)力時，在巷道兩側(cè)的圍巖中出現(xiàn)了較高的剪切應(yīng)變增量，見圖4(a)、(b)。與方案2 相比，方案1 的剪切應(yīng)變增量的高值區(qū)嚴(yán)格位于巷道的兩側(cè)，而在巷道的頂部及底部的圍巖中剪切應(yīng)變增量較低。當(dāng)模型的四周的壓應(yīng)力相等時（見圖4(c)）圍巖中的剪切應(yīng)變增量高值區(qū)的等值線圖呈圓環(huán)狀，位于巷道附近。圖4(d)的結(jié)果與圖4(a)剛好相反。在圖4(d)中，圍巖中的剪切應(yīng)變增量的高值區(qū)位于巷道頂部和底部，長條形的剪切應(yīng)變增量等值線基本上都在水平方向上。

計算表明，方案1～4 的最大剪切應(yīng)變增量分別為 2.412 1×10-3、1.891 3×10-3、1.837 6×10-3及 2.456 1×10-3。可以發(fā)現(xiàn)，方案1、4 的最大剪切應(yīng)變增量的值較高，而方案3 的值較小。另外，在圖4 中，為了能更清晰地顯示出長條形的剪切應(yīng)變增量等值線圖，未顯示大于2.5×10-4的剪切應(yīng)變增量。

圖4 不同方案的剪切應(yīng)變增量（當(dāng)時步=4×104） Fig.4 Contour maps of shear strain increment at different schemes (timesteps=4×104)

3.2 垂直方向應(yīng)力的分布規(guī)律

圖5 給出了4 方案的垂直方向應(yīng)力的高值區(qū)分布的等值線圖，該圖僅顯示了大于0 的垂直方向的應(yīng)力分布。在FLAC 中，以壓應(yīng)力為負，而以拉應(yīng)力為正。為了能更清晰地顯示出垂直方向的應(yīng)力為拉應(yīng)力的單元的拉應(yīng)力分布規(guī)律，在圖中未顯示大于1.5 MPa 的垂直方向應(yīng)力。在圖中，只有單元的垂直方向應(yīng)力是拉應(yīng)力才被顯示出來?？梢园l(fā)現(xiàn)，在圖5(a)～(c)中，垂直方向應(yīng)力為拉應(yīng)力的單元的分布規(guī)律具有一定類似性。隨著水平方向應(yīng)力的增加，垂直方向應(yīng)力為拉應(yīng)力單元的數(shù)目有增多的趨勢；垂直方向應(yīng)力的高值區(qū)在巷道的頂部及底部形成的V 形坑越來越明顯。圖5(d)的計算結(jié)果與圖5(a)～(c)差別較大，如果顯示方案4 的水平方向的應(yīng)力分布（若只顯示大于0 部分，且旋轉(zhuǎn)90°后再比較），就不應(yīng)與圖5(a)有大的差別。

圖5 不同方案的垂直方向應(yīng)力（當(dāng)時步=4×104） Fig.5 Vertical stress distributions at different schemes (timesteps=4×104)

由計算結(jié)果可以發(fā)現(xiàn)，方案1～4 大于0 的垂直方向應(yīng)力的最大值分別為2.89、2.62、3.03、3.03 MPa。這些值都與所施加的最大壓應(yīng)力值（3.00 MPa）相近。這說明，盡管在模型的四周施加的是壓應(yīng)力，但由于巖石被處理為等效連續(xù)介質(zhì)，由于顆粒之間的相對運動規(guī)律較為復(fù)雜，因此，在模型中誘發(fā)出了與所施加的壓應(yīng)力接近的拉應(yīng)力。眾所周知，通常巖石的抗拉強度是抗壓強度的1/10～1/4，因此，巖石在壓應(yīng)力作用下產(chǎn)生的如此高的拉應(yīng)力是易于將巖石拉壞的。筆者的數(shù)值計算經(jīng)驗表明，采用均質(zhì)或非均質(zhì)模型，在巖石內(nèi)部不可能產(chǎn)生如此大的拉應(yīng)力。盡管在壓縮載荷作用下，在巖石內(nèi)部將誘發(fā)出拉破壞在巖石力學(xué)界被廣泛承認，但試圖模擬出拉應(yīng)力或者超過一定值的拉應(yīng)力一直是困難的。本文的方法從等效連續(xù)介質(zhì)的角度出發(fā)，得到了不同于均質(zhì)和非均質(zhì)模型的計算結(jié)果。

3.3 最小主應(yīng)力的分布規(guī)律

圖6 給出了不同方案的最小主應(yīng)力的高值區(qū)分布的等值線圖。應(yīng)當(dāng)指出，為了使各方案的最小主應(yīng)力分布規(guī)律顯示得更清晰，在圖6 中僅顯示了那些小于0 的最小主應(yīng)力，而且未顯示超過10 MPa的壓應(yīng)力。在FLAC 中，最小主應(yīng)力一般為絕對值較大的負值，代表很高的壓應(yīng)力?？梢园l(fā)現(xiàn)，圖6的計算結(jié)果與圖4 中的剪切應(yīng)變增量分布的計算結(jié)果比較類似，圖6(a)、(d)的結(jié)果也剛好相反。計算表明，方案1～4 的最小主應(yīng)力的最小值（很高的壓應(yīng)力）分別為-17.01、-15.36、-18.79、-24.29 MPa，而最小主應(yīng)力的最大值（較高的拉應(yīng)力）分別為3.34、3.02、2.80、2.84 MPa。由此可以發(fā)現(xiàn)，模型中的最大壓應(yīng)力遠超過模型邊界上施加的壓應(yīng)力（最大值為3 MPa），可以高達8 倍，例如，在方案4 中，最小主應(yīng)力的最小值為-24.29 MPa，其大小是3.00 MPa 的8.1 倍。在如此高的壓應(yīng)力作用下，模型內(nèi)部將不可避免地發(fā)生損傷。模型中的最大壓應(yīng)力與邊界上所施加的壓應(yīng)力的比值必然與材料的微觀結(jié)構(gòu)有關(guān)。如果將本來是非連續(xù)介質(zhì)視為連續(xù)介質(zhì)，縱然在設(shè)計中考慮了一定的安全系數(shù)，也可能無法避免存在的安全隱患。

圖6 不同方案的最小主應(yīng)力（當(dāng)時步=4×104） Fig.6 The minor principal stress distributions at different schemes (timesteps=4×104)

另外，方案1～4 的最小主應(yīng)力的最大值在3 MPa附近，也近似等于在模型的邊界上所施加的最大的壓應(yīng)力的值。

3.4 最大主應(yīng)力的分布規(guī)律

圖7 給出了不同方案的最大主應(yīng)力的高值區(qū)分布的等值線圖。為了使各方案的最大主應(yīng)力分布規(guī)律顯示得更加清晰，圖中僅顯示了那些大于0 的最大主應(yīng)力，而且未顯示大于1.5 MPa 的最大主應(yīng)力。在FLAC 中，最大主應(yīng)力通常為正值，代表較高的拉應(yīng)力。可以發(fā)現(xiàn)，與最小主應(yīng)力的高值區(qū)呈環(huán)向分布不同（見圖6），最大主應(yīng)力的正高值區(qū)基本上都位于徑向。計算表明，方案1～4 的最大主應(yīng)力的最大值分別為4.69、4.25、3.54、4.35 MPa，而最小值分別為-4.15、-3.90、-3.59、-3.78 MPa。還可以發(fā)現(xiàn)，對于方案2、3，最大主應(yīng)力的正高值區(qū)都位于巷道附近，而對于方案1、4，最大主應(yīng)力的正高值區(qū)已經(jīng)深入到圍巖內(nèi)部。這是由于模型邊界上所受的2 個方向上的壓應(yīng)力之差較大的緣故。

圖7 不同方案的最大主應(yīng)力（當(dāng)時步=4×104） Fig.7 The maximum principal stress distributions at different schemes(timesteps=4×104)

3.5 一行單元應(yīng)力的分布規(guī)律

圖8給出了方案1～4的一行單元的水平和垂直方向應(yīng)力的分布規(guī)律。這一行單元選在巷道左側(cè)的圍巖中，而且取在通過巷道直徑的水平方向上。這一行單元僅是實體單元，即包括界面和顆粒單元，不包括孔隙。圖8 的橫坐標(biāo)為從左到右的單元數(shù)，第1 個單元表示離模型左邊界最近的單元，而第139個單元表示離巷道表面最遠的單元。因此，水平方向應(yīng)力相當(dāng)于徑向應(yīng)力rσ ，而垂直方向應(yīng)力相當(dāng)于環(huán)向應(yīng)力θσ 。由圖可以發(fā)現(xiàn)，rσ 、θσ 的分布規(guī)律和以往的均質(zhì)模型得到的解析解或數(shù)值解有重大的差別。在一些位置的θσ 的值非常高，因而，二者的分布曲線具有多個峰值，既不光滑，也不連續(xù)。

圖8 不同方案的巷道左側(cè)一行單元應(yīng)力分布 （當(dāng)時步=4×104） Fig.8 Stress distributions of a row of elements in surrounding rock at left side of a circular tunnel at different schemes (timesteps=4×104)

由圖8 可以發(fā)現(xiàn)，在巷道的表面，rσ 為0，而在模型的左邊界上，rσ 與所施加的應(yīng)力相等。這些結(jié)果的正確性是顯然的。在方案1～3 中，隨著Ph的增加，小于0 的θσ （表示壓）的峰值的絕對值下降，而小于0 的rσ 的峰值有上升的趨勢。方案4 的rσ 、θσ 的分布與方案1～3 差別較大，這是由于對于方案4 而言，在模型左、右邊界上施加的壓應(yīng)力的值大于在模型上、下邊界上施加的壓應(yīng)力值。如果顯示方案4 的位于拱頂?shù)囊涣袉卧乃胶痛怪狈较驊?yīng)力的分布，則不會與方案1 相比有大的差別。另外，由圖可以發(fā)現(xiàn)，在離巷道附近的圍巖中，有許多單元的rσ 為正值，這表明這些單元受到的徑向應(yīng)力為拉應(yīng)力，而且具有多個峰值；與方案2、3相比，方案1 中的拉應(yīng)力的范圍最大，而且值最高。拉應(yīng)力的分布具有多個峰值的特征應(yīng)該與圍巖發(fā)生層裂有關(guān)。

3.6 討論

（1）圖9(a)、(b)給出了2 個方向不等壓條件下，巷道圍巖的破壞形式[1－2，5]。這與本文方案1、2 的計算結(jié)果較為一致。與圖4(a)、(b)、(d)中的剪切應(yīng)變增量的高值區(qū)的分布相比，圖6(a)、(b)、(d)中的最小主應(yīng)力高值區(qū)的分布的結(jié)果與這些文獻給出的結(jié)果更為一致。環(huán)向的壓應(yīng)力高值區(qū)和剪切應(yīng)變增量高值區(qū)似乎具有一定的間距，該間距應(yīng)該與顆粒的平均尺寸有關(guān)，這一點需要進一步研究。

（2）眾所周知，在單向壓縮條件下，巖石易于劈裂，圖8 的結(jié)果表明，在巷道表面附近，一些單元的徑向應(yīng)力居然還是拉應(yīng)力，這更加劇了單元發(fā)生沿上述兩種高值區(qū)的劈裂。這應(yīng)該是層裂現(xiàn)象的根源，圖9(c)給出了一個巖石試件在一向壓應(yīng)力及另一向拉應(yīng)力作用下發(fā)生拉破壞的示意圖。

（3）在本文中，由于將顆粒和界面都剖分成了尺寸相同的單元，因此，采用的模擬等效連續(xù)介質(zhì)的方法計算量較大。一般認為，這種方法不適于工程實際問題的分析，更適于破壞機制的分析。然而，以離散元為基礎(chǔ)的顆粒流方法也常見被用于工程實際問題的分析之中。一些作者所選取的顆粒尺寸都普遍大于巖石顆粒的真實尺寸。以砂巖為例，其平均顆粒半徑可以達到0.25 mm[9]。文獻[10]在數(shù)值計算中選取的顆粒半徑為0.2 m；文獻[11]在數(shù)值計算中選取的顆粒半徑為0.01 m，這與本文接近；文獻

[12]在數(shù)值計算中選取的顆粒半徑約為1.5 cm。在計算中，采用大顆粒進行模擬，一方面是不得已而為之，另一方面或許可以在某種程度上近似巖石的非連續(xù)性，或者認為，大顆粒與節(jié)理巖體中的巖塊具有某種等效性。從這一方面講，利用比真實顆粒尺寸大得多的顆粒進行模擬，是可以用于工程尺度問題的定性分析的。此時，顆粒僅是巖體中完整部分的某種近似，而顆粒之間的接觸、錯動以及孔隙在某種程度上代表了巖體中的非連續(xù)性的某些側(cè)面。由于本文中顆粒的尺寸大于砂巖的實際顆粒的尺寸，因此，似乎應(yīng)該認為本文模型是節(jié)理巖體的某種等效或近似。盡管本文中顆粒的尺寸大于砂巖顆粒的真實尺寸，但從圖6(a)、(b)、(d)給出的最小主應(yīng)力高值區(qū)的分布和圖9(a)、(b)給出的破壞區(qū)域可以發(fā)現(xiàn)，二者的范圍大體相當(dāng)。

圖9 圍巖層裂的一些結(jié)果[1－2, 5]及巖樣受拉壓應(yīng)力 發(fā)生拉破壞的示意圖 Fig.9 Some results about the exfoliation phenomenon of the surrounding rock [1－2, 5] and schematic about the tensile failure of a rock specimen under combination of compressive and tensile stresses

（4）通過計算可以發(fā)現(xiàn)，模型中個別單元的應(yīng)力或應(yīng)變的值會很高，這與這些單元附近的微結(jié)構(gòu)的細節(jié)有密切的關(guān)系，例如孔隙度及顆粒的排列方式等，但其余的絕大部分單元的應(yīng)力都在某一范圍之內(nèi)，這應(yīng)該得到重點的關(guān)注。如果將所有的單元的應(yīng)力都顯示在一張圖片上，個別的高值就將主流的與高值相比的低值掩蓋了，因此，在圖4～7 種舍棄了個別的高值。

（5）限于本文篇幅，未能研究下列問題：界面單元強度的影響，孔隙存在及填充后的計算結(jié)果的差別，顆粒尺寸及排列方式的影響等。在文獻[13]中比較詳細地介紹了顆粒的堆積算法，研究了圓形洞室尺寸效應(yīng)問題，關(guān)于其他形狀洞室的研究也有必要開展。

4 結(jié) 論

（1）將巖石視為等效連續(xù)介質(zhì)后，所得的數(shù)值計算結(jié)果與均質(zhì)和非均質(zhì)模型的計算結(jié)果相比，具有重大的差異。本文的計算結(jié)果顯示，在模型中可以誘發(fā)出與所施加的壓應(yīng)力大小相當(dāng)?shù)睦瓚?yīng)力，模型內(nèi)部的壓應(yīng)力遠高于在模型邊界上所施加的壓應(yīng)力的值。上述結(jié)果無論從哪一方面講，過去的均質(zhì)和非均質(zhì)模型的計算結(jié)果都可能存在潛在的風(fēng)險。因此，非連續(xù)性有必要在設(shè)計中予以考慮。

（2）模型中的最小主應(yīng)力和剪切應(yīng)變增量的分布具有一定的類似性。當(dāng)水平方向和垂直方向的壓應(yīng)力不相等時，二者的高值區(qū)位于巷道的兩側(cè)或頂部和底部，而在靜水壓力條件下，最小主應(yīng)力和剪切應(yīng)變增量的高值區(qū)的分布近似呈環(huán)形，主要位于巷道的周圍。最小主應(yīng)力和剪切應(yīng)變增量的高值區(qū)分布基本平行，這些位置是將來裂紋出現(xiàn)的位置。圍巖層裂現(xiàn)象的原因被認為是環(huán)向的高壓應(yīng)力和徑向的高拉應(yīng)力的共同作用。

[1] MARTIN C D. The effect of cohesion loss and stress path on brittle rock strength[J]. Canadian Geotechnical Journal， 1997， 34(5)： 698－725.

[2] BA?ANT Z P， LIN F B， LIPPMANN H. Fracture energy release and size effect in borehole breakout[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 1993， 17(1)： 1－14.

[3] 馮濤， 潘長良. 洞室?guī)r爆機理的層裂屈曲模型[J]. 中國有色金屬學(xué)報， 2000， 10(2)： 287－290. FENG Tao， PAN Chang-liang. Lamination spallation buckling model for formation mechanism of rock burst[J]. The Chinese Journal of Nonferrous Metals， 2000， 10(2)： 287－290.

[4] 王敏強， 侯發(fā)亮. 板狀破壞的巖體巖爆判別的一種方法[J]. 巖土力學(xué)， 1993， 14(3)： 53－59. WANG Min-qiang， HOU Fa-liang. A method to determine rock burst in rock mass with plate-form damage[J]. Rock and Soil Mechanics， 1993， 14(3)： 53－59.

[5] KLERCK P A， SELLERS E J， OWEN D R J. Discrete fracture in quasi-brittle materials under compressive and tensile stress states[J]. Computer Methods in Applied Mechanics Engineering， 2004， 193： 3035－3056.

[6] TANG C A， THAM L G， LEE P K K， et al. Numerical studies of the influence of microstructure on rock failure in uniaxial compression-Part II： Constraint， slenderness and size effect[J]. International Journal of Rock Mechanics and Mining Sciences， 2000， 37(4)： 571－583.

[7] FANG Z， HARRISON J P. Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks[J]. International Journal of Rock Mechanics and Mining Sciences， 2002， 39(4)： 443－457.

[8] 張凱， 周輝， 胡大偉， 等. 彈塑性條件下巖土孔隙介質(zhì)有效應(yīng)力系數(shù)理論模型[J]. 巖土力學(xué)， 2010， 31(4)： 1035－1041. ZHANG Kai， ZHOU Hui， HU Da-wei， et al. Theoretical model of effective stress coefficient for rock/soil-like porous media under elastoplastic condition[J]. Rock and Soil Mechanic， 2010， 31(4)： 1035－1041.

[9] HAIMSON B C. Borehole breakouts in Berea sandstone reveal a new fracture mechanism[J]. Pure and Applied Geophysics， 2003， 160(5－6)： 813－831.

[10] 王濤， 呂慶， 李楊， 等. 顆粒離散元方法中接觸模型的開發(fā)[J]. 巖石力學(xué)與工程學(xué)報， 2009， 28(增刊2)： 4040－4045. WANG Tao， Lü 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.

[11] 周健， 鄧益兵， 賈敏才， 等. 基于顆粒單元接觸的二維離散-連續(xù)耦合分析方法[J]. 巖土工程學(xué)報， 2010， 32(10)： 1479－1484. ZHOU Jian， DENG Yi-bing， JIA Min-cai， 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.

[12] HAZZARD J F， YOUNG R P. Moment tensors and micromechanical models[J]. Tectonophysics， 2002， 356： 181－197.

[13] 伍小林， 王學(xué)濱， 潘一山. 含不同半徑孔洞的顆粒體模型的力學(xué)行為數(shù)值模擬[J]. 中國地質(zhì)災(zāi)害與防治學(xué)報， 2011， 22(1)： 107－113. WU Xiao-lin， WANG Xue-bin， PAN Yi-shan. Numerical simulation of mechanical behavior of the granular model with a hole with different radii[J]. The Chinese Journal of Geological Hazard and Control， 2011， 22(1)： 107－113.