孔隙煤岩损伤破坏行为的数值模拟

彭瑞东,张玉军,杨永明,刘坚志

(1.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室,北京 100083;2.中国矿业大学(北京)力学与建筑工程学院,北京 100083)

孔隙煤岩损伤破坏行为的数值模拟

彭瑞东1,2,张玉军2,杨永明2,刘坚志2

(1.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室,北京 100083;2.中国矿业大学(北京)力学与建筑工程学院,北京 100083)

利用ANSYS有限元软件,构建了孔隙煤岩细观单元的二维数值计算模型;提出了一种针对孔隙煤岩复杂几何模型的网格优化措施,实现了整体映射网格与局部自由网格相结合的划分方法,提高了生成单元的质量以及求解效率;借助单元生死技术模拟了孔隙煤岩的脆性损伤破坏过程,研究了孔隙率及孔隙形态对煤岩整体力学行为的影响。研究表明,孔隙煤岩的损伤破坏行为和孔隙率、孔隙位置、孔隙形态等几何参数有关。孔隙集中的区域会出现应力或应变集中的现象,随着孔隙率及孔隙长短轴比增大,煤岩细观单元损伤现象更加明显,弹性模量下降,煤岩强度急剧下降,且损伤阈值也有所下降。

煤岩;孔隙;损伤;单元生死

随着我国煤炭开采深度及难度的不断发展,以及煤层气开发的逐渐发展,与此相关的岩石力学问题也变得更加棘手。科学准确地掌握煤岩的非线性力学响应已成为采矿工程中亟需解决的关键基础问题之一。作为一种复杂的地质材料,煤岩具有非均质的不连续的多相复合结构,其内部大量存在不同尺度的不规则的孔隙,这些孔隙往往又是各种伴生气体或水的贮藏空间及运移通道。在外载作用下,煤岩内部原有的各种孔隙将发生变化,局部还可能产生开裂而形成新的孔隙。伴随孔隙结构的变化发展,某些区域的孔隙将逐渐贯通并形成宏观裂缝,进而导致整体的失稳破坏。正因为如此,煤岩的变形破坏过程实质上是一个十分复杂的损伤演化过程,表现为内部孔隙结构的不断演化。研究揭示煤岩损伤破坏行为的规律特点已成为岩石力学研究的重点、难点和热点问题。

数值模拟在工程计算中已经体现出巨大的优势, ANSYS,FLAC等软件的应用对于煤炭开采中的应力场、变形场、渗流场分析计算发挥了重要作用[1-4]。为了更加合理的反映岩石的变形破坏行为,如何在数值模拟中实现岩石损伤软化的模拟成为关注的焦点问题。一些损伤本构模型被提出并引入到数值模拟软件中,从而模拟出了岩石、混凝土等材料的峰后力学行为[5-9]。一些能量判别准则被提出并引入到数值模拟中,从而实现了岩石中损伤及裂纹扩展的模拟[10-13]。唐春安等将岩石材料性质的统计分布特点引入到数值模拟中,采用损伤软化本构模型,并实现岩石细观单元的破裂处理,从而研发了RFPA软件,可以实现岩石破坏过程的全程模拟[14-16]。还有一些基于离散元的颗粒流模型[17-18]、DDA模型[19-20]也被用来模拟岩石的变形破坏行为。不仅如此,数值计算软件的应用也为煤岩变形破坏机理的研究提供了有力支持,如裂隙扩展规律等可以借助数值模拟进行深入分析研究[21-22]。

进一步通过模拟煤岩内部的孔隙结构特征并对其进行数值计算分析,可以帮助揭示煤岩孔隙结构的损伤演化过程及其对煤岩最终破坏的影响规律。事实上这一思路和方法已在一些相对比较规则的多孔材料中得到应用。例如Chen等[23-24]利用ABAQUS梁单元模拟了金属蜂窝材料的胞壁并对其力学响应进行了预测,借助有限元方法研究了刚性夹杂和孔洞对规则六边形蜂窝材料的弹性模量和屈服强度的影响。Chung和Waas[25]用有限元对聚碳酸醋脂蜂窝材料在面内受压缩的情况进行了数值模拟。鞠杨等[26-28]通过砂岩及泡沫混凝土的CT扫描实验研究了岩石孔隙的几何特征与分布规律,建立了岩石孔隙结构的统计模型并结合有限元软件进行了数值模拟研究。岳中琦和陈沙[29-30]综合数字图像处理理论、几何矢量转换技术与有限元网格自动生成原理,建立了岩土工程材料的数字图像有限元分析模型。此外还有一些离散元模型也被采用[31-32]。这些研究有力的推动了对孔隙煤岩变形破坏行为及机理的研究。但目前的有限元模拟研究还主要局限于对含有等轴圆孔的高孔隙率岩石材料的研究[33],这对于砂岩等富含孔隙的岩石是适用的。但还有一些岩石孔隙率较低,且在地质作用过程中因受压而以狭长孔隙形式存在。较之圆形孔隙,这些具有狭长孔隙且孔隙率较低的岩石材料给数值模拟带来一些新的困难,因此有必要进一步对此展开研究。

