基于两种颗粒体模型的巷道围岩应力应变分析

伍小林，王学滨，潘一山

（辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000）

0 引言

洞室围岩的破坏及失稳现象广泛存在于各种岩石工程中，例如水电洞室岩爆、油井崩裂坍塌及交通隧洞塌方，这些灾害严重地威胁着人员和设备的安全。因此，研究洞室围岩的破坏区分布、应力及应变分布具有重要的意义。

在数值分析方面，目前主要采用两种模型对洞室围岩的力学行为进行分析。其一是连续介质模型［1－6］，即将岩石视为连续介质，采用有限元或有限差分等数值方法进行分析。这种方法易于获得围岩的应力、应变、位移等力学量的时空分布规律，一般不直接考虑岩石中存在的小尺度的裂隙、缺陷、孔洞等因素，仅通过将岩石的力学性能降低，以间接地考虑这些因素的存在。其二是非连续介质模型，即将岩石视为各种非连续介质［7－9］，例如颗粒体或块体，采用离散元方法或不连续分析方法进行分析。这些方法一般并不容易获得围岩中的应力、应变的时空分布规律。文献［10－11］最近提出了一种非连续介质模型，此模型是以连续介质模型作为基础的。这种模型将岩石视为由颗粒和在颗粒之间的界面构成的颗粒体材料，将颗粒和界面都离散成尺寸相同的单元，进而可采用适于模拟连续介质力学行为的数值方法进行分析，可以方便地获得围岩中应力、应变、塑性区、声发射及能量释放率的分布及演变规律。在文献［10］中，研究了含圆形孔洞的颗粒体材料在单向压缩条件下的应力及应变的分布规律，在孔洞的顶部和底部发现了呈三角形的受拉区，这不同于连续介质模型的结果。在文献［11］中，研究了洞室的尺寸对围岩中的剪切应变增量、最大主应力及最小主应力的分布规律的影响。研究发现随着洞室半径的增大，呈圆环形的剪切应变增量与最小主应力的高值区的圈数、呈辐射状最大主应力的高值区的延伸范围都呈先慢后快的趋势增长。这说明洞室的破坏现象具有一定的尺寸效应，这可在定性的程度上解释文献［12－13］观察到的岩爆的尺寸效应。

本文在文献［11］的基础上，将颗粒体模型中的除了颗粒和界面之外的部分（即孔隙）用不同强度的材料进行填充，以等效模拟砂岩中强度弱的充填体（例如一些胶结物、杂质、碎屑等）及孔隙的综合效果，采用有限差分方法模拟了洞室围岩中的环向和径向应力以及塑性区、剪切应变增量的分布规律，着重研究了应力与传统模型（均质模型）的结果及与文献［11］中模型（含孔隙的颗粒体材料模型）的结果之间的差别，研究了应变高值区的分布规律，并将本文模型的结果与文献［11］的结果进行了对比。本文的研究对于进一步认清围岩中应力、应变的分布规律会产生积极的意义。

1 研究方法介绍

建立包含界面的颗粒堆积体模型包括下列4个步骤：生成颗粒堆积体、切割模型、确定界面的尺寸与位置及剖分单元。

首先，在某一区域A内随机生成圆形颗粒堆积体模型（详细过程参见［11］）。其基本原理是模拟圆形颗粒（例如，砂颗粒）在重力作用下的自然堆积过程。颗粒的半径服从均匀分布。颗粒的最大、最小半径分别为0.0125m及0.00625m。

其次，在上述区域A中切割出一个尺寸为1m×1m的正方形区域 B，详细的切割过程见文献［10］。在切割过程中，不允许一个颗粒被切成两半，这是由于颗粒往往强度高，不易在切割机上被切成两半，而它们往往完好地从母体中脱落下来。区域B内的颗粒数目为2795个。

