巴西圆盘劈裂二维及三维数值模拟研究

工程科学学报 Chinese Journal of Engineering 巴西圆盘劈裂二维及三维数值棋拟研究 乔兰刘建李庆文赵国彦 Numerical study of the Brazilian tensile test:2D and 3D simulations QIAO Lan,LIU Jian,LI Qing-wen,ZHAO Guo-yan 引用本文： 乔兰，刘建，李庆文，赵国彦.巴西圆盘劈裂二维及三维数值模拟研究J.工程科学学报，2022,44(1)：131-142.doi: 10.13374j.issn2095-9389.2020.09.28.006 QIAO Lan,LIU Jian,LI Qing-wen,ZHAO Guo-yan.Numerical study of the Brazilian tensile test:2D and 3D simulations[J].Chinese Journal of Engineering,.2022,441131-142.doi:10.13374.issn2095-9389.2020.09.28.006 在线阅读View online::htps/ldoi.org10.13374/.issn2095-9389.2020.09.28.006 您可能感兴趣的其他文章 Articles you may be interested in 不同荷载条件下低孔隙率砂岩巴西劈裂试验的声发射特性 Acoustic emission characteristics of Brazilian test for low-porosity sandstone under different load conditions 工程科学学报.2020.42(8：988 https:/1doi.org10.13374.issn2095-9389.2019.08.12.004 考虑岩石交界面方向效应的巴西劈裂试验研究 Experimental study of directivity effect of rock interface under Brazilian splitting 工程科学学报.2017,399%：1295 https:/doi.org10.13374.issn2095-9389.2017.09.001 三维软硬互层边坡的破坏模式与稳定性研究 Numerical analysis of the failure modes and stability of 3D slopes with interbreeding of soft and hard rocks 工程科学学报.2017,392182 https::/1doi.org/10.13374j.issn2095-9389.2017.02.003 盐岩巴西劈裂损伤愈合特性实验研究 Experimental study of the self-healing property of damaged salt rock by Brazilian splitting 工程科学学报.2020,42(5：570 https:/1doi.org10.13374.issn2095-9389.2019.06.04.001 劈裂荷载下的岩石声发射及微观破裂特性 Acoustic emission and micro-rupture characteristics of rocks under Brazilian splitting load 工程科学学报.2019.41(11：1422 https::/doi.org10.13374.issn2095-9389.2018.11.29.005 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 Three-dimensional microscopic model reconstruction of basalt and numerical direct tension tests 工程科学学报.2019,41(8：997 https:1doi.org10.13374.issn2095-9389.2019.08.005

巴西圆盘劈裂二维及三维数值模拟研究 乔兰 刘建 李庆文 赵国彦 Numerical study of the Brazilian tensile test: 2D and 3D simulations QIAO Lan, LIU Jian, LI Qing-wen, ZHAO Guo-yan 引用本文: 乔兰, 刘建, 李庆文, 赵国彦. 巴西圆盘劈裂二维及三维数值模拟研究[J]. 工程科学学报, 2022, 44(1): 131-142. doi: 10.13374/j.issn2095-9389.2020.09.28.006 QIAO Lan, LIU Jian, LI Qing-wen, ZHAO Guo-yan. Numerical study of the Brazilian tensile test: 2D and 3D simulations[J]. Chinese Journal of Engineering, 2022, 44(1): 131-142. doi: 10.13374/j.issn2095-9389.2020.09.28.006 在线阅读 View online: https://doi.org/10.13374/j.issn2095-9389.2020.09.28.006 您可能感兴趣的其他文章 Articles you may be interested in 不同荷载条件下低孔隙率砂岩巴西劈裂试验的声发射特性 Acoustic emission characteristics of Brazilian test for low-porosity sandstone under different load conditions 工程科学学报. 2020, 42(8): 988 https://doi.org/10.13374/j.issn2095-9389.2019.08.12.004 考虑岩石交界面方向效应的巴西劈裂试验研究 Experimental study of directivity effect of rock interface under Brazilian splitting 工程科学学报. 2017, 39(9): 1295 https://doi.org/10.13374/j.issn2095-9389.2017.09.001 三维软硬互层边坡的破坏模式与稳定性研究 Numerical analysis of the failure modes and stability of 3D slopes with interbreeding of soft and hard rocks 工程科学学报. 2017, 39(2): 182 https://doi.org/10.13374/j.issn2095-9389.2017.02.003 盐岩巴西劈裂损伤愈合特性实验研究 Experimental study of the self-healing property of damaged salt rock by Brazilian splitting 工程科学学报. 2020, 42(5): 570 https://doi.org/10.13374/j.issn2095-9389.2019.06.04.001 劈裂荷载下的岩石声发射及微观破裂特性 Acoustic emission and micro-rupture characteristics of rocks under Brazilian splitting load 工程科学学报. 2019, 41(11): 1422 https://doi.org/10.13374/j.issn2095-9389.2018.11.29.005 玄武岩三维细观孔隙模型重构与直接拉伸数值试验 Three-dimensional microscopic model reconstruction of basalt and numerical direct tension tests 工程科学学报. 2019, 41(8): 997 https://doi.org/10.13374/j.issn2095-9389.2019.08.005

工程科学学报.第44卷，第1期：131-142.2022年1月 Chinese Journal of Engineering,Vol.44,No.1:131-142,January 2022 https://doi.org/10.13374/j.issn2095-9389.2020.09.28.006;http://cje.ustb.edu.cn 巴西圆盘劈裂二维及三维数值模拟研究 乔兰”，刘建)，李庆文)@，赵国彦) 1)北京科技大学土木与资源工程学院.北京1000832)中南大学资源与安全工程学院.长沙410083 区通信作者，E-mail:qingwenli@ustb.edu.cn 摘要为揭示巴西圆盘起裂模式的变化规律及其破裂演化过程，运用连续介质弹塑性分析开展巴西圆盘劈裂二维及三维 数值模拟研究.通过开展二维模拟研究，探究压拉比及加载角对试样起裂破坏模式的影响：通过三维模拟研究，探究圆盘试 样三维破裂面的形成及扩展过程.二维数值模拟结果表明，接触加载角及压拉比越大，巴西圆盘试样越容易发生中心起裂： 端部起裂由剪切破坏引起，而劈裂裂纹进一步扩展则由张拉破坏驱动.三维数值模拟结果表明，初始起裂点位于三维圆盘端 面，随加载角增大其逐渐向端面圆心移动；当圆盘发生端面中心起裂时，三维破裂面以弧形边界向试样内部发散扩展.无论 圆盘试样发生中心起裂还是端部起裂，由于三维效应巴西劈裂试验可能都会低估岩石的抗拉强度. 关键词巴西劈裂：数值模拟：压拉比：加载角：三维破裂面 分类号TU452 Numerical study of the Brazilian tensile test:2D and 3D simulations QIAO Lan,LIU Jian,LI Qing-wen.ZHAO Guo-yan 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)School of Resources and Safety Engineering,Central South University,Changsha 410083,China Corresponding author,E-mail:qingwenli@ustb.edu.cn ABSTRACT The Brazilian splitting test is widely used to determine the tensile strength of rocks and rock-like materials due to its easy sample preparation and an easier compressive test setup as an indirect testing method compared with performing a direct uniaxial tensile test.However,the accuracy of this method has also been criticized for a long time in the literature since its introduction.This paper carried out two-dimensional (2D)/three-dimensional (3D)numerical simulations of the Brazilian tensile test using a continuum elastoplastic analysis to reveal the variation of fracture modes of the Brazilian disk and its fracture evolution process.The effect of compression-tension ratios and contact loading angles on the fracture modes of the disk specimens was studied through 2D simulations. Through 3D simulations,the initiation and expansion processes of the 3D fracture under different loading angles were explored.The simulated results of failure modes,stress distributions,and calculated tensile strengths were analyzed.The 2D numerical results show that the larger the contact loading angle and the compression-tension ratio,the more likely the Brazilian disk specimens crack first at the disk center.The fracture initiation under the loading rims is caused by shear failure,but further propagation of the split fracture is driven by tension failure.The 3D numerical simulation results show that the crack initiation point is always located on the end face of the disk and gradually moves to the center from the loading ends as the loading angle increases.When the central tensile cracking appears,the 3D fracture expanded toward the inside of the specimen with an arc boundary.Regardless of whether the disk specimen starts to fracture initially at the disk center or the loading points,the Brazilian tensile test may underestimate the tensile strength of rocks due to the 3D effect. KEY WORDS Brazilian test;numerical simulation;compression-tension ratio;loading angle;3D fracture 收稿日期：2020-09-28 基金项目：国家自然科学基金委与山东联合基金重点资助项目(U1806209):北京科技大学青年教师学科交叉研究资助项目(FRF-IDRY- 19-002):中央高校基本科研业务费专项资金资助项目(FRF-TP.19-021A3)

巴西圆盘劈裂二维及三维数值模拟研究 乔    兰1)，刘    建1)，李庆文1) 苣，赵国彦2) 1) 北京科技大学土木与资源工程学院，北京 100083    2) 中南大学资源与安全工程学院，长沙 410083 苣通信作者， E-mail: qingwenli@ustb.edu.cn 摘    要    为揭示巴西圆盘起裂模式的变化规律及其破裂演化过程，运用连续介质弹塑性分析开展巴西圆盘劈裂二维及三维 数值模拟研究. 通过开展二维模拟研究，探究压拉比及加载角对试样起裂破坏模式的影响；通过三维模拟研究，探究圆盘试 样三维破裂面的形成及扩展过程. 二维数值模拟结果表明，接触加载角及压拉比越大，巴西圆盘试样越容易发生中心起裂； 端部起裂由剪切破坏引起，而劈裂裂纹进一步扩展则由张拉破坏驱动. 三维数值模拟结果表明，初始起裂点位于三维圆盘端 面，随加载角增大其逐渐向端面圆心移动；当圆盘发生端面中心起裂时，三维破裂面以弧形边界向试样内部发散扩展. 无论 圆盘试样发生中心起裂还是端部起裂，由于三维效应巴西劈裂试验可能都会低估岩石的抗拉强度. 关键词    巴西劈裂；数值模拟；压拉比；加载角；三维破裂面 分类号    TU452 Numerical study of the Brazilian tensile test: 2D and 3D simulations QIAO Lan1) ，LIU Jian1) ，LI Qing-wen1) 苣 ，ZHAO Guo-yan2) 1) School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China 2) School of Resources and Safety Engineering, Central South University, Changsha 410083, China 苣 Corresponding author, E-mail: qingwenli@ustb.edu.cn ABSTRACT    The Brazilian splitting test is widely used to determine the tensile strength of rocks and rock-like materials due to its easy sample preparation and an easier compressive test setup as an indirect testing method compared with performing a direct uniaxial tensile test. However, the accuracy of this method has also been criticized for a long time in the literature since its introduction. This paper carried  out  two-dimensional  (2D)/three-dimensional  (3D)  numerical  simulations  of  the  Brazilian  tensile  test  using  a  continuum elastoplastic  analysis  to  reveal  the  variation  of  fracture  modes  of  the  Brazilian  disk  and  its  fracture  evolution  process.  The  effect  of compression-tension ratios and contact loading angles on the fracture modes of the disk specimens was studied through 2D simulations. Through 3D simulations, the initiation and expansion processes of the 3D fracture under different loading angles were explored. The simulated results of failure modes, stress distributions, and calculated tensile strengths were analyzed. The 2D numerical results show that the larger the contact loading angle and the compression–tension ratio, the more likely the Brazilian disk specimens crack first at the disk center. The fracture initiation under the loading rims is caused by shear failure, but further propagation of the split fracture is driven by tension failure. The 3D numerical simulation results show that the crack initiation point is always located on the end face of the disk and gradually moves to the center from the loading ends as the loading angle increases. When the central tensile cracking appears, the 3D fracture expanded toward the inside of the specimen with an arc boundary. Regardless of whether the disk specimen starts to fracture initially at the disk center or the loading points, the Brazilian tensile test may underestimate the tensile strength of rocks due to the 3D effect. KEY WORDS    Brazilian test；numerical simulation；compression–tension ratio；loading angle；3D fracture 收稿日期: 2020−09−28 基金项目: 国家自然科学基金委与山东联合基金重点资助项目（U1806209）；北京科技大学青年教师学科交叉研究资助项目（FRF-IDRY- 19-002）；中央高校基本科研业务费专项资金资助项目（FRF-TP-19-021A3） 工程科学学报，第 44 卷，第 1 期：131−142，2022 年 1 月 Chinese Journal of Engineering, Vol. 44, No. 1: 131−142, January 2022 https://doi.org/10.13374/j.issn2095-9389.2020.09.28.006; http://cje.ustb.edu.cn