本文将借助ANSYS有限元软件,探讨如何建立具有随机分布孔隙的煤岩细观单元,如何优化网格划分建立合适的数值计算模型,并通过单元生死技术模拟孔隙煤岩的损伤破坏行为,研究孔隙结构对煤岩损伤破坏行为的影响。本文研究主要解决复杂孔隙煤岩数值模拟方面的一些共性问题,包括几何建模、单元划分、损伤模拟等,研究结果有助于揭示孔隙分布对煤岩整体力学行为的影响。

1 孔隙煤岩二维数值模型的构建

煤岩的孔隙结构直接影响着其力学行为特点。构建能反映煤岩孔隙结构特征的细观代表性体积单元(RVE)并进行网格划分生成有限元模型成为孔隙煤岩数值模拟研究的关键基础问题。

1.1 孔隙煤岩的细观单元

孔隙煤岩是由固相矿物和通过固相形成的孔隙所组成的复合体,它区别于普通密实固体材料的最显著特点是具有孔隙。因此,孔隙煤岩最基本的结构参量应是直接表征其孔隙性状的指标,如孔隙率、孔隙形貌、孔隙尺寸及其分布等。大量CT及SEM实验表明,孔隙分布多为均匀分布,而孔径一般多服从指数分布和正态分布[26,33]。而且除了等轴圆孔,岩石中还赋存有大量狭长孔隙。尤其是对于具有较低孔隙率的岩浆岩、变质岩及部分致密的沉积岩,由于经受长期地质作用,其中的孔隙往往以狭长孔隙形式存在。在工程尺度上的裂隙岩体其中的裂隙也是狭长型的。根据以前的实验结果,本文构建了二维煤岩细观单元,单元边长为100μm,其中含有随机分布的圆形或椭圆孔隙,孔隙位置采用均匀分布,孔径服从指数分布,并考虑了不同的孔隙率、不同长短轴比的影响。

通过FORTRAN或Matlab等进行简单的程序设计,即可按特定的概率统计规律生成所需的孔隙大小、倾角及位置坐标,进而利用APDL程序可以在ANSYS中将其作为圆心坐标及长短轴建立一系列狭长椭圆孔,再通过布尔运算将这些狭长椭圆孔从正方形区域中减去即可得到复杂多材料的二维模型。为考虑孔隙率的影响,分别建立了5种孔隙煤岩细观单元,其中孔隙率分别为0.130%,1.102%,5.132%, 9.543%及19.580%。针对孔隙率为0.130%的低孔隙率情况,还考虑了长短轴比分别为1∶1,2∶1,5∶1,10∶1的不同孔隙形态。图1为孔隙率0.130%、孔隙长短轴比为10∶1的煤岩细观单元。

图1 孔隙率0.13%的二维煤岩细观单元Fig.1 2Dmesoscopic rock unitwith 0.13%porosity

1.2 细观单元的网格优化

在建立了孔隙煤岩细观单元的几何模型之后,需要选择合适的方法将几何模型分割为一系列有限元单元的组合体,即进行网格划分。对于有限元分析来说,网格划分的好坏直接影响到求解的精度和速度。在ANSYS中提供了自由网格划分和映射网格划分2种方法。对于复杂多孔材料模型,由于各种孔隙的随机分布,使得几何模型十分复杂,因此难以进行映射网格划分,直接划分网格只能采用自由网格划分。孔隙煤岩二维细观单元的自由网格划分情况如图2所示。可以明显看出,由于结构中狭长孔隙的存在,导致孔隙之间的网格划分极不规则,而且随着孔隙率的增大,网格划分将越来越不规则。这样自由划分的网格,不仅单元的数目多,而且单元质量差,尤其是在几何形状突变处易形成过于细密和畸形的网格,这不仅会增加计算时间,还会造成模型刚度和质量矩阵奇异,对整个模型的求解非常不利,尤其是在非线性计算中甚至会导致计算不收敛。因此有必要寻找一种改进的网格划分方法。

考虑到复杂多孔材料的应力应变主要在孔隙附近发生较大变化,笔者只希望在孔隙结构附近区域划分比较细的网格,而在孔隙结构以外的区域能有相对粗疏一些的网格来表示,这样就能降低计算规模。因此可建立比孔隙直径稍大的一个矩形区域将孔隙包裹起来,该区域内孔隙的边上划分较密的网格,区域外围的边与煤岩基体相连接,划分为较疏的网格。该区域也正好反映了孔隙附近的应力应变集中区域的大小。将各个孔隙周围的这些区域去除后,基体就成为一个规则的区域,可以进行映射网格划分,生成质量较好的四边形单元。而这些区域内部采用自由网格划分,生成由疏到密逐渐过渡的三角形单元。

图2 孔隙煤岩细观单元有限元网格的自由划分Fig.2 Freemeshing for finite elements of porous rock unit

图3 孔隙煤岩二维模型网格优化方法Fig.3 Freemeshing for finite elements of porous rock unit