然后，确定界面的尺寸与位置。在颗粒堆积体中，若两个颗粒发生接触或者它们之间的距离非常近，则需要在它们之间设置一个矩形界面。设置界面的目的是为了模拟颗粒之间的接触和滑动。

图1 考虑基体的圆形颗粒堆积体模型Fig.1 Themodelofcirclegranular materialconsideringmatrix

最后，将整个区域B剖分成若干正方形单元（图1（a））。图1（a）中的巷道是在本文第3节中开挖的。图1（b）是图1（a）的右下角的局部放大图。在图1（b）中，白色单元为基体单元，灰色单元为颗粒单元，黑色单元为界面单元。颗粒、界面及基体均被离散为正方形单元构成的集合体。位于圆形颗粒内的单元被称之为颗粒单元，位于矩形界面内的单元被称之为界面单元，除颗粒与界面单元之外的单元被称之为基体单元。本文采用众多正方形的颗粒单元模拟圆形的砂颗粒，采用基体单元模拟胶结物、孔隙及强度低的碎屑和杂质的综合效果。本文中正方形单元的边长为0.002m，模型共包含500×500个单元。

2 计算模型、方案及本构关系

本文共采用了两种模型进行计算：第一种为考虑基体的连续介质模型（方案1～7）；第二种为不考虑基体的连续介质模型（方案8），在该模型中将原来的基体单元删除，因此是非连续介质模型中的一种，与文献［10－11］中的模型相同。方案2～7的基体的强度与方案1的不同（分别为方案1的75％、50％、25％、10％、1％及0.1％），而其它参数均相同。方案8中没有基体，而代之以孔隙。在各个方案中，模型四周所施加的压应力均为3MPa。

颗粒单元被视为各向同性的线弹性材料，弹性模量取为26.5179GPa，泊松比取为0.21。界面与基体单元在破坏之前，其本构关系与颗粒单元一致，而在破坏之后被视为莫尔－库仑材料。界面单元的弹性模量和泊松比与颗粒单元的相同，内聚力取为2MPa，内摩擦角取为30°。方案1～7中所选取的基体单元的参数见表1。

表1 方案1～7的基体单元的参数Table1 Valuesforparametersofmatrix elementsinschemes1to7

计算步骤如下：首先，给定本构关系、边界条件及加载条件，为了使模型尽快达到静力平衡条件，采用了与文献［10－11］不同的处理方式：按平面应变条件初始化模型内部单元的各个方向的压应力值。在本文中，若最大失衡力小于0.05N，则认为模型已经达到了静力平衡条件。然后，利用编写的 FISH函数［5－6］删除巷道内部的单元。图1（a）中位于模型中部的圆形区域表示巷道，巷道的直径为0.25m。由于开挖卸荷产生了失衡力，因此模型不再处于静力平衡状态。最后，对开挖后的模型进行计算，阻尼由FLAC自行施加，一直计算到3万时步为止。从失衡力－时步曲线可以发现，此时，各方案均再次达到了静力平衡状态。本文给出的所有结果均为3万时步时的结果。

图2 方案1～8的破坏区的分布规律Fig.2 Failurezonedistributioninschemes1to8

3 计算结果及分析

3.1 破坏区的分布规律

图2给出了方案1～8的破坏区的分布规律。图2中的黑色区域为破坏区（方案1～3及8的黑色边框除外）。由图2（a）可以发现，方案1的破坏区的范围很小，仅出现在巷道表面上。

当基体强度分别降为方案1的75％（方案2，见图2（b））和50％（方案3，见图2（c））时，破坏区的范围增大，延伸到巷道表面附近。

当基体强度分别降为方案1的25％、10％、1％及0.1％（方案4～7，见图2（d－g））时，可以明显地观察到破坏区遍布于整个模型中。在离巷道一定范围内可以发现圆形颗粒（即圆形颗粒被暴露出来），这说明此范围内的界面与基体几乎全都发生了破坏，而在较远处，界面很少发生破坏，基体基本上都发生了破坏。此外，还可以发现，随着基体强度的降低，暴露出的圆形颗粒的范围在逐渐增大。