·132 工程科学学报，第44卷，第1期 岩石、混凝土和陶瓷等脆性材料的抗拉强度 进行连续介质弹塑性分析，一方面开展巴西劈裂 远低于其抗压强度及抗剪强度，岩石的宏观剪切 二维模拟研究，探究圆盘试样压拉比及加载角对 破坏一般源于其内部大量张拉裂纹的相互作用及 试样起裂破坏模式的影响；另一方面开展三维模 贯通-习，Hoek和Martin)也曾指出深部岩体的破 拟研究，探究圆盘试样三维破裂面的形成及扩展 坏主要是应力驱动的张拉破坏，因此岩石类脆性 过程 材料的抗拉强度是人们尤其关注的重要参数.相 1 巴西圆盘应力解析解 对于直接拉伸试验，以巴西劈裂试验为代表的抗 拉强度间接测试方法因其操作简单方便而在国内 Hondros!2o推导出平面问题中弧形均布载荷 外得到广泛应用.国际岩石力学学会(ISRM于 作用下巴西圆盘内部任意一点应力分布的解析 1978年将巴西劈裂试验推荐为测定岩石材料抗拉 解，如图1所示，则加载直径上的应力分布为（以 强度的建议方法.与此同时，自巴西劈裂试验方 拉应力为正)： 法于1943年被提出以来，关于其有效性的探讨一 P 直是国内外学者的研究热点刀根据Griffith强度 理论，只有在满足中心起裂的条件下通过巴西劈 裂试验测得的岩石抗拉强度才是有效的阿，并且从 破坏结果来看其破裂面至少应表现为与张拉破坏 对应的单一径向破裂面.但是，Fairhurst指出对 于压拉比较低的脆性材料采用小加载角时起裂点 R 可能由于应力集中效应而位于岩样端部.Erarslan 0☐ 等⑧-)、Komurlu等o-采用ISRM推荐的弧形压 模开展巴西劈裂试验，发现圆盘试样发生多破裂 压碎破坏.Jaeger和Hoskins2提出采用弧形均布 T 载荷开展巴西劈裂试验，通过增大加载角以减小 图1巴西圆盘弧形均布载荷加载示意图 Fig.I A disk subjected to diametrically distributed loads in the 加载端的应力集中程度从而促进中心起裂 Brazilian test Erarslan等⑧-列认为加载角为20°或30时可以获得 可靠的抗拉强度.对于混凝土巴西圆盘试样，Garcia 2P 1-(r/R)2 sin2a 0r= 等]建议加载角应大于12°.Gutierrez-Moizant等4 πDta 1-2(r/R)2cos2o+(r/R)+ 则认为开展巴西劈裂试验最合适加载角为20° tan-11+(r/R)2 1-(r/R)2 tana (1) Bahaaddini等]通过开展数值模拟研究发现，对于 岩石类材料加载角为30°或更高时才能确保中心 2P [1-(r/R)2 sin2a 起裂.目前，鲜少有学者探究压拉比对试样起裂破 0g= Dta1-2(r/R)2cos2a+(r/R)4 坏模式的影响.此外，国内外学者比较关注圆盘试 样受压直径上的起裂点位置及破裂发展方式，较 1+(r/R)2 tan tang (2) 少有学者探究巴西劈裂三维破裂面的起裂及扩展 1-(r/R)2 过程.由于张拉破裂面沿受压直径向圆盘内部扩 Tre=0 (3) 展，受限于目前的观测手段，以高速摄影仪6及数 式中：Ors Od,Tre分别为受压直径上的径向应力、周 字图像相关技术刀为代表的表面变形观测技术便 向应力及剪应力，P为施加的外部载荷，D为巴西 不再适用，而数值模拟方法便显示出其独特的优 圆盘直径，R为圆盘半径，1为圆盘厚度，2a为加载 势.喻勇向、Li和Wong!7及邓华锋等u8通过开展 角度，r为受压直径上一点到圆心的距离.在满足 三维弹性分析探讨圆盘内部的应力分布规律及起 中心起裂的基础上，岩石抗拉强度的计算公式为 裂点位置，但未模拟三维破裂面的形成过程.许学 2p sin2a-1 2P πDta }D (4) 良运用PFCD开展巴西劈裂三维模拟，以细观 裂纹集聚表征张拉破裂面的扩展过程，模拟结果 当接触加载角小于15时，采用近似公式计算 中破裂面仍然从加载端开始起裂，并且起裂阶段初 岩石抗拉强度引起的误差不大于2.3%.此外，ISRM 始裂纹全部为张拉裂纹.鉴于此，本文借助FLAC3D 推荐的弧形压模与巴西圆盘接触在试样破坏时产