以前述二维孔隙煤岩细观单元模型为例,模型的基体是一个矩形结构,该基体由2μm×2μm的小正方形用AGLUE命令黏接而成(图3(a)),在有孔隙存在的位置建立边长比孔隙直径大的正方形,用布尔运算中ASBA命令在基体上减去这些正方形(图3(b)),然后再逐个创建含椭圆孔隙的这些小矩形,用布尔运算中的AGLUE命令将所有面粘贴在一起成为一个实体(图3(c))。这样,小矩形外的部分可以用映射方法划分网格,得到比较规整的网格,同时可通过定义孔隙每个边的份数实现对网格的控制,优化应力集中区域(孔隙周围)网格划分(图3(d))。由于采用了AGLUE命令,各矩形区域内外共用一条边,通过边上的节点确保两侧网格相连,物理场可以有效传递。

图3(d)为采用该方法生成的孔隙率为0.130%的煤岩细观单元模型的网格划分结果,对比图2可以看到,与前面将模型作为一个整体自动划分的的网格情况相比,改进后的方法划分的网格规整的多,而且单元质量比较好。这样不仅可以减少计算时间,也大大提高了最后计算结果的精度。通过加载计算对2种模型进行比较,前者计算时间较长,而且不能在塑性范围内进行求解,而后者则在弹性和弹塑性范围内都能够求解,而且求解速度较快,求解精度也比较高。这就说明这种改进的方案是可行的。

1.3 单孔模型ANSYS数值模拟结果的验证

在有限元数值计算中,网格划分疏密程度对结果有着重要影响。为了确定验证孔隙煤岩细观单元模型,并确定合适的网格划分密度,对只含一个椭圆孔隙的单孔模型进行了理论分析和数值计算。细观单元边长为100μm,椭圆孔位于单元中心,短轴为0.4μm,长轴为4μm。有限元网格采用平面183单元。模型左侧简支约束,右侧施加均匀拉应力q= 1 MPa。计算得到的第1主应力云图如图4所示。

图4 单孔模型第1主应力云图Fig.4 First principal stress contour of singe poremodel

据文献[34]可知,对于距离边界较远处的椭圆形孔口,其长轴和短轴分别为2a和2b,孔边不受面力,当在与长轴成α角的方向受到均匀拉应力q时,孔边应力可按下式计算:

式中,α为x轴与作用力的交角,x轴和椭圆长轴重合;σθ为孔边应力;θ为极坐标系中的极角;m为描述孔形状的参数,取m=(a-b)/(a+b)。当α=0时,即水平椭圆孔的情况,最大主应力位于θ=±π/2处,大小为(1+2b/a)q。当α=π/2时,即水垂直椭圆孔的情况,最大主应力位于θ=±π/2处,大小为(1+2a/ b)q。

表1列出了长短轴之比为10的狭长椭圆孔隙端部的最大主应力的理论计算结果和数值计算结果。对比可以看出,最大应力值、最大应力位置的数值解与解析解吻合得相当好。

表1 ANSYS数值模拟解与理论解析解的对比Table 1 Theoretical solution and ANSYS numerical solution

2 孔隙煤岩的损伤破坏行为模拟

2.1 基于单元生死的孔隙煤岩损伤破坏模拟方法

对于比较致密的孔隙煤岩这类脆性材料,在出现应力集中后一般不会发生塑性变形而产生应变集中现象,也没有塑性强化现象,而是会在高应力区发生开裂,表现出损伤软化现象,从而使煤岩弹模、强度等下降。为了模拟材料的损伤,在ANSYS有限元分析中需借助单元生死技术。

每级加载结束后,需要根据有限元单元的应力应变等参数判断是否杀死该单元,可以采用最大拉应力或Mises应力进行判别,这些应力值可在ANSYS中直接得到。但孔隙煤岩的失效破坏通常采用莫尔强度理论来判别,为此需要借助单元表来计算等效应力。据文献[35]可知,基于莫尔理论,煤岩的失效应力为

式中,rM13第3主应力;[σt]为煤岩在单轴拉伸时的许用拉应力;[σc]为煤岩在单轴压缩时的许用拉应力。

于是可应用＜ETABLE＞命令建立第1主应力表以及第3主应力表,再通过＜SADD＞命令按照式(2)建立莫尔等效应力表。当某个有限元单元的第1主应力或莫尔等效应力大于临界值时,就认为这个有限元单元失效,将其杀死。基于单元表进行选择操作,就可将全部失效单元选出并杀死。然后再进行下一级的加载,直至整个细观单元失效破坏。

下面分别针对前述的几种孔隙率煤岩细观单元进行模拟分析。有限元网格采用平面183单元。现只考虑脆性损伤,煤岩基体采用线弹性本构模型,取弹性模量E为20 GPa,泊松比ν为0.3,最大许用拉应力及莫尔等效应力的临界值为30 MPa。模型左侧简支约束,右侧逐级施加位移载荷d模拟单轴压缩情形。