由方案8的计算结果（图2（h））可以发现，发生破坏的单元都是界面单元，在巷道表面附近发生破坏的单元较多，而在其它区域则较少。

3.2 一行单元应力的分布规律

3.2.1 一行单元水平方向应力的分布规律

图 3（a）和（b）给出了方案 1、3、4、7 及 8 中一行单元水平方向应力（与径向应力基本相等）的分布规律。为了便于分清不同方案的结果，图3中未给出其它方案的结果。这一行单元选在巷道左侧的围岩中，且取在通过巷道直径的水平方向上（图1（a））。图3的横坐标为上述一行单元的单元编号，第1号单元为离模型左边界最近的单元，第188号单元为离巷道表面最近的单元。在 FLAC中，以压应力为负值，而以拉应力为正值。

由图3可以发现，方案1的结果是光滑且连续的。在巷道的表面，水平方向的应力值为零，而在模型的左边界上，其值与所施加的应力值（3MPa）相等。这与传统的连续介质模型所得的结果［14］是类似的。

当基体强度降为方案1的50％（方案3）时，其结果呈现出一定的波动性。在第127号单元的左侧，方案3的结果在方案1的附近波动，而在其右侧，水平方向应力的绝对值几乎都比方案1的低。

当基体强度降为方案1的25％（方案4）时，其结果与方案3的有一定的类似性，但它的波动幅度更大。特别需要指出的是，在第158号到183号单元（共计25个单元）范围内竟然出现了拉应力区（水平方向应力大于零）。此拉应力区的范围与图2（d）中暴露出的圆形颗粒的范围大致相同。

当基体强度降为方案1的0.1％（方案7）时，其结果的波动幅度与方案4的相比进一步增大。方案7的结果的拉应力区的范围是从第136到188号单元（共计52个单元），是方案4的两倍多，且与图2（g）中暴露出的圆形颗粒的范围也大致相同。

图3 方案 1、3、4、7及8的一行单元水平方向应力的分布规律Fig.3 Horizontalstressdistributionofa rowofelementsinschemes1，3，4，7and8

方案8的结果与方案7的具有一定的类似性，但是经过仔细观察可以发现与方案3和4的也具有一定的类似性，水平方向应力绝对值较低的位置大致都相同。由此推测，这是由于细观结构（颗粒的排列方式）在其中起了重要的作用。

3.2.2 一行单元垂直方向应力的分布规律

图 4（a）和（b）给出了方案 1、3、4、7 及 8 中一行单元垂直方向应力（与环向应力基本相等）的分布规律。图4中横坐标的含义同第3.2.1节所述。由图4可以发现，方案1的结果是光滑的，在巷道表面附近出现一个峰值，离峰值越远其绝对值越低，直到趋向于一个定值（3.1MPa）。

当基体强度降为方案1的50％（方案3）时，其结果在方案1附近波动。

当基体强度分别降为方案1的25％（方案4）和0.1％（方案7）时，方案4和7的结果是比较类似的，其波动幅度都比方案3的大。

由方案8的结果可以发现，其波动幅度比方案3、4及7的都大，且垂直方向应力的绝对值普遍都低。很多单元的垂直方向应力的值为零，因为这些单元对应于孔隙，而离孔隙较近单元的垂直方向应力绝对值也都较低。此外，还可以发现，这些孔隙的位置与方案4和7的垂直方向应力的绝对值较低的位置大致相同。这是由于基体的强度低（其承载能力自然低），因而应力转移到了其附近的界面或颗粒单元上。

图4 方案1、3、4、7及8的一行单元垂直方向应力的分布规律Fig.4 Verticalstressdistributionofa rowofelementsinschemes1，3，4，7and8

3.3 剪切应变增量高值区的分布规律