岩石、混凝土和陶瓷等脆性材料的抗拉强度 远低于其抗压强度及抗剪强度，岩石的宏观剪切 破坏一般源于其内部大量张拉裂纹的相互作用及 贯通[1−2] ，Hoek 和 Martin[3] 也曾指出深部岩体的破 坏主要是应力驱动的张拉破坏，因此岩石类脆性 材料的抗拉强度是人们尤其关注的重要参数. 相 对于直接拉伸试验，以巴西劈裂试验为代表的抗 拉强度间接测试方法因其操作简单方便而在国内 外得到广泛应用. 国际岩石力学学会 (ISRM) 于 1978 年将巴西劈裂试验推荐为测定岩石材料抗拉 强度的建议方法[4] . 与此同时，自巴西劈裂试验方 法于 1943 年被提出以来，关于其有效性的探讨一 直是国内外学者的研究热点[5−7] . 根据 Griffith 强度 理论，只有在满足中心起裂的条件下通过巴西劈 裂试验测得的岩石抗拉强度才是有效的[5] ，并且从 破坏结果来看其破裂面至少应表现为与张拉破坏 对应的单一径向破裂面. 但是，Fairhurst[5] 指出对 于压拉比较低的脆性材料采用小加载角时起裂点 可能由于应力集中效应而位于岩样端部. Erarslan 等[8−9]、Komurlu 等[10−11] 采用 ISRM 推荐的弧形压 模开展巴西劈裂试验，发现圆盘试样发生多破裂 压碎破坏. Jaeger 和 Hoskins[12] 提出采用弧形均布 载荷开展巴西劈裂试验，通过增大加载角以减小 加 载 端 的 应 力 集 中 程 度 从 而 促 进 中 心 起 裂 . Erarslan 等[8−9] 认为加载角为 20°或 30°时可以获得 可靠的抗拉强度. 对于混凝土巴西圆盘试样，García 等[13] 建议加载角应大于 12°. Gutiérrez-Moizant 等[14] 则认为开展巴西劈裂试验最合适加载角为 20°. Bahaaddini 等[15] 通过开展数值模拟研究发现，对于 岩石类材料加载角为 30°或更高时才能确保中心 起裂. 目前，鲜少有学者探究压拉比对试样起裂破 坏模式的影响. 此外，国内外学者比较关注圆盘试 样受压直径上的起裂点位置及破裂发展方式，较 少有学者探究巴西劈裂三维破裂面的起裂及扩展 过程. 由于张拉破裂面沿受压直径向圆盘内部扩 展，受限于目前的观测手段，以高速摄影仪[16] 及数 字图像相关技术[17] 为代表的表面变形观测技术便 不再适用，而数值模拟方法便显示出其独特的优 势. 喻勇[6]、Li 和 Wong[7] 及邓华锋等[18]通过开展 三维弹性分析探讨圆盘内部的应力分布规律及起 裂点位置，但未模拟三维破裂面的形成过程. 许学 良[19] 运用 PFC3D 开展巴西劈裂三维模拟，以细观 裂纹集聚表征张拉破裂面的扩展过程，模拟结果 中破裂面仍然从加载端开始起裂，并且起裂阶段初 始裂纹全部为张拉裂纹. 鉴于此，本文借助 FLAC3D 进行连续介质弹塑性分析，一方面开展巴西劈裂 二维模拟研究，探究圆盘试样压拉比及加载角对 试样起裂破坏模式的影响；另一方面开展三维模 拟研究，探究圆盘试样三维破裂面的形成及扩展 过程. 1    巴西圆盘应力解析解 Hondros[20] 推导出平面问题中弧形均布载荷 作用下巴西圆盘内部任意一点应力分布的解析 解，如图 1 所示，则加载直径上的应力分布为（以 拉应力为正）： P R σr σθ 2α 图 1    巴西圆盘弧形均布载荷加载示意图 Fig.1     A  disk  subjected  to  diametrically  distributed  loads  in  the Brazilian test σr =− 2P πDtα    [ 1−(r/R) 2 ] sin 2α 1−2(r/R) 2 cos 2α+(r/R) 4 + tan−1 [ 1+(r/R) 2 1−(r/R) 2 tanα ]} （1） σθ = 2P πDtα    [ 1−(r/R) 2 ] sin 2α 1−2(r/R) 2 cos 2α+(r/R) 4 − tan−1 [ 1+(r/R) 2 1−(r/R) 2 tanα ]} （2） τrθ = 0 （3） 式中：σr , σθ , τrθ 分别为受压直径上的径向应力、周 向应力及剪应力，P 为施加的外部载荷，D 为巴西 圆盘直径，R 为圆盘半径，t 为圆盘厚度，2α 为加载 角度，r 为受压直径上一点到圆心的距离. 在满足 中心起裂的基础上，岩石抗拉强度的计算公式为 σt = 2P πDt { sin 2α α −1 } ≈ 2P πDt （4） 当接触加载角小于 15°时，采用近似公式计算 岩石抗拉强度引起的误差不大于 2.3%. 此外，ISRM 推荐的弧形压模与巴西圆盘接触在试样破坏时产 · 132 · 工程科学学报，第 44 卷，第 1 期

乔兰等：巴西圆盘劈裂二维及三维数值模拟研究 133· 生的加载角为10°左右 及塑性拉应变累计至临界值后，单元内聚力及抗 拉强度弱化至残余强度并保持不变.为探究圆盘 2二维数值模拟 试样压拉比及加载角对试样起裂破坏模式的影 二维巴西劈裂数值模拟采用平面应变模型， 响，分别考虑6种压拉比()及加载角(2a):=20、 数值试样直径为50mm,厚度为0.5mm.模型网格 18、16、14、12和10及2a=4.58°、6.88°、9.18°、11.48°、 采用六面体立方体单元，单元尺寸为0.5mm,整个 13.78°和16.10°.针对每种压拉比，保持单轴抗压 模型共包含7864个单元，16130个节点.加载过程 强度为200MPa不变.数值模拟结果如表2所示. 采用位移控制，边界条件为：固定模型y方向的变 形，在模型顶端和底端施加：方向的速度边界条 件，加载速率为每步1×100m.二维加载示意图如 图2所示.采用弹塑性应变软化本构模型及Mohr- Coulomb强度准则，参考Iddefjord花岗岩的力 学参数，数值模型参数如表1所示.当模型单元发 生屈服后，单元内聚力及抗拉强度随塑性剪应变 图2二维巴西劈裂加载示意图 及塑性拉应变累计逐渐线性软化；当塑性剪应变 Fig.2 Loading setup of the 2D numerical Brazilian test 表1数值模型参数 Table 1 Material properties of the numerical model Elastic modulus/GPa Poisson's ratio Compressive strength/MPa Tensile strength/MPa Cohesion/MPa 60 0.25 200 38.39 Friction angle/() Residue cohesion/MPa Residue tensile strength/MPa Plastic shear strain limit Plastic tension strain limit 48 1.92 0.4 5×10 2×10 表2不同压拉比及加载角巴西试样的破坏模式 Table 2 Failure modes of the numerical Brazilian disks with different compression-tension ratios()and contact loading angles(2a) Failure modes 2a=16.10° 2a=13.78° 2a=11.48° 2a=9.18° 2a=6.88° 2a=4.58 20