2.2 荷载作用下孔隙煤岩中的应力分布特征

图5 不同孔隙率时煤岩细观单元的损伤破坏过程Fig.5 Damage and failure process of rock units at different porosity

图5所示为不同孔隙率的煤岩细观单元在采用单元生死技术后的损伤破坏过程。图中标注出了当前施加在右侧的位移大小d,由于不同孔隙率下峰值载荷不同,加载步长d也不相同。图示为莫尔应力分布云图,煤岩基体的应力绝大部分也是一致的,只是孔边附近出现应力集中,但由于损伤开裂,应力值并不会超过30 MPa。且由于孔隙之间的相互作用,微裂隙的演化形式比较复杂。当孔隙率较低时,微裂纹只是在孔隙的两端发展,最终汇聚成一条主裂纹导致整体丧失承载能力。当孔隙率较高时,可能会出现微裂纹与孔隙的汇聚,产生复杂的裂隙网络,最终形成多条裂纹并贯通,从而引起整体丧失承载能力。

图6为不同孔隙形态下煤岩细观单元在采用单元生死技术后的损伤破坏过程。图示为莫尔应力分布云图,如前所述,煤岩基体的应力绝大部分也是一致的,只是孔边附近出现应力集中。从图6可看出,随着孔隙长短轴比增大,孔隙形态由等轴圆形向狭长缝形变化,应力集中程度更加明显,应力梯度更加剧烈,这也意味着狭长孔隙导致更为严重的应力集中。

图6 不同孔隙长短轴比时煤岩细观单元的损伤破坏过程Fig.6 Damage and failure process of rock units at different pore aspect ratio

3 孔隙煤岩的损伤破坏规律分析

在不同的孔隙率下,煤岩细观单元的力学响应是不同的。通过数值计算,不仅可以模拟得到煤岩内部的应力应变分布云图,还可进一步计算得到煤岩整体的刚度、强度指标,并通过统计杀死单元数目定量计算煤岩的损伤情况,从而可以定量描述不同孔隙特征下煤岩的整体力学响应及损伤破坏规律。

3.1 孔隙煤岩的应力-应变曲线

当施加位移载荷为Δl时,煤岩细观单元的名义应变为

其中,l0为细观单元的边长100μm。相应的名义应力为

其中,Ri为煤岩细观单元左侧边上各结点的约束反力,即水平坐标x=0的各结点的约束反力。因所分析的是二维模型,δ为细观单元的厚度,即在ANSYS中的二维平面单元183的厚度,一般取为单位厚度1。图7(a)为各种孔隙率的煤岩细观单元在不同载荷步下计算得到的名义应力-名义应变曲线,图7(b)为不同孔隙形态下的煤岩细观单元名义应力-名义应变曲线。

图7 不同孔隙率和孔隙长短轴比时煤岩细观单元的名义应力-名义应变曲线Fig.7 Nominal stress-strain curves of rock units at different porosity and pore aspect ratio

当煤岩内含有孔隙时,内部的应力场应变场是不均一的,在整体名义应力不大的情况下,局部也会出现较高应力或应变,并在这些局部区域发生开裂,如前述图5,6所示。随着微裂纹的扩展,煤岩损伤逐渐加剧,加载刚度不断下降,直至达到某一临界值后,煤岩的变形已经完全可以通过内部弹性能的释放引发微裂纹扩展来实现,在名义应力-名义应变曲线上表现为峰值载荷后的软化下降段,且随着裂纹汇聚贯通,损伤愈加严重,卸载刚度急剧增大,如图7所示。

3.2 孔隙煤岩的弹性特征

煤岩细观单元的名义弹性模量可根据名义应力-名义应变曲线的线性段斜率求的,即

表2列出了各种孔隙率的煤岩细观单元的弹性模量。由于在线弹性阶段没有损伤,无论是否采用单元生死技术都得到相同的弹模计算结果。

表2 不同孔隙率煤岩细观单元的力学性能指标Table 2 Mechanical p roperties of rock unitswith d ifferent porosity

从表2可以看出,不同孔隙率的煤岩细观单元的名义弹性模量是不同的,孔隙率的变化对煤岩的名义弹性模量有较大影响。随着孔隙率的增大,煤岩的名义弹性模量减小,如图8(a)所示。

图8(b)为相同孔隙率下不同孔隙形态煤岩细观单元的名义弹性模量,可见孔隙长短轴比对煤岩细观单元弹性模量的影响不大。尤其当孔隙率较低时,这种影响几乎可以忽略。当孔隙率较大时,在相同孔隙率下,相比具有圆形孔隙的煤岩,具有狭长孔隙的煤岩弹模略低一点。结合图7(b)可以看出,这种变化与具有狭长孔隙的煤岩更易发生损伤有关,其应力应变曲线上的刚度下降段出现得更早。

3.3 孔隙煤岩的强度特征

从图7(a)和表2可以看出,不同孔隙率下煤岩细观单元的极限强度是不同的,不含孔隙(即孔隙率为0)的煤岩强度最大,随着孔隙率增大煤岩强度急剧下降。可见,这含量甚微的孔隙对煤岩的力学特性有着重要影响。