图5给出了方案1～8的剪切应变增量高值区的分布规律。图5中颜色越深表示剪切应变增量的值越高（方案1～4的黑色边框除外）。为了更清晰地显示出剪切应变增量高值区的分布规律，图5中未显示大于5×10－4的剪切应变增量。由图5（a）可以发现，方案1的结果呈圆环形分布，且离巷道表面越近其值越高。这与传统的连续介质模型所得的结果［14］类似。

由图5（b－g）可以发现，随着基体强度的降低，方案2～7的剪切应变增量高值区的范围增大。方案2和3的结果的最大值都比方案1的高，且其结果已开始变得斑驳，这说明剪切应变增量在环向已经变得不均匀。这种现象出现的原因是由于基体的强度低，容易发生破坏，剪切应变增量会较高，因而在模型中出现了局部的高值区（即斑驳的现象）。方案5～7的剪切应变增量高值区已延伸到了模型的边界上。在方案4～7中，可以观察到类似于塑性力学中的滑移线（对数螺线），这些滑移线相互交织形成了滑移线网。这种现象出现的原因是由于基体非常弱，会使高的剪切应变增量集中于此，而颗粒始终作为线弹性材料，其剪切应变增量的值一直都较低。斑驳现象的出现实际上是滑移线网的初步发展。滑移线网的出现应该与由坚硬的弹性颗粒单元和软弱的基体单元构成的材料的非均质性有关。

由方案8的结果［11］（图5（h））可以发现，多个环向的应变集中区出现在巷道表面附近，这与方案1～7的结果有显著的差异。

图5 方案1～8剪切应变增量高值区的分布规律Fig.5 Shearstrainincrementdistribution inschemes1to8

4 结论

（1）当基体强度参数取为与界面的相等时，计算出的结果与传统的连续介质模型所得的结果类似，破坏区仅发生在巷道的表面上。随着基体强度参数的降低，巷道围岩中的径向和环向应力呈现的波动幅度越来越大；破坏区的范围逐渐增加直到基体单元几乎全都发生了破坏；巷道表面附近的界面与基体单元几乎全都发生破坏的范围增加。考虑基体的连续介质模型的结果与不考虑基体的连续介质模型的结果有明显差异，但是也具有一定类似性，它们应力绝对值低的位置大致相同，而不考虑基体的连续介质模型应力的波动幅度更大。

（2）当基体强度参数取为与界面的相等时，计算出的结果与传统的连续介质模型所得的结果类似，剪切应变场呈圆环形分布，且离巷道表面越近，其值越高。随着基体强度参数的降低，计算结果由光滑的剪切应变场变得斑驳，直到形成相互交织的滑移线网，且网的范围逐渐增加。滑移线网的出现应该与由坚硬的弹性颗粒单元和软弱的基体单元构成的材料的非均质性有关。这与不考虑基体的连续介质模型的结果有明显的差异（多个环向的应变集中区）。

［1］MartinCD.Theeffectofcohesionlossandstresspathon brittlerockstrength［J］.CanadianGeotechnicalJournal，1997，34（5）：698－725.

［2］MartinCD，KaiserPK，MccreathDR.Hoek-Brown parametersforpredicting the depth ofbrittle failure aroundtunnels［J］. Canadian GeotechnicalJournal，1999，36：136－151.

［3］祝方才，潘长良，曹平.基于各向异性损伤的厚壁圆筒应力分布有限元分析［J］.岩石力学与工程学报，2003，22（5）：838－842.ZHUFangcai，PANChangliang，CAOPing. Stresses analysisofsurroundingrocksbasedonthickcylinderand anisotropicdamagetheorybyfiniteelementmethod［J］.ChineseJournalofRockMechanicsandEngineering，2003，22（5）：838－842.

［4］陈祥，孙进忠，张杰坤，等.岩爆的判别指标和分级标准及可拓综合判别方法［J］.土木工程学报，2009，42（9）：82－88.CHENXiang，SUNJinzhong，ZHANGJiekun，etal.Judgmentindexesandclassificationcriteriaofrock-burst withtheextensionjudgmentmethod［J］.ChinaCivil EngineeringJournal，2009，42（9）：82－88.