生的加载角为 10°左右[4] . 2    二维数值模拟 二维巴西劈裂数值模拟采用平面应变模型， 数值试样直径为 50 mm，厚度为 0.5 mm. 模型网格 采用六面体立方体单元，单元尺寸为 0.5 mm，整个 模型共包含 7864 个单元，16130 个节点. 加载过程 采用位移控制，边界条件为：固定模型 y 方向的变 形，在模型顶端和底端施加 z 方向的速度边界条 件，加载速率为每步 1×10−10 m. 二维加载示意图如 图 2 所示. 采用弹塑性应变软化本构模型及 Mohr￾Coulomb 强度准则，参考 Iddefjord 花岗岩[21] 的力 学参数，数值模型参数如表 1 所示. 当模型单元发 生屈服后，单元内聚力及抗拉强度随塑性剪应变 及塑性拉应变累计逐渐线性软化；当塑性剪应变 及塑性拉应变累计至临界值后，单元内聚力及抗 拉强度弱化至残余强度并保持不变. 为探究圆盘 试样压拉比及加载角对试样起裂破坏模式的影 响，分别考虑 6 种压拉比 (λ) 及加载角 (2α)：λ=20、 18、16、14、12 和10 及2α=4.58°、6.88°、9.18°、11.48°、 13.78°和 16.10°. 针对每种压拉比，保持单轴抗压 强度为 200 MPa 不变. 数值模拟结果如表 2 所示. 2α 图 2    二维巴西劈裂加载示意图 Fig.2    Loading setup of the 2D numerical Brazilian test 表 1 数值模型参数 Table 1   Material properties of the numerical model Elastic modulus/GPa Poisson’s ratio Compressive strength/MPa Tensile strength/MPa Cohesion/MPa 60 0.25 200 8 38.39 Friction angle/(°) Residue cohesion/MPa Residue tensile strength/MPa Plastic shear strain limit Plastic tension strain limit 48 1.92 0.4 5 × 10−4 2 × 10−4 表 2 不同压拉比及加载角巴西试样的破坏模式 Table 2   Failure modes of the numerical Brazilian disks with different compression−tension ratios (λ) and contact loading angles (2α) λ Failure modes 2α = 16.10° 2α = 13.78° 2α = 11.48° 2α = 9.18° 2α = 6.88° 2α = 4.58° 20 18 16 14 12 10 乔    兰等： 巴西圆盘劈裂二维及三维数值模拟研究 · 133 ·

134 工程科学学报，第44卷，第1期 从表2横向分析，接触加载角越大，巴西圆盘 鹏等关于页岩巴西劈裂的声发射定位结果类 试样越容易发生中心起裂，因为增大加载角可以 似.值得指出的是，无论是均匀试样还是非均匀试 减小圆盘试样顶底端的应力集中程度.以=14 样，端部起裂均由剪切破坏引起，而劈裂裂纹进一 为例，当2a≤9.18时，巴西试样发生端部起裂： 步扩展则由张拉破坏驱动 当2a≥11.48时，巴西试样发生中心起裂.沿表2 3三维数值模拟 纵向分析，压拉比越大，巴西圆盘试样越容易发 生中心起裂.以2a=11.48°为例，当1≥14时，试样 三维巴西圆盘数值模型直径为50mm,厚度 发生中心起裂；而当1≤12时，圆盘试样发生端 为25mm;采用六面体网格进行划分，对受压直径 部起裂.较大的压拉比，一方面表示岩样的抗拉 附近的网格进行适当加密，模型共包含69125个单 强度相对较低，因此试样更容易发生张拉破坏； 元，74360个节点，模型网格如图4所示.采用应力 另一方面则表示岩样的抗压强度相对较高，抵抗 边界条件，在受压区域内沿模型表面施加法向正 剪切破坏的能力强.因此，巴西圆盘试样的起裂 应力，加载速率为每步500Pa.三维数值模拟仍采 破坏模式受内部因素与外部因素共同控制，总体 用表1中的力学参数；为模拟不同接触加载角下 上，圆盘试样压拉比越大，外部加载角越大，试样 巴西圆盘三维破裂面的起裂及扩展过程，考虑以 破坏越容易产生中心起裂，反之则越容易发生端 下8种加载角度：2a=29.7°、25.7°、21.7°、17.9°、14.4°、 部起裂 10.9°、7.6°、4.4° 图3(a)~(b)为均匀巴西圆盘试样的加载破 图5为加载角2a=29.7°、25.7°、14.4°和4.4°时 坏过程，包括峰值载荷(Fpak)前及峰后阶段，接 巴西圆盘试样三维破裂面的起裂及扩展过程，由 触加载角为9.18°，压拉比分别为1=20、10.当 于篇幅限制，其他4种加载角下的数值模拟结果 =20时，圆盘试样发生中心起裂，随外载进一步增 没有给出，但破裂面发展规律基本类似.当2a= 加，径向张拉裂纹自中心起裂点逐渐向两端扩展， 29.7时，初始起裂点位于圆盘端面中心，随外部载 并在加载点附近趋于停止，这与Erarslan等8通 荷增加，破裂面以弧形边界向试样内部扩展，在圆 过实验观察到的裂纹扩展模式相似，陈沙等P21也 盘厚度方向首先贯通后继续沿受压直径向两端扩 通过数值模拟获得过类似的裂纹扩展过程.而当 展并在加载点附近趋于停止.值得注意的是，贯通 =10时，两条起裂裂纹从圆盘受力区域边缘向圆 后的破裂面沿受压直径方向的扩展边界也并非直 盘内部扩展，初始起裂裂纹主要发生剪切破坏，向 线，这与平面应变假设下贯通破裂面的直线扩展 试样中心扩展后则转变为张拉破坏，最终形成两 情形相矛盾.当2=25.7°、14.4时，起裂点仍位于 条贯通的劈裂缝，十分明显这是由于各向同性圆盘 圆盘端面受压直径上，但却偏离端面中心，同时出 试样及外载关于受压直径的对称性引起.图3(c)~ 现4个对称的起裂点：加载角越小，起裂点越靠近 (ε)为非均匀巴西圆盘试样的加载破坏过程，加载 上下加载端，并且破裂面扩展边界线仍为弧形 角仍为9.18°，压拉比分别为=20、10、8.运用刘 当2=4.4时，起裂点位于圆盘端面加载区域边缘， 建等2)提出的岩石介质细观物理力学参数统计分 即破裂没有在受压中心面上发展，这与二维模拟 布模型，采用对数正太分布描述细观单元弹性模 中端部起裂的情况类似.事实上，8种接触加载角 量及单轴抗压强度的空间变异性24-2，试样均质 下圆盘试样的起裂点均位于圆盘端面，而不是圆 度参数设为0.3.当=20时，数值试样同样产生中 盘内部中心.喻勇网及邓华锋等P0通过对线荷载 心起裂，受模型非均质性的影响起裂点偏离圆盘 作用下的巴西圆盘试样进行三维弹性分析指出， 正中心，同时在圆盘内部受压直径周围分布有离 圆盘横截面上的应力值沿试样厚度方向是有变化 散的张拉破裂单元.这是由于一方面材料的非均 的，越靠近两端面，水平拉应力越大，因此试样破 质性改变了圆盘内部的应力分布，引起局部应力 坏首先从端面开始 集中；另一方面，可能是由于软弱介质单元的过早 4解析验证及应力分析 破坏引起.而当=10、8时，与均匀条件下的模拟 结果明显不同的是，圆盘试样发生端部起裂后最 4.1二维模拟结果分析 终仅形成一条劈裂缝；在劈裂缝形成之前，试样两 图6为=20时不同加载角下二维巴西圆盘受 端出现剪切破坏单元，而试样内部基本无破坏单 压直径上切向应力与径向应力的数值模拟结果与 元出现，压拉比为8的数值试样更明显，这与杨志 理论计算结果的对比示意图.由图可知，除在加载