图8 不同孔隙率和孔隙长短轴比时煤岩细观单元的名义弹性模量Fig.8 Nom inal elastic modulus of rock units at different porosity and pore aspect ratio

当煤岩不含孔隙时,内部应力场均一,仅当整个模型的第1主应力或莫尔等效应力达到最大许用应力时,煤岩整体破坏,应力降为0,表现出明显的脆断特征。当煤岩内含有孔隙时,内部的应力场不再均一,在整体名义应力不大的情况下,局部也会出现较高应力,并在这些局部区域发生开裂,形成损伤。借助单元生死技术,可以模拟出煤岩损伤后的应力重分布情况,从而得到具有峰值点及软化下降段的名义应力-应变曲线,如前述图7所示。从图中可以看出,在不同孔隙率下,这种力学响应过程是类似的。不过随着孔隙率的增大,损伤程度也将更为严重,因此极限强度和加载刚度的下降、卸载刚度的增加也更为明显。

表2列出了各种孔隙率的煤岩细观单元的极限强度,其对应应变值为极限应变。不采用单元生死技术时,煤岩细观单元的名义应力将随加载的进行不断增大。表2也列出了不采用单元生死技术时煤岩细观单元在达到极限应变时的名义应力值,不难发现该值要大于考虑损伤时的强度值。由此可知,损伤会导致煤岩强度的下降。

不同孔隙率下煤岩细观单元的强度也在图9(a)中给出,可见随着孔隙率的增大,煤岩强度急剧下降。图9(b)为相同孔隙率下不同孔隙形态煤岩细观单元的强度,可知,具有狭长孔隙的煤岩要比具有圆形孔隙的煤岩强度低,因为具有狭长孔隙的煤岩更易发生损伤。

3.4 孔隙煤岩的损伤演化规律

为了定量研究孔隙煤岩损伤破坏过程中的演化规律,基于损伤力学理论定义损伤变量为

其中,A0为细观单元的原始面积;l0为细观单元的边长100μm;φ为孔隙率;AD为细观单元中损伤区域的面积;ADi为损伤失效的有限元单元的面积。图10(a)为各种孔隙率的煤岩细观单元在不同载荷步下计算得到的损伤变量演化曲线,图10(b)为不同孔隙形态下的煤岩细观单元损伤变量演化曲线。

图9 不同孔隙率和孔隙长短轴比时煤岩细观单元的强度Fig.9 Strength of rock units at different porosity and pore aspect ratio

由于孔隙煤岩初始时便存在孔隙,因此计算出的初始损伤变量并不为0,可称之为初始损伤,这反映了煤岩孔隙率的高低。孔隙率越大,初始损伤变量也就越大,如图10(b)所示。

由图10可以看出,对于含有孔隙的煤岩,随着应变的增加,孔隙煤岩的损伤面积在不断增大,而且损伤面积的增长速率越来越大。这就表明随着变形的发展,孔隙煤岩的损伤程度由慢到快不断加剧。

根据图7(a)所示不同孔隙率下煤岩的极限强度,可以对应确定出相应的损伤变量大小,称之为损伤阈值,计算结果列在表2中。损伤阈值是煤岩发生自主裂纹扩展的临界值。一旦煤岩的损伤变量值大于损伤阈值,煤岩就将发生失稳破坏。从图10(a)和表2可以看出,损伤阈值也对应着损伤演化曲线的一个拐点。且随着孔隙率的增大,损伤阈值下降,这表明损伤更易引发煤岩的失稳破坏。

从图10(b)可以看出,随着孔隙长短轴比的增大,损伤阈值下降,这也就表明狭长孔隙较之圆形孔

隙更易导致损伤,因此如前所述,狭长孔隙更易导致煤岩弹性模量、强度降低。在工程实际中,各种节理可模拟为狭长孔隙,他们对煤岩整体力学行为的影响更加剧烈,因此需加倍关注。

4 结 论

(1)在对孔隙煤岩模型进行网格划分时,可建立比孔隙直径稍大的一个矩形或正六面体区域将孔隙包裹起来,该区域内孔隙的边上划分较密的网格,区域外围的边与材料基体相连接,划分为较疏的网格,该区域也正好反映了孔隙附近的应力应变集中区域的大小,藉此就可对区域外的基体进行映射网格划分,生成质量较好的四边形或六面体单元,而对区域内部采用自由网格划分,生成由疏到密逐渐过渡的三角形或四面体单元,这种网格优化措施大大提高了网格划分的质量。

(2)孔隙煤岩的损伤破坏行为与孔隙率、孔隙位置等几何参数有关,对于煤岩类脆性材料,孔隙集中的区域会出现应力集中的现象。

(3)煤岩细观单元的力学响应反映了煤岩中孔隙的影响,研究表明,随着孔隙率增大,损伤现象更加明显,煤岩强度急剧下降,且损伤阈值也有所下降。

