工程科学学报，第44卷，第X期 多边形之间形成的接触定义为沿晶接触；(3)通过 表2矿物品粒物理、力学参数 连接多边形颗粒的中心和角点对多边形进行进一 Table 2 Physical and mechanical parameters of grains 步离散，内部新产生的接触定义为晶内接触，为穿 Grain type Percentage/%plg.cm) S/GPa B/GPa 晶破坏提供了可能的路径.通过赋予沿晶接触和 Q 27 2.650 44.0 37.0 0.08 品内接触不同的力学性质就可以实现岩石材料的 B 5 2,850 12.4 41.1 0.36 沿晶及穿晶破坏的模拟.图3给出了含有4类矿 K 58 2,560 27.2 53.7 0.28 物成分的数值试件（未考虑原生裂隙和孔洞；其 10 2.630 29.3 50.8 0.26 中，①代表钾长石穿晶开裂，②代表黑云母沿晶 Note:Q,B,K and P represent quartz,biotite,K-feldspar and plagioclase 开裂))，并给出了穿晶、沿晶模拟结果，其微观开 respectively;o represents Poisson's ratio;BS and p represent bulk modulus,shear plane modulus and density respectively 裂结果和真实开裂一致，矿物具体成分和参数见 表2.同时，选择库伦滑移模型作为接触的本构模 线弹性阶段的法向和切向应力-位移之间的关 型.当接触上的剪力或者张力超过对应的强度极 系为刃： 限时，接触会产生相应的剪切或张拉破坏，块体采 △cn=-kn△ln (1) 用弹性模型，图4详细描述了接触本构关系，其在 △Ts=-ks△G (2) 100mm Quartz Trans-granular cracking lagioclas -feldspar Intergranular cracking Biotite 图3模型配置 Fig.3 Model configuration 触的剪应力达到其对应抗剪强度tmax=cp+tanp 时，接触剪应力修正为：T=c+tang;c和cp分别代 ar failure 表接触的残余粘聚力及其峰值；和分别代表残 Tension failure 余摩擦角及其峰值；tmax和t分别代表峰值剪应力 和残余剪应力. 此外，离散单元法中的参数和物理世界的参数 有一定差距，需要通过数值试验反复校核，以达到 模拟真实岩石宏观响应（弹性模量、泊松比、单轴或 Constant o 多轴强度及起裂应力等)的目的.本文参考前面描 u.(u) 述的花岗岩试件的单轴试验结果进行大量的微观 Overlap 参数校核，其校核后的参数和结果见表2~5，数值 和物理单轴应力-应变曲线见图5.可见，数值得到 图4本构关系 的宏观力学、强度参数同物理试验结果基本一致 Fig.4 Constitutive relationship 就应力-应变曲线而言，数值和物理试验的结果具 式中，kn、k、△cn、△rs、△和△分别为法向刚度、 有一定的差异，前者并没有加载前期的压缩阶段， 切向刚度、法向应力增量、切向应力增量、法向位 直接呈现线弹性特征，这已在众多研究成果中被发 移增量和切向位移增量.当接触的法向应力σ超 现0-刘以上差异是由于本文未考虑岩石的原生裂 过其抗拉极限T时，法向应力σ，则减小为0：当接 隙、孔洞，且块体被赋予线弹性本构.应力-应变曲多边形之间形成的接触定义为沿晶接触；(3) 通过 连接多边形颗粒的中心和角点对多边形进行进一 步离散，内部新产生的接触定义为晶内接触，为穿 晶破坏提供了可能的路径. 通过赋予沿晶接触和 晶内接触不同的力学性质就可以实现岩石材料的 沿晶及穿晶破坏的模拟. 图 3 给出了含有 4 类矿 物成分的数值试件（未考虑原生裂隙和孔洞；其 中，①代表钾长石穿晶开裂，②代表黑云母沿晶 开裂）），并给出了穿晶、沿晶模拟结果，其微观开 裂结果和真实开裂一致，矿物具体成分和参数见 表 2. 同时，选择库伦滑移模型作为接触的本构模 型. 当接触上的剪力或者张力超过对应的强度极 限时，接触会产生相应的剪切或张拉破坏，块体采 用弹性模型，图 4 详细描述了接触本构关系，其在 线弹性阶段的法向和切向应力–位移之间的关 系为[17] ： ∆σn = −kn∆un （1） ∆τs = −ks∆u e s （2） mm022 200 mm 100 mm Quartz Plagioclase Biotite ① Intergranular cracking Trans-granular cracking K-feldspar ② 图 3    模型配置 Fig.3    Model configuration Tension failure Shear failure 1 T 1 Constant σn Overlap Compression τmax τr ks kn kn ks un (us ) 图 4    本构关系 Fig.4    Constitutive relationship kn ks ∆σn ∆τs ∆un ∆u e s σn σn 式中， 、 、 、 、 和 分别为法向刚度、 切向刚度、法向应力增量、切向应力增量、法向位 移增量和切向位移增量. 当接触的法向应力 超 过其抗拉极限 T 时，法向应力 则减小为 0；当接 τmax = cp +tanφp τr = cr +tanφr cr cp φr φp τmax τr 触的剪应力达到其对应抗剪强度 时，接触剪应力修正为： ； 和 分别代 表接触的残余粘聚力及其峰值； 和 分别代表残 余摩擦角及其峰值； 和 分别代表峰值剪应力 和残余剪应力. 此外，离散单元法中的参数和物理世界的参数 有一定差距，需要通过数值试验反复校核，以达到 模拟真实岩石宏观响应（弹性模量、泊松比、单轴或 多轴强度及起裂应力等）的目的. 本文参考前面描 述的花岗岩试件的单轴试验结果进行大量的微观 参数校核，其校核后的参数和结果见表 2～5，数值 和物理单轴应力–应变曲线见图 5. 可见，数值得到 的宏观力学、强度参数同物理试验结果基本一致. 就应力–应变曲线而言，数值和物理试验的结果具 有一定的差异，前者并没有加载前期的压缩阶段， 直接呈现线弹性特征，这已在众多研究成果中被发 现[20–21] . 以上差异是由于本文未考虑岩石的原生裂 隙、孔洞，且块体被赋予线弹性本构. 应力–应变曲 表 2    矿物晶粒物理、力学参数 Table 2    Physical and mechanical parameters of grains Grain type Percentage/% ρ/(g·cm–3) Sh /GPa Bu /GPa υ Q 27 2,650 44.0 37.0 0.08 B 5 2,850 12.4 41.1 0.36 K 58 2,560 27.2 53.7 0.28 P 10 2,630 29.3 50.8 0.26 Note: Q, B, K and P represent quartz, biotite, K-feldspar and plagioclase respectively; υ represents Poisson's ratio; Bu , Sh and ρ represent bulk modulus, shear plane modulus and density respectively. · 4 · 工程科学学报，第 44 卷，第 X 期