从表 2 横向分析，接触加载角越大，巴西圆盘 试样越容易发生中心起裂，因为增大加载角可以 减小圆盘试样顶底端的应力集中程度. 以 λ=14 为例，当 2α ≤ 9.18°时，巴西试样发生端部起裂； 当 2α ≥ 11.48°时，巴西试样发生中心起裂. 沿表 2 纵向分析，压拉比越大，巴西圆盘试样越容易发 生中心起裂. 以 2α=11.48°为例，当 λ ≥ 14 时，试样 发生中心起裂；而当 λ ≤ 12 时，圆盘试样发生端 部起裂. 较大的压拉比，一方面表示岩样的抗拉 强度相对较低，因此试样更容易发生张拉破坏； 另一方面则表示岩样的抗压强度相对较高，抵抗 剪切破坏的能力强. 因此，巴西圆盘试样的起裂 破坏模式受内部因素与外部因素共同控制，总体 上，圆盘试样压拉比越大，外部加载角越大，试样 破坏越容易产生中心起裂，反之则越容易发生端 部起裂. 图 3（a）～（b）为均匀巴西圆盘试样的加载破 坏过程，包括峰值载荷 (Fpeak) 前及峰后阶段，接 触加载角 为 9.18°，压拉比分别 为 λ=20、 10. 当 λ=20 时，圆盘试样发生中心起裂，随外载进一步增 加，径向张拉裂纹自中心起裂点逐渐向两端扩展， 并在加载点附近趋于停止，这与 Erarslan 等[8−9] 通 过实验观察到的裂纹扩展模式相似，陈沙等[22] 也 通过数值模拟获得过类似的裂纹扩展过程. 而当 λ=10 时，两条起裂裂纹从圆盘受力区域边缘向圆 盘内部扩展，初始起裂裂纹主要发生剪切破坏，向 试样中心扩展后则转变为张拉破坏，最终形成两 条贯通的劈裂缝，十分明显这是由于各向同性圆盘 试样及外载关于受压直径的对称性引起. 图 3（c）～ （e）为非均匀巴西圆盘试样的加载破坏过程，加载 角仍为 9.18°，压拉比分别为 λ=20、10、8. 运用刘 建等[23] 提出的岩石介质细观物理力学参数统计分 布模型，采用对数正太分布描述细观单元弹性模 量及单轴抗压强度的空间变异性[24−25] ，试样均质 度参数设为 0.3. 当 λ=20 时，数值试样同样产生中 心起裂，受模型非均质性的影响起裂点偏离圆盘 正中心，同时在圆盘内部受压直径周围分布有离 散的张拉破裂单元. 这是由于一方面材料的非均 质性改变了圆盘内部的应力分布，引起局部应力 集中；另一方面，可能是由于软弱介质单元的过早 破坏引起. 而当 λ=10、8 时，与均匀条件下的模拟 结果明显不同的是，圆盘试样发生端部起裂后最 终仅形成一条劈裂缝；在劈裂缝形成之前，试样两 端出现剪切破坏单元，而试样内部基本无破坏单 元出现，压拉比为 8 的数值试样更明显，这与杨志 鹏等[26] 关于页岩巴西劈裂的声发射定位结果类 似. 值得指出的是，无论是均匀试样还是非均匀试 样，端部起裂均由剪切破坏引起，而劈裂裂纹进一 步扩展则由张拉破坏驱动. 3    三维数值模拟 三维巴西圆盘数值模型直径为 50 mm，厚度 为 25 mm；采用六面体网格进行划分，对受压直径 附近的网格进行适当加密，模型共包含 69125 个单 元，74360 个节点，模型网格如图 4 所示. 采用应力 边界条件，在受压区域内沿模型表面施加法向正 应力，加载速率为每步 500 Pa. 三维数值模拟仍采 用表 1 中的力学参数；为模拟不同接触加载角下 巴西圆盘三维破裂面的起裂及扩展过程，考虑以 下8 种加载角度：2α=29.7°、25.7°、21.7°、17.9°、14.4°、 10.9°、7.6°、4.4°. 图 5 为加载角 2α=29.7°、25.7°、14.4°和 4.4°时 巴西圆盘试样三维破裂面的起裂及扩展过程，由 于篇幅限制，其他 4 种加载角下的数值模拟结果 没有给出，但破裂面发展规律基本类似. 当 2α= 29.7°时，初始起裂点位于圆盘端面中心，随外部载 荷增加，破裂面以弧形边界向试样内部扩展，在圆 盘厚度方向首先贯通后继续沿受压直径向两端扩 展并在加载点附近趋于停止. 值得注意的是，贯通 后的破裂面沿受压直径方向的扩展边界也并非直 线，这与平面应变假设下贯通破裂面的直线扩展 情形相矛盾. 当 2α=25.7°、14.4°时，起裂点仍位于 圆盘端面受压直径上，但却偏离端面中心，同时出 现 4 个对称的起裂点；加载角越小，起裂点越靠近 上下加载端，并且破裂面扩展边界线仍为弧形. 当 2α=4.4°时，起裂点位于圆盘端面加载区域边缘， 即破裂没有在受压中心面上发展，这与二维模拟 中端部起裂的情况类似. 事实上，8 种接触加载角 下圆盘试样的起裂点均位于圆盘端面，而不是圆 盘内部中心. 喻勇[8] 及邓华锋等[20] 通过对线荷载 作用下的巴西圆盘试样进行三维弹性分析指出， 圆盘横截面上的应力值沿试样厚度方向是有变化 的，越靠近两端面，水平拉应力越大，因此试样破 坏首先从端面开始. 4    解析验证及应力分析 4.1    二维模拟结果分析 图 6 为 λ=20 时不同加载角下二维巴西圆盘受 压直径上切向应力与径向应力的数值模拟结果与 理论计算结果的对比示意图. 由图可知，除在加载 · 134 · 工程科学学报，第 44 卷，第 1 期