运用数值模拟的方法来研究孔隙煤岩的损伤破坏行为是一种比较新颖的方法,这可以充分发挥数值模拟的优势,代替部分试验,本文对此进行了初步的尝试和探讨,但进一步的研究还有待深入,需要通过大量模型的计算,统计分析孔隙率、孔隙分布特征、基体材料属性等对孔隙煤岩损伤破坏行为的影响,从而为揭示孔隙煤岩损伤破坏的机理提供依据。

[1] 武 雄,汪小刚,段庆伟,等.导水断裂带发育高度的数值模拟[J].煤炭学报,2008,33(6):609-612.

Wu Xiong,Wang Xiaogang,Duan Qingwei,etal.Numericalmodeling about developing high ofwater flowing fractured zone[J].Journal of China Coal Society,2008,33(6):609-612.

[2] 赵 鹏,谢凌志,熊 伦.无煤柱开采条件下煤岩体支承压力的数值模拟[J].煤炭学报,2011,36(12):2029-2034.

Zhao Peng,Xie Lingzhi,Xiong Lun.Numerical simulation of abutment pressure in coal for non-pillar mining[J].Journal of China Coal Society,2011,36(12):2029-2034.

[3] 张鹏海,杨天鸿,郑 超,等.基于采动应力场与微震活动性的岩体稳定性分析[J].煤炭学报,2013,38(2):183-188.

Zhang Penghai,Yang Tianhong,Zheng Chao,et al.Analysis of surrounding rock stability based on mining stress field and microseismicity[J].Journal of China Coal Society,2013,38(2):183-188.

[4] 李建功,康建宁,刘 红,等.地震波在巷道弹塑性围岩中传播规律的数值模拟研究[J].煤炭学报,2011,36(S2):282-286.

Li Jiangong,Kang Jianning,Liu Hong,et al.Numerical simulation on propagation laws of seismic wave in elastoplastic tunnel rockmedium [J].Journal of China Coal Society,2011,36(S2):282-286.

[5] 曹 鹏,冯德成,沈新普,等.基于ABAQUS平台的塑性损伤子程序开发及其稳定性研究[J].工程力学,2012,29(S2):101-106.

Cao Peng,Feng Decheng,Shen Xinpu,etal.Developmentofplasticity-damagemodel based on abaqus and algorithm stability analysis [J].Engineering Mechanics,2012,29(S2):101-106.

[6] 龙渝川,李正良.基于能量耗散机制的混凝土受压损伤模型[J].工程力学,2010,27(S2):171-177.

Long Yuchuan,Li Zhengliang.A compressive damagemodel of concrete based on energy dissipation mechanism[J].Engineering Mechanics,2010,27(S2):171-177.

[7] 刘智光,陈健云,白卫峰.基于随机损伤模型的混凝土轴拉破坏过程研究[J].岩石力学与工程学报,2009,28(10):2048-2058.Liu Zhiguang,Chen Jianyun,BaiWeifeng.Research on concrete failure process under uniaxial tension based on stochastic damagemodel [J].Chinese Journal of Rock Mechanics and Engineering,2009, 28(10):2048-2058.

[8] 赵吉坤,张子明,刘仲秋,等.大理岩破坏过程的三维细观弹塑性损伤模拟研究[J].岩土工程学报,2008,30(9):1309-1315.

Zhao Jikun,Zhang Ziming,Liu Zhongqiu,et al.3D numerical simulation of elasto-plastic damage and failure process of marble rock [J].Chinese Journal of Geotechnical Engineering,2008,30(9): 1309-1315.

[9] 杨伯源,李运军,汪忠明.基于非局部模型的混凝土损伤的数值模拟[J].合肥工业大学学报(自然科学版),2007,30(4):472-474.

Yang Boyuan,Li Yunjun,Wang Zhongming.Numerical simulation of concretematerial based on the non-local damagemodel[J].Journal of Hefei University of Technology,2007,30(4):472-474.

[10] 李树忱,冯现大,李术才,等.深部岩体分区破裂化现象数值模拟[J].岩石力学与工程学报,2011,30(7):1337-1344.