［5］王学滨，潘一山，陶帅.不同尺寸的圆形隧洞剪切应变局部化过程模拟［J］.中国地质灾害与防治学报，2009，20（4）：101－108.WANGXuebin，PANYishan，TAOShuai.Numerical simulationofshearstrainlocalization processesalong circulartunnelswithdifferentdiameters［J］.TheChinese JournalofGeologicalHazardandControl，2009，20（4）：101－108.

［6］王学滨，潘一山.不同侧压系数条件下的圆形巷道岩爆过程模拟［J］.岩土力学，2010，31（6）：1937－1942.WANGXuebin，PANYishan.Numericalsimulationof rockburstprocessesofacirculartunnelatdifferentlateral pressurecoefficients［J］. RockandSoilMechanics，2010，31（6）：1937－1942.

［7］王涛，吕庆，李杨，等.颗粒离散元方法中接触模型的开发［J］.岩石力学与工程学报，2009，28（增2）：4040－4045.WANGTao，LUQing，LIYang，etal.Developmentof contactmodelinparticlediscreteelementmethod［J］.ChineseJournalofRockMechanicsandEngineering，2009，28（S2）：4040－4045.

［8］周健，邓益兵，贾敏才，等.基于颗粒单元接触的二维离散 －连续耦合分析方法［J］.岩土工程学报，2010，32（10）：1479－1484.ZHOUJian，DENGYibing，JIAMincai，etal.Coupling methodoftwo-dimensionaldiscontinuum-continuumbased oncontactbetweenparticleandelement［J］.Chinese JournalofGeotechnicalEngineering，2010，32（10）：1479－1484.

［9］Hazzard J F，Young R P. Momenttensors and micromechanicalmodels［J］. Tectonophysics，2002，356：181－197.

［10］伍小林，王学滨，潘一山.单向压缩条件下圆形颗粒体含孔洞试样的力学模拟［J］.水利水运工程学报，2010（3）：40－44.WUXiaolin，WANGXuebin，PANYishan.Numerical simulation of mechanical behavior of a specimen composedofcirculargranularmaterial［J］. Hydro-ScienceandEngineering，2010（3）：40－44.

［11］伍小林，王学滨，潘一山.含不同半径孔洞的颗粒体模型的力学行为数值模拟［J］.中国地质灾害与防治学报，2011，22（1）：107－114.WUXiaolin，WANGXuebin，PANYishan.Numerical simulationofmechanicalbehaviorofthegranularmodel withholeswithdifferentradii［J］.TheChineseJournalof GeologicalHazardandControl，2011，22（1）：107 －114.

［12］孟陆波，赵建壮.岩爆的洞室效应［J］.中国地质灾害与防治学报，2008，19（2）：59－62.MENGLubo，ZHAOJianzhuang.Thechambereffectsof rockburst［J］.TheChineseJournalofGeologicalHazard andControl，2008，19（2）：59－62.

［13］孟陆波，李天斌，周春宏，等.某水电站洞室岩爆的尺寸效应［J］.水力发电，2008，34（1）：29－31.MENGLubo，LITianbin，ZHOUChunhong，etal.Rock burst'sdimensionaleffectonthetunnelsofahydropower station［J］.WaterPower，2008，34（1）：29－31.

［14］陶帅，王学滨，潘一山，等.线性及非线性屈服函数交界处临界应力对隧洞围岩力学行为的影响［J］.水资源与水工程学报，2011，22（2）：31－36.TAOShuai，WANGXuebin，PANYishan，etal.Effects ofcriticalstresscorrespondingtothetangentialpoint betweenlinearand nonlinearyield functionson the mechanicalbehavioroftunnelsurroundingrock［J］.JournalofWaterResources＆ WaterEngineering，2011，22（2）：31－36.