乔兰等：巴西圆盘劈裂二维及三维数值模拟研究 135· 83%0Fpek gocoFpea 99.29%Fet 100%F 85%Fe或 (a) 83% 100% 98%Fk 76% 68% (b) 99%F 1009%Fe (c) 95%oFpeak 100%Fpea 989%Fe* 92%oF peak 81% (d) 97% 100%Fe 96%Fpe 86%Fa 59% (e) □None Shear-n,tension-n,shear-p,tension-p Tension-n,shear-p,tension-p Shear-n,shear-p Shear-p Tension-n,tension-p Shear-n,shear-p,tension-p Shear-p,tension-p Tension-p 图3不同压拉比时巴西圆盘试样的破坏过程.(a)均匀试样-20：(b)均匀试样=10：(c)非均匀试样=20：(d)非均匀试样=10：(e)非均匀试样 1=8 Fig.3 Failure processes of the numerical disks with different compression-tension ratios (under the loading angle of 9.18:(a)homogeneous disk= 20.(b)homogeneous disk=10.(c)heterogeneous disk=20.(d)heterogeneous disk=10.(e)heterogeneous disk=8 端附近数值模拟结果与解析解存在较小误差外， 网格划分，因此在模型顶底端为直线边界而非弧 受压直径上大部分区域具有较好的吻合性，说明 形边界，并且施加的速度边界条件与受压直径平 了二维数值模拟的有效性.需要指出的是，由于二 行而非沿圆盘半径指向圆盘中心，如图3所示，这 维数值模型采用正六面体单元逼近圆盘试样进行 是造成上述误差的原因.在加载端附近，径向应力

端附近数值模拟结果与解析解存在较小误差外， 受压直径上大部分区域具有较好的吻合性，说明 了二维数值模拟的有效性. 需要指出的是，由于二 维数值模型采用正六面体单元逼近圆盘试样进行 网格划分，因此在模型顶底端为直线边界而非弧 形边界，并且施加的速度边界条件与受压直径平 行而非沿圆盘半径指向圆盘中心，如图 3 所示，这 是造成上述误差的原因. 在加载端附近，径向应力 83%Fpeak 100%Fpeak 98%Fpeak 76%Fpeak 68%Fpeak (b) 94%Fpeak 98%Fpeak 99%Fpeak 100%Fpeak 84%Fpeak (c) 95%Fpeak 100%Fpeak 98%Fpeak 92%Fpeak 81%Fpeak (d) 97%Fpeak 100%Fpeak 96%Fpeak 86%Fpeak 59%Fpeak (e) 99%Fpeak 99.2%Fpeak 100%Fpeak 85%Fpeak (a) 83%Fpeak None Shear-n, shear-p Shear-n, shear-p, tension-p Shear-n, tension-n, shear-p, tension-p Shear-p Shear-p, tension-p Tension-n, shear-p, tension-p Tension-n, tension-p Tension-p 图 3    不同压拉比时巴西圆盘试样的破坏过程. （a）均匀试样 λ=20；（b） 均匀试样 λ=10；（c）非均匀试样 λ=20；（d）非均匀试样 λ=10；（e）非均匀试样 λ=8 Fig.3    Failure processes of the numerical disks with different compression–tension ratios (λ) under the loading angle of 9.18°: (a) homogeneous disk λ = 20; (b) homogeneous disk λ = 10; (c) heterogeneous disk λ = 20; (d) heterogeneous disk λ = 10; (e) heterogeneous disk λ = 8 乔    兰等： 巴西圆盘劈裂二维及三维数值模拟研究 · 135 ·

·136 工程科学学报，第44卷，第1期 与切向应力均为压应力，靠近圆盘中心，两者快速 减小，并且切向应力由压应力逐渐转变为拉应力 随着接触加载角增大，受压直径上张拉区域逐渐 减小，加载端附近压缩区域逐渐增大，径向与切向 压应力也逐渐减小，这说明端部应力集中程度在 减弱 表3为不同加载角及压拉比下二维巴西圆盘 试样抗拉强度的计算结果.由表3可知，对于发生 中心起裂的圆盘试样，通过公式(4)计算的抗拉强 图4三维巴西圆盘模型网格示意图 度与圆盘试样的真实抗拉强度比较接近，相对误 Fig.4 Grid model of the 3D numerical Brazilian disk 差在3%以内（括号内数据为相对误差）：而对于发 (a) (b) (d) 图5不同加载角时巴西圆盘试样三维破裂面的起裂及扩展过程.(a)2a=29.7°：(b)2a=25.7°：(c)2a=14.4°：(d)2a=4.4° Fig.5 Fracture initiation and propagation processes of the 3D Brazilian disks with different contact loading angles:(a)2a=29.7;(b)2a=25.7;(c)2a= 14.4°，(d)2a=4.4°

与切向应力均为压应力，靠近圆盘中心，两者快速 减小，并且切向应力由压应力逐渐转变为拉应力. 随着接触加载角增大，受压直径上张拉区域逐渐 减小，加载端附近压缩区域逐渐增大，径向与切向 压应力也逐渐减小，这说明端部应力集中程度在 减弱. 表 3 为不同加载角及压拉比下二维巴西圆盘 试样抗拉强度的计算结果. 由表 3 可知，对于发生 中心起裂的圆盘试样，通过公式（4）计算的抗拉强 度与圆盘试样的真实抗拉强度比较接近，相对误 差在 3% 以内（括号内数据为相对误差）；而对于发 Z X 图 4    三维巴西圆盘模型网格示意图 Fig.4    Grid model of the 3D numerical Brazilian disk Y Y Y Y Z X Y Z X Z Z Z Z Z Z X X X X X X Y Z X Y Z X Y Y Z X Z X Z X Y Z X Y Z X Y Z X (a) (b) (c) (d) 图 5    不同加载角时巴西圆盘试样三维破裂面的起裂及扩展过程. （a）2α=29.7°；（b）2α=25.7°；（c）2α=14.4°；（d）2α=4.4° Fig.5    Fracture initiation and propagation processes of the 3D Brazilian disks with different contact loading angles: (a) 2α = 29.7°; (b) 2α = 25.7°; (c) 2α = 14.4°; (d) 2α = 4.4° · 136 · 工程科学学报，第 44 卷，第 1 期

乔兰等：巴西圆盘劈裂二维及三维数值模拟研究 137 25 Numerical solution:o .Numerical solution: Hon山os'solu山on:Ge Hond山ros'solut山an:oe 15 Numerical solution:o Numerical solution:o Hondros'solution: Hondros'solution: 5 5 0 0 -15-2015-10-50 51015 2 -15-20-15-10-50 101520 -5 110002 0 10 -15 -20 -20 (a) (b) Principle stresses normalized by 2P/Dt Principle stresses normalized by 2P/aDt 20 Numerical solution: Numerical solution: Hond山0 s solution:aa ,Hod山os'solution:ca 15 Numerieal solution a Numerical solution:a 10 Hondros'solution: Hondros'solution: 10 5 0 0 -15-20-15-10-5 51015 -15-20-15-10-5 510 15 -10 10 -15 -15 -20 -20 (c) (d) Principle stresses normalized by 2P/aDt Principle stresses normalized by 2P/aDt 25 25 20 .Numerical solution: .Numerical solution: Hondros'solution: 20 Hondros'solution: 15 Numerical solution:o, 15 Numerical solution:o, Hondros'solution: -Hondros'solution: 10 5 0 0 -5-2015-10-5 510 15 -5-20F15-10-5 1015 20 -5 -5 10 10 20 55 (e) (0 Principle stresses normalized by 2P/Dr Principle stresses normalized by 2P/Dr 图6不同加载角时二维巴西圆盘受压直径上应力分布的数值模拟结果与理论计算结果对比示意图.(a)2a=4.58°：(b)2a=6.88°：(c)2a=9.18°： (d)2a=11.48°：(e)2a=13.78°：(f)2a=16.10° Fig.6 Comparison of normalized stresses along the compressed diameter between the numerical results and Hondros'solutions with different contact loading angles:(a)2a=4.58°，(b)2a=6.88，°(c)2a=9.18°，(d)2a=11.48°，(e)2a=13.78°，(02a=16.10° 生端部起裂的圆盘试样来说，其计算抗拉强度明 径上径向与切向应力的数值模拟结果与Hondros 显低于试样的真实抗拉强度，并且压拉比及加载 解析解的计算结果十分吻合，而端面受压直径上 角越小，两者相对误差也越大.上述分析表明只有 应力分布的数值模拟结果与Hondros解析结果 在满足中心起裂的基础上依据峰值荷载计算出的 存在一定误差，这是由于三维圆盘试样及外载关 抗拉强度才是有效的，而端部起裂条件下巴西劈 于中心剖面对称，因此中心剖面更贴近平面应变 裂实验会低估岩石的抗拉强度 条件.圆盘端面受压直径上切向拉应力大于中心 4.2三维模拟结果分析 剖面受压直径上的切向拉应力，2a=25.7°、29.7时 图7为接触加载角2c=4.4°、14.4°、25.7°、29.7° 尤其明显.图8为8种加载角下三维巴西圆盘中 时三维巴西圆盘试样端面(Y=0.025)及中心剖面 心剖面及端面受压直径上切向应力分布示意图， (Y=0.0125)受压直径上径向应力与切向应力的分 中心剖面受压直径上最大切向拉应力始终位于 布示意图，图中同时给出了Hondros解析解的计 剖面中心；而在端面受压直径上，当加载角比较 算结果.由图可知，4种加载角下中心剖面受压直 小时，最大切向拉应力靠近圆盘顶底端，随着加

生端部起裂的圆盘试样来说，其计算抗拉强度明 显低于试样的真实抗拉强度，并且压拉比及加载 角越小，两者相对误差也越大. 上述分析表明只有 在满足中心起裂的基础上依据峰值荷载计算出的 抗拉强度才是有效的，而端部起裂条件下巴西劈 裂实验会低估岩石的抗拉强度. 4.2    三维模拟结果分析 图 7 为接触加载角 2α=4.4°、14.4°、25.7°、29.7° 时三维巴西圆盘试样端面 (Y=0.025) 及中心剖面 (Y=0.0125) 受压直径上径向应力与切向应力的分 布示意图，图中同时给出了 Hondros 解析解的计 算结果. 由图可知，4 种加载角下中心剖面受压直 径上径向与切向应力的数值模拟结果与 Hondros 解析解的计算结果十分吻合，而端面受压直径上 应力分布的数值模拟结果与 Hondros 解析结果 存在一定误差，这是由于三维圆盘试样及外载关 于中心剖面对称，因此中心剖面更贴近平面应变 条件. 圆盘端面受压直径上切向拉应力大于中心 剖面受压直径上的切向拉应力，2α=25.7°、29.7°时 尤其明显. 图 8 为 8 种加载角下三维巴西圆盘中 心剖面及端面受压直径上切向应力分布示意图， 中心剖面受压直径上最大切向拉应力始终位于 剖面中心；而在端面受压直径上，当加载角比较 小时，最大切向拉应力靠近圆盘顶底端，随着加 Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr (a) (b) (c) (d) (e) (f) Z coordinate of radius/mm Z coordinate of radius/mm −25 20 −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 20 −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 20 −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 20 −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 20 −20 −15 −10 −5 0 5 10 15 −25 −20 −15 −10 −5 0 5 10 15 20 25 −25 20 Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Z coordinate of radius/mm Z coordinate of radius/mm Principle stresses normalized by 2P/πDt Z coordinate of radius/mm Z coordinate of radius/mm Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr Numerical solution: σθ Hondros’ solution: σθ Numerical solution: σr Hondros’ solution: σr 图 6    不同加载角时二维巴西圆盘受压直径上应力分布的数值模拟结果与理论计算结果对比示意图. （a） 2α=4.58°；（b）2α=6.88°；（c）2α=9.18°； （d）2α=11.48°；（e）2α=13.78°；（f）2α=16.10° Fig.6    Comparison of normalized stresses along the compressed diameter between the numerical results and Hondros ’ solutions with different contact loading angles: (a) 2α = 4.58°; (b) 2α = 6.88°; (c) 2α = 9.18°; (d) 2α = 11.48°; (e) 2α = 13.78°; (f) 2α = 16.10° 乔    兰等： 巴西圆盘劈裂二维及三维数值模拟研究 · 137 ·

·138 工程科学学报，第44卷，第1期 表3不同加载角及压拉比下二维巴西圆盘试样抗拉强度计算结果 Table 3 Calculated tensile strength of 2D numerical Brazilian disks with different contact loading angles and compression-tension ratios Calculated tensile strength/MPa Real strength/MPa 2a=16.10° 2a=13.78° 2a=11.48° 2a=9.18° 2a=6.88° 2a=4.58° 20 10.0 10.3(2.6%) 10.2(2.2%) 102(1.7%) 10.1(1.4%) 10.1(1.2%) 8.7(13.3%) 18 11.1 11.4(2.5%) 11.3(2.0%) 11.3(1.6%) 11.3(1.3%) 10.7(3.3%) 8.6(22.6%) 16 12.5 12.8(2.5%) 12.7(1.9%) 12.7(1.5%) 12.1(3.4%) 11.3(9.7%) 8.7(30.8%) 14 14.3 14.6(2.2%) 14.5(1.6%) 14.5(1.2%) 11.9(16.7%) 11.6(19.0%) 8.7(39.5%) 12 16.7 16.8(1.0%) 16.3(2.5%) 14.7(12.1%) 11.7(29.9%) 10.7(35.6%) 8.6(48.3%) 10 20.0 17.0(14.8%) 16.6(17.0%) 14.6(26.8%) 12.9(35.7%) 10.7(46.7%) 8.6(56.9%) Note:The data in parentheses are relative errors. 20 Surface: Surface: :!: ce:d, Surface:o, 15 Middle:o Middle:o. Middle:o. 10 Honrros: 10 Honrros: 50 Hondros: Hondros:o 0 -24 -16 -8 8 16 24 6 -12 -8 4 8 - 10 -10 Maximum tangehtial Maximum tangentia 15 tensile stress -15 tensile stress -20 (a) (b) Principle stresses nommalized by 2P/Dt Principle stresses normalized by 2P/tDt 2 ◆-Surface:oe Surface:o -Surface:o, Surface: 15 -Middle:o. Middle:o. 10 Honrros:0 10 Honrros:0 -Hondros:o, Hondros:o, 0 Q 0 4 2 0 6 -10 Maximum tangentia -10 tensile stress -15 Maximum taneentia -15 tensile 2 (c) (d) 15 2 Principle stresses normalized by 2P/Dt Principle stresses normalized by 2P/nDt 图7不同加载角下三维巴西圆盘受压直径上应力分布的数值模拟结果与理论计算结果对比示意图.(a)2a=4.4°：(b)2a=14.4°：(c)2a=25.7°： (d)2a=29.7° Fig.7 Comparison of normalized stresses along the compressed diameter between the numerical results and Hondros'solutions with different contact loading angles:(a)2a=4.4°：(b)2a=14.4°，(c)2a=25.7,(d2a=29.7 载角增大，最大切向拉应力逐渐向端面圆心移 越小，水平张拉变形越大，因此张拉应力就相对越 动，当2=29.7时，最大切向拉应力出现在端面圆 大，所以三维圆盘试样的最先起裂点出现在圆盘 心处 端面.而在圆盘端面上，除2a-4.4外，其他7种加 图9为2a=29.7°、25.7°、144°、7.6°、4.4时巴 载角下，起裂点位置就是端面受压直径上最大切 西圆盘轴向受压平面上(X=O)切向拉应力分布云 向拉应力出现的点.本文运用考虑张拉截断的 图.十分明显，圆盘内部切向拉应力明显小于圆盘 Mohr-Coulomb强度准则进行模拟分析，当圆盘端 外部切向拉应力，这是由于各质点受到周围质点 面受压直径上某一点的最大切向拉应力大于试样 的约束程度不一样引起的.从圆盘中心剖面(Y= 抗拉强度时就会发生张拉破坏，其本质为最大拉 0.0125)越靠近圆盘端面，质点受到的水平约束也 应力强度准则.Li和Wong!7.运用Staceyl27张拉应

载角增大，最大切向拉应力逐渐向端面圆心移 动，当 2α=29.7°时，最大切向拉应力出现在端面圆 心处. 图 9 为 2α=29.7°、25.7°、14.4°、7.6°、4.4°时巴 西圆盘轴向受压平面上 (X=0) 切向拉应力分布云 图. 十分明显，圆盘内部切向拉应力明显小于圆盘 外部切向拉应力，这是由于各质点受到周围质点 的约束程度不一样引起的. 从圆盘中心剖面 (Y= 0.0125) 越靠近圆盘端面，质点受到的水平约束也 越小，水平张拉变形越大，因此张拉应力就相对越 大，所以三维圆盘试样的最先起裂点出现在圆盘 端面. 而在圆盘端面上，除 2α=4.4°外，其他 7 种加 载角下，起裂点位置就是端面受压直径上最大切 向拉应力出现的点. 本文运用考虑张拉截断的 Mohr-Coulomb 强度准则进行模拟分析，当圆盘端 面受压直径上某一点的最大切向拉应力大于试样 抗拉强度时就会发生张拉破坏，其本质为最大拉 应力强度准则. Li 和 Wong[7] 运用 Stacey[27] 张拉应 Z coordinate of radius/mm Z coordinate of radius/mm Z coordinate of radius/mm Z coordinate of radius/mm (a) (b) (c) (d) −32 −24 −16 −8 0 8 16 24 −25 −20 −15 −10 −5 0 5 10 15 20 25 Maximum tangential tensile stress −16 −12 −8 −4 0 4 8 12 −25 −20 −15 −10 −5 0 5 10 15 20 25 Maximum tangential tensile stress −8 −6 −4 −2 0 2 4 6 −25 −20 −15 −10 −5 0 5 10 15 20 25 Maximum tangential tensile stress −8 −6 −4 −2 0 2 4 6 −25 −20 −15 −10 −5 0 5 10 15 20 25 Maximum tangential tensile stress Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Principle stresses normalized by 2P/πDt Hondros: σr Honrros: σθ Surface: σθ Surface: σr Middle: σθ Middle: σr Hondros: σr Honrros: σθ Surface: σθ Surface: σr Middle: σθ Middle: σr Hondros: σr Honrros: σθ Surface: σθ Surface: σr Middle: σθ Middle: σr Hondros: σr Honrros: σθ Surface: σθ Surface: σr Middle: σθ Middle: σr 图 7    不同加载角下三维巴西圆盘受压直径上应力分布的数值模拟结果与理论计算结果对比示意图. （a）2α=4.4°；（b）2α=14.4°；（c）2α=25.7°； （d）2α=29.7° Fig.7    Comparison of normalized stresses along the compressed diameter between the numerical results and Hondros ’ solutions with different contact loading angles: (a) 2α = 4.4°; (b) 2α = 14.4°; (c) 2α = 25.7°; (d) 2α = 29.7° 表 3 不同加载角及压拉比下二维巴西圆盘试样抗拉强度计算结果 Table 3 Calculated tensile strength of 2D numerical Brazilian disks with different contact loading angles and compression–tension ratios λ Real strength/MPa Calculated tensile strength/MPa 2α = 16.10° 2α = 13.78° 2α = 11.48° 2α = 9.18° 2α = 6.88° 2α = 4.58° 20 10.0 10.3 (2.6%) 10.2 (2.2%) 10.2 (1.7%) 10.1 (1.4%) 10.1 (1.2%) 8.7 (13.3%) 18 11.1 11.4 (2.5%) 11.3 (2.0%) 11.3 (1.6%) 11.3 (1.3%) 10.7 (3.3%) 8.6 (22.6%) 16 12.5 12.8 (2.5%) 12.7 (1.9%) 12.7 (1.5%) 12.1 (3.4%) 11.3 (9.7%) 8.7 (30.8%) 14 14.3 14.6 (2.2%) 14.5 (1.6%) 14.5 (1.2%) 11.9 (16.7%) 11.6 (19.0%) 8.7 (39.5%) 12 16.7 16.8 (1.0%) 16.3 (2.5%) 14.7 (12.1%) 11.7 (29.9%) 10.7 (35.6%) 8.6 (48.3%) 10 20.0 17.0 (14.8%) 16.6 (17.0%) 14.6 (26.8%) 12.9 (35.7%) 10.7 (46.7%) 8.6 (56.9%) Note: The data in parentheses are relative errors. · 138 · 工程科学学报，第 44 卷，第 1 期

乔兰等：巴西圆盘劈裂二维及三维数值模拟研究 139· (a) (b) 0 -20=4.4 1 20=760 2a=7.6° 10.9 sassans =109 =14.49 799 170 79 =29.7° Lutua =29.7° -25-20-15-10-505101520 色 25 5 -20-15-10-50510152025 Z coordinate of radius/mm Z coordinate of radius/mm 图8不同加载角下三维巴西圆盘试样中心剖面及端面受压直径上切向应力分布示意图.()中心剖面(=0.0125)切向应力：(b)端面(=0.025)切 向应力 Fig.8 Normalized tangential stresses along the compressed diameter of the surface and middle section of 3D disks with different contact loading angles: (a)middle section (Y=0.0125);(b)surface (Y=0.025) 03 0.6 .6 (a (b) (c) (d) ⑧ 国9不同加载角时三维巴西圆盘试样轴向受压平面上(X=0)切向应力分布云图.(a)2a=29.7:(b)2a=25.7°：(c)2a=14.4°：(d)2a=7.6°：(e)2a=4.4° Fig.9 Contour plots of normalized tangential stresses (on the compressed middle section plane (Y=0)of the 3D disks with different contact loading angles:(a2a=29.7,(b)2a=25.7,(c)2a=14.4°，(d2a=7.6°：(e)2a=4.4° 变准则对三维巴西圆盘起裂点位置进行分析，当 置还受圆盘材料的泊松比影响，较小的泊松比会 圆盘受压直径上某一点的最大张拉应变超过某一 促进中心起裂 限值时就会发生开裂，其本质为最大拉应变强度 表4为9种加载角下三维巴西圆盘试样抗拉 准则.喻勇指出三维巴西圆盘的应力应变分布 强度计算结果，除上述8种加载角外增加2a=33.7° 规律与材料泊松比有关，而与材料弹性模量无关 这一情况，表中Distance表示起裂点至端面圆心的 Li和Wong刀通过数值模拟研究发现，三维圆盘端 距离.当2a=44时，圆盘试样发生端部起裂，峰值 面圆心处的拉应变随泊松比线性增大，而端面受 载荷为5427N,计算抗拉强度最小，而相对误差最 压直径上的最大拉应变随泊松比呈现非线性增大 大为65.5%：随着加载角增大，起裂点逐渐向圆盘 趋势，但是其未探讨受压直径上最大张拉应变位 端面圆心移动，峰值载荷逐渐增大，计算抗拉强度 置的变化规律.图10为泊松比=0.15、0.20、0.25、 逐渐接近试样的真实抗拉强度.但是，当2=29.7° 0.30、0.35时最大切向拉应力点及最大切向拉应变 及33.7时，圆盘试样发生端面中心起裂，此时计算 点位置随接触加载角的变化情况.由图可知，随泊 抗拉强度分别为7.17MPa、7.16MPa,仍低于试样 松比减小，受压直径上最大切向拉应力点及最大 的真实抗拉强度，相对误差均为10.4%.上述计算 切向拉应变点位置均会向端面圆心移动.上述模 及分析结果表明，即便在保证中心起裂的条件下， 拟结果表明，对于三维巴西圆盘，其端面起裂点位 由于三维效应巴西劈裂试验可能仍会低估岩石试

变准则对三维巴西圆盘起裂点位置进行分析，当 圆盘受压直径上某一点的最大张拉应变超过某一 限值时就会发生开裂，其本质为最大拉应变强度 准则. 喻勇[6] 指出三维巴西圆盘的应力应变分布 规律与材料泊松比有关，而与材料弹性模量无关. Li 和 Wong[7] 通过数值模拟研究发现，三维圆盘端 面圆心处的拉应变随泊松比线性增大，而端面受 压直径上的最大拉应变随泊松比呈现非线性增大 趋势，但是其未探讨受压直径上最大张拉应变位 置的变化规律. 图 10 为泊松比 μ=0.15、0.20、0.25、 0.30、0.35 时最大切向拉应力点及最大切向拉应变 点位置随接触加载角的变化情况. 由图可知，随泊 松比减小，受压直径上最大切向拉应力点及最大 切向拉应变点位置均会向端面圆心移动. 上述模 拟结果表明，对于三维巴西圆盘，其端面起裂点位 置还受圆盘材料的泊松比影响，较小的泊松比会 促进中心起裂. 表 4 为 9 种加载角下三维巴西圆盘试样抗拉 强度计算结果，除上述 8 种加载角外增加 2α=33.7° 这一情况，表中 Distance 表示起裂点至端面圆心的 距离. 当 2α=4.4°时，圆盘试样发生端部起裂，峰值 载荷为 5427 N，计算抗拉强度最小，而相对误差最 大为 65.5%；随着加载角增大，起裂点逐渐向圆盘 端面圆心移动，峰值载荷逐渐增大，计算抗拉强度 逐渐接近试样的真实抗拉强度. 但是，当 2α=29.7° 及 33.7°时，圆盘试样发生端面中心起裂，此时计算 抗拉强度分别为 7.17 MPa、7.16 MPa，仍低于试样 的真实抗拉强度，相对误差均为 10.4%. 上述计算 及分析结果表明，即便在保证中心起裂的条件下， 由于三维效应巴西劈裂试验可能仍会低估岩石试 2 −25 3 1 −20 2 0 −15 1 −1 −10 0 −2 −5 −1 0 −2 5 10 15 20 25 −25 −20 −15 −10 −5 0 5 10 15 20 25 Z coordinate of radius/mm Z coordinate of radius/mm Tangential stresses normalized by 2 P/πDt Tangential stresses normalized by 2 P/πDt (a) (b) Z X Z X Y 2α=4.4° 2α=7.6° 2α=10.9° 2α=14.4° 2α=17.9° 2α=21.7° 2α=25.7° 2α=29.7° 2α=4.4° 2α=7.6° 2α=10.9° 2α=14.4° 2α=17.9° 2α=21.7° 2α=25.7° 2α=29.7° 图 8    不同加载角下三维巴西圆盘试样中心剖面及端面受压直径上切向应力分布示意图. （a）中心剖面（Y=0.0125）切向应力；（b）端面（Y=0.025）切 向应力 Fig.8    Normalized tangential stresses along the compressed diameter of the surface and middle section of 3D disks with different contact loading angles: (a) middle section (Y=0.0125); (b) surface (Y=0.025) (a) (b) (e) 0 0.3 0 0.3 0.6 0 0.3 0.6 0 0.3 0.6 0.3 0.6 0.6 0 0.3 0.6 0 0.3 0.6 0.3 0.6 0.6 0.8 0.6 0.8 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 1 1 1 1 (c) (d) 图 9    不同加载角时三维巴西圆盘试样轴向受压平面上（X=0）切向应力分布云图. （a）2α=29.7°；（b）2α=25.7°；（c）2α=14.4°；（d）2α=7.6°；（e）2α=4.4° Fig.9    Contour plots of normalized tangential stresses (σxx) on the compressed middle section plane (X = 0) of the 3D disks with different contact loading angles: (a) 2α = 29.7°; (b) 2α = 25.7°; (c) 2α = 14.4°; (d) 2α = 7.6°; (e) 2α = 4.4° 乔    兰等： 巴西圆盘劈裂二维及三维数值模拟研究 · 139 ·