Li Shuchen,Feng Xianda,Li Shucai,et al.Numerical simulation of zonal disintegration for deep rock mass[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(7):1337-1344.

[11] 孙 倩,李树忱,冯现大,等.基于应变能密度理论的岩石破裂数值模拟方法研究[J].岩土力学,2011,32(5):1575-1582.

Sun Qian,LiShuchen,Feng Xianda,et al.Study ofnumerical simulation method of rock fracture based on strain energy density theory [J].Rock and Soil Mechanics,2011,32(5):1575-1582.

[12] 蓝 航,潘俊锋,彭永伟.煤岩动力灾害能量机理的数值模拟[J].煤炭学报,2010,35(S1):10-14.

Lan Hang,Pan Junfeng,Peng Yongwei.Numerical simulation for energymechanism of underground dynamic disaster[J].Journal of China Coal Society,2010,35(S1):10-14.

[13] 王来贵,赵 娜,何 峰,等.岩石蠕变损伤模型及其稳定性分析[J].煤炭学报,2009,34(1):64-68.

Wang Laigui,Zhao Na,He Feng,et al.Rock creep damage model and its stability analysis[J].Journal of China Coal Society,2009, 34(1):64-68.

[14] 唐春安,芮勇勤,刘红元,等.含瓦斯“试样”突出现象的RFPA2D数值模拟[J].煤炭学报,2000,25(5):501-505.

Tang Chun’an,Rui Yongqin,Liu Hongyuan,et al.Numerical simulation to outburst mechanism of coal or rock containing gas with FRPA2Dsystem[J].Journal of China Coal Society,2000,25(5): 501-505.

[15] 梁正召,唐春安,张永彬,等.岩石三维破裂过程的数值模拟研究[J].岩石力学与工程学报,2006,25(5):931-936.

Liang Zhengzhao,Tang Chunan,Zhang Yongbin,et al.3D numerical simulation of failure process of rock[J].Chinese Journal of Rock Mechanics and Engineering,2006,25(5):931-936.

[16] 朱万成,唐春安,杨天鸿,等.岩石破裂过程分析用(RFPA2D)系统的细观单元本构关系及验证[J].岩石力学与工程学报, 2003,22(1):24-29.

Zhu Wanchcng,Tang Chun’an,Yang Tianhong,et al.Constitutive relationship ofmesoscopic elements used in RFPA2Dand its validations[J].Chinese Journal of Rock Mechanics and Engineering, 2003,22(1):24-29.

[17] 黄 达,岑夺丰.单轴静-动相继压缩下单裂隙岩样力学响应及能量耗散机制颗粒流模拟[J].岩石力学与工程学报,2013, 32(9):1926-1936.

Huang Da,Cen Duofeng.Mechanical responses and energy dissipationmechanism of rock specimen with a single fissure under static and dynamic uniaxial compression using particle flow code simulations[J].Chinese Journal of Rock Mechanics and Engineering, 2013,32(9):1926-1936.

[18] 马 刚,周 伟,常晓林,等.堆石体三轴剪切试验的三维细观数值模拟[J].岩土工程学报,2011,33(5):746-753.

Ma Gang,Zhou Wei,Chang Xiaolin,et al.3Dmesoscopic numerical simulation of triaxial shear tests for rockfill[J].Chinese Journal of Geotechnical Engineering,2011,33(5):746-753.

[19] 高亚楠,程红梅,李玺茹,等.单一煤层钻孔割缝卸压增透的非连续变形分析(DDA)模拟[J].煤炭学报,2011,36(12):2068-2073.

Gao Yanan,Cheng Hongmei,Li Xiru,et al.The numerical modelling of slotted drilling in single coal seam based on DDA[J].Journal of China Coal Society,2011,36(12):2068-2073.

[20] 左建平,陈忠辉,王怀文,等.深部煤矿采动诱发断层活动规律[J].煤炭学报,2009,34(3):305-309.

Zuo Jianping,Chen Zhonghui,Wang Huaiwen,et al.Experimental investigation on fault activation pattern under deep mining[J].Journal of China Coal Society,2009,34(3):305-309.

[21] 郭德勇,吕鹏飞,裴海波,等.煤层深孔聚能爆破裂隙扩展数值模拟[J].煤炭学报,2012(2):274-278.

Guo Deyong,LüPengfei,Pei Haibo,et al.Numerical simulation on crack propagation of coal bed deep-hole cumulative blasting[J].Journal of China Coal Society,2012(2):274-278.

[22] 付金伟,朱维申,曹冠华,等.岩石中三维单裂隙扩展过程的试验研究和数值模拟[J].煤炭学报,2013,37(3):411-417.

Fu Jinwei,Zhu Weishen,Cao Guanhua,et al.Experimental study and numerical simulation of propagation and coalescence process of a single three-dimensional flaw in rocks[J].Journal of China Coal Society,2013,37(3):411-417.

[23] Chen C,Lu T J,Fleck N A.Effect of imperfections on the yielding of two-dimensional foams[J].Journalof the Mechanics and Physics of Solids,1999,47(11):2235-2272.

[24] Chen C,Lu T J,Fleck N A.Effect of inclusions and holes on the stiffness and strength of honeycombs[J].International Journal of Mechanical Sciences,2001,43(2):487-504.

[25] Chung J,Waas A M.Compressive response of circular cell polycarbonate honeycombs under inp lane biaxial static and dynamic loading-Part II:simulations[J].International Journal of Impact Engineering,2002,27(10):1015-1047.

[26] 鞠 杨,杨永明,宋振铎,等.岩石孔隙结构的统计模型[J].中国科学E辑:技术科学,2008,38(7):1026-1041.

Ju Yang,Yang Yongming,Song Zhenduo,et al.A statistical model for porous structure of rocks[J].Science in China Series E:Technological Sciences,2008,38(7):1026-1041.

[27] 鞠 杨,杨永明,毛彦喆,等.孔隙介质中应力波传播机制的实验研究[J].中国科学E辑:技术科学,2009,39(5):904-918.

Ju Yang,Yang Yongming,Mao Yanzhe,et al.Laboratory investigation on mechanisms of stress wave propagations in porous media [J].Science in China Series E:Technological Sciences,2009, 39(5):904-918.

[28] 鞠 杨,王会杰,杨永明,等.应力波作用下岩石类孔隙介质变形破坏与能量耗散机制的数值模拟研究[J].中国科学E辑:技术科学,2010,40(6):711-726.

Ju Yang,Wang Huijie,Yang Yongming,et al.Numerical simulation ofmechanisms of deformation,failure and energy dissipation in porous rockmedia subjected towave stresses[J].Science in China Series E:Technological Sciences,2010,40(6):711-726.

[29] 岳中琦,陈 沙,郑 宏,等.岩土工程材料的数字图像有限元分析[J].岩石力学与工程学报,2004,23(6):889-897.

Yue Zhongqi,Chen Sha,Zheng Hong,et al.Digital image proceeding based on finite element method for geomaterials[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(6):889-897.

[30] 陈 沙,岳中琦,谭国焕.基于数字图像的非均质岩土工程材料的数值分析方法[J].岩土工程学报,2005,27(8):956-964.

Chen Sha,Yue Zhongqi,Tan Guohuan.Digital image based numericalmodeling method for heterogeneous geomaterials[J].Chinese Journal of Geotechnical Engineering,2005,27(8):956-964.

[31] 王宝善,陈 颙,葛洪魁,等.高孔隙岩石变形的离散单元模型[J].地球物理学报,2005,48(6):1336-1342.

Wang Baoshan,Chen Yu,Ge Hongkui,et al.A discrete elementmodel of porous rock deformation[J].Chinese J.Geophys., 2005,48(6):1336-1342.

[32] Baxevanis T,Dufour F,Gilles P C.Interface crack propagation in porous and time-dependentmaterials analyzed with discretemodels [J].Int.J.Fractal,2006(141):561-571.

[33] 彭瑞东,凌天龙,杨永明,等.基于ANSYS的复杂多孔材料数值模拟方法研究[J].中国科技论文在线精品论文,2010,3(15): 1584-1593.

Peng Ruidong,Ling Tianlong,Yang Yongming,et al.Research on the numerical simulation of complex porousmaterials based on ANSYS[J].Highlights of Sciencepaper Online,2010,3(15):1584-1593.

[34] 徐芝纶.弹性力学[M].北京:高等教育出版社,1978.

Xu Zhilun.Elastic mechanic[M].Beijing:Higher Education Press,1978.

[35] 孙训方,方孝淑,关来泰.材料力学[M].北京:高等教育出版社,1979.

Sun Xunfang,Fang Xiaoshu,Guan Laitai.Mechanics of materials [M].Beijing:Higher Education Press,1979.

Num erical sim ulation of porous rock damage and failure

PENG Rui-dong1,2,ZHANG Yu-jun2,YANG Yong-ming2,LIU Jian-zhi2

(1.State Key Laboratory ofCoal Resources and Safe Mining,China University ofMining and Technology(Beijing),Beijing 100083,China;2.School ofMechanics and Civil Engineering,China University ofMining and Technology(Beijing),Beijing 100083,China)

A serial of 2D numericalmodels representing mesoscopic units of porous rocks were constructed and then damage and failure behavior of rockswas simulated.The effect of porosity and pore shape on the overallmechanical behavior of porous rockswas discussed.An improved meshing approach for geometricmodels of porous rockswere proposed to put globalmappingmeshing and local free meshing into together,and consequently element quality and solving efficiency were improved.Element birth and death were introduced to simulate the brittle damage and failure processes of porous rocks.Itwas indicated that damage and failure of porous rocks depends on the geometric parameters such as porosity,pore position and pore shape.Stress or strain concentration appears in those area surrounded with more pores.With the increase of porosity of porous rocks and aspect ratio of pores,the damage of rocks becomesmore severe,the elastic modulus of rocks decreases,the strength of rocks sharply drops and the critical damage variable also declines.

rocks;porous;damage;element birth and death

TD315

A

0253-9993(2014)06-1039-010

彭瑞东,张玉军,杨永明,等.孔隙煤岩损伤破坏行为的数值模拟[J].煤炭学报,2014,39(6):1039-1048.

10.13225/j.cnki.jccs.2013.1536

Peng Ruidong,Zhang Yujun,Yang Yongming,et al.Numerical simulation of porous rock damage and failure[J].Journal of China Coal Society,2014,39(6):1039-1048.doi:10.13225/j.cnki.jccs.2013.1536

2013-10-20 责任编辑:常 琛

国家自然科学基金青年科学基金资助项目(10802092);教育部新世纪优秀人才支持计划资助项目(NCET-12-0966);中央高校基本科研业务费专项基金资助项目(2009QM03)

彭瑞东(1974—),男,山西忻州人,副教授。Tel:010-62331253,E-mail:prd@cum tb.edu.cn