煤层本构关系与数值模拟实战：从选型到参数标定

1. 本构关系到底是啥——先讲清楚这个“工程地基”

搞采矿、做岩石力学数值模拟的人，几乎都绕不开“本构关系”这个词。但说实话，不少同行对这个概念的理解停留在“应力应变曲线拟合”这个层面，要么从教材上抄一段Drucker-Prager参数应付评审，要么直接从文献里扒一套Mohr-Coulomb参数用到底。这么做其实挺危险的，因为煤层本构关系是整个数值模型的“地基”，地基歪了，上面盖的楼再漂亮也白搭。

所谓本构关系，通俗讲就是材料“受力之后怎么变形、怎么破坏”的数学描述。你给煤样加压，它先压缩、再弹性回弹、然后塑性流动、最后峰后软化直至残余强度，这一整条应力应变曲线背后的规律，就是本构关系要刻画的内容。数值模拟里，网格节点上的位移、应变、应力全靠本构方程把“外力”和“变形”关联起来。可以说，本构关系选得对不对，直接决定了你模拟出的巷道变形量是3毫米还是3米，决定了你判断的冲击危险等级是“无”还是“强”。

拿煤炭开采来说，我们关心的问题非常多：采掘工作面的围岩稳定性怎么评价？煤柱尺寸留多大才安全？冲击地压什么时候可能发生，能不能提前预警？瓦斯抽采钻孔的塌孔风险怎么预测？这些问题听起来五花八门，但落到数值计算层面，统统都要靠本构关系来回答。所以这篇文章我打算从工程应用视角出发，把煤层本构关系的来龙去脉、模型选型、参数标定、数值实现到工程落地的完整链路都梳理一遍。无论是正在写论文的研究生，还是做采场设计、冲击地压防治的工程师，这篇文章应该都能给你一些用得上的参考。

需要说明的是，本文中很多具体参数和操作细节来自于我个人的项目实践和行业通用做法，不同矿区煤岩性质差异很大，绝不能生搬硬套，但思路和流程是通用的，可以照着这个框架去搭建自己的方案。

2. 为什么煤的本构关系不好搞——从煤的“性格”说起

2.1 煤不是普通岩石，它是“浑身毛病”的复合材料

我在做第一个煤矿数值模型时，天真的想法是用经典的弹塑性模型把煤层一包，参数拿实验室单轴抗压强度反算就完事了。结果算出来的采动应力分布跟现场微震监测数据差得离谱，后来跟导师聊完才明白，煤这种材料跟砂岩、灰岩完全不是一个量级的东西。

煤的“性格”极其复杂。首先它是非均质的，煤体里面散布着镜煤、亮煤、暗煤、丝炭这些不同显微组分，成分差异导致力学性质在毫米尺度上就变化很大。其次煤是强各向异性的，因为煤层是沉积形成的，层理面、节理裂隙非常发育，平行层理方向的抗压强度往往只有垂直层理方向的一半甚至更低。再加上煤体内天然裂隙（割理系统）纵横交错，这些弱面在受力过程中会率先张开、滑移、贯通，最终形成宏观破坏面。

更麻烦的是，煤层还经常和夹矸层互层，一整套“煤-矸-煤”组合体的变形破坏规律远超单一岩性的描述范畴。所以做数值模拟时，如果我们拿一个均匀的、各向同性的线弹性模型去描述煤层，那基本等于用“一个标准身材的模特”去代表“所有体型的人”，误差大是必然的。

2.2 峰后软化——煤最容易被忽略的“致命伤”

我接触的很多数值模型里，大家用理想弹塑性模型（比如Mohr-Coulomb）模拟煤层，峰值强度一过，应力就保持不变。这样做对“弹性区判定”或许凑合够用，但一旦涉及巷道大变形、冲击地压启动、煤柱失稳这些问题，理想弹塑性就完全罩不住了。

煤层是典型的峰后应变软化材料。什么意思呢？就是煤样在达到峰值抗压强度之后，承载能力并不会维持不变，而是随着应变继续增加快速下降，降到一定程度后保持一个较低的残余强度。这个峰后软化段的斜率（软化模量）和残余强度的取值，对模拟结果的影响极其敏感。我做过的采动应力演化模拟里，软化参数调个20%，采场塑性区宽度能从15米跳到25米，这个差异足以改变支护方案和煤柱尺寸设计。

深层原因在于煤的破坏本质是内部微裂隙的萌生、扩展和贯通。峰后阶段微裂隙大量发展，试件的有效承载面积急剧减少，宏观表现就是“强度垮塌式衰减”。如果本构模型不考虑这个过程，等于假设煤体破坏后还能继续扛住应力，这与现场巷道两帮大面积片帮、煤炮频发的实际现象完全矛盾。

2.3 更多让人挠头的特殊行为

除了峰后软化，煤层本构关系还得考虑几件“糟心事”。第一是应变率效应，同样是煤，静态加载和冲击加载下的强度可以差30%到50%甚至更多。冲击地压是动力现象，加载速率极高，这时候拿静态参数模拟动态过程，结果可靠性就打了折扣。第二是围压效应，深部煤层处于三向应力状态，围压对煤的强度、塑性变形能力和破坏模式影响显著，本构参数必须随围压水平做相应调整。第三是瓦斯和水的耦合影响，瓦斯吸附会降低煤的有效应力，水的存在会弱化煤的强度，这些因素在某些场景下（比如高瓦斯矿井的卸压抽采）不能忽略。

说到这儿你应该明白了，煤层的“本构关系”不是一个简单的数学模型，而是要把煤这种天然缺陷材料的非线性、非均质、各向异性、软化、率相关、环境敏感等行为尽可能多地包容进来的一套系统性描述体系。

3. 常用本构模型怎么选——从实战场景出发的参数对照

3.1 各模型适用场景速查

做数值模拟这些年，我接触过的煤层本构模型不下十种，但日常工程应用中真正扛大梁的其实就那么几类。我直接按应用场景给你做个对照表，方便快速选型：

3.2 为什么我偏爱“应变软化+M-C”组合

对不同工程问题，我自己的选型习惯是这样：如果只是想快速摸一下采场应力分布、塑性区范围，用带抗拉截断的Mohr-Coulomb就够了，它简单稳定，不容易翻车。但如果涉及煤柱留设尺寸论证、冲击危险区域圈定，我强烈建议上应变软化模型，具体做法是让粘聚力和内摩擦角在峰后随塑性应变线性折减，直到残余值。这种改良方案保留了Mohr-Coulomb的简洁框架，同时抓住了煤体峰后软化破坏的核心机制，计算代价不大，工程可靠性却高了一个档次。

有人可能会问，为什么不直接上损伤模型或者离散元？我的经验是，损伤模型的理论虽然漂亮，但损伤变量的演化方程怎么定、怎么标定，学界到现在也没有统一标准，工程应用容易“杀鸡用牛刀”。PFC这类离散元软件做煤岩破坏机理研究确实厉害，但参数需要靠“试错标定”来匹配宏观力学响应，一个双轴压缩实验就要调半天，想做大型采场模型，计算量也吃不消。所以工程模拟我一般优先考虑连续介质框架下的“应变软化弹塑性模型”，性价比最高。

3.3 Drucker-Prager其实很少单独用在煤层上

这里提醒大家一个误区。很多做三维数值模拟的朋友习惯把Drucker-Prager（D-P）模型当万能弹塑性模型用，因为它在ABAQUS里内嵌得不错，而且屈服面光滑，数值收敛性好。但煤层本身就是层理面控制的剪切破坏和拉伸破坏并存，D-P模型的屈服面在π平面是圆形，无法反映煤体拉压强度不等和各向异性的特性，在某些应力路径下给出的破坏模式会和实际差很多。所以我的建议是：能用M-C就不轻易用D-P；用D-P时一定要通过参数换算公式，让D-P屈服面与M-C屈服面在主应力空间的关键点（如单轴拉、单轴压）保持一致，否则计算结果会和设计规范严重脱节。

4. 参数标定与获取——数值模拟的“良心工程”

4.1 室内实验与数据处理全流程

参数标定是整个本构模拟里最繁琐、也最考验功夫的环节。数值模拟圈流传一句话：“参数不准，算出来的就是高级垃圾。”这话糙理不糙。

标准流程第一步是搞室内实验。最基础的是单轴压缩实验，能得到弹性模量E、泊松比ν、单轴抗压强度σ_c；三轴压缩实验在不同围压（比如2MPa、4MPa、6MPa、8MPa）下得到一组峰值强度和对应的塑性变形数据，用于标定粘聚力c和内摩擦角φ；直接拉伸或巴西劈裂实验用于标定抗拉强度。如果做应变软化模型，还得多做几个循环加卸载实验，用来观测峰后承载力退化规律，从而确定软化参数。

数据处理这儿有几个细节经验。做三轴实验的时候，试样两端一定要打磨平，不然端部效应带来的假强度数据会直接污染后续标定；加载速率控制在0.05mm/min到0.1mm/min之间比较合适，太快的速率下煤样容易“假脆性”，强度偏高；煤样含水率要记录清楚，因为饱和煤的强度比干燥煤可以低30%，这个不记录，后面做对比分析就说不清了。

4.2 三轴数据怎么换算成模型参数——详细步骤

用Mohr-Coulomb模型举例。你手头有一组不同围压下的峰值强度数据(σ₃, σ₁)，处理方法如下：

第一步，绘制摩尔圆。每个围压σ₃对应一个主应力差（σ₁-σ₃），画出一组摩尔圆。

第二步，做摩尔包络线。把所有摩尔圆的公切线画出来（一般采用最小二乘法拟合线性包络线）。

第三步，计算粘聚力和内摩擦角。包络线在τ轴上的截距就是粘聚力c，包络线与σ轴的夹角就是内摩擦角φ。具体公式是：

τ = c + σ·tanφ

其中σ是正应力，τ是剪应力。如果用主应力参数表示，峰值强度线可以写成：

σ₁ = σ_c + k·σ₃

这里的σ_c是单轴抗压强度，k是围压影响系数，它们和c、φ的换算关系为：

σ_c = 2c·cosφ / (1 - sinφ)

k = (1 + sinφ) / (1 - sinφ)

我项目里某矿煤样的三轴实验数据大概是这样：单轴抗压强度12.8MPa，围压4MPa时峰值强度38.6MPa，围压8MPa时峰值强度52.3MPa。用上面公式反算下来，c大约是3.6MPa，φ约29°，弹性模量E约2.4GPa，泊松比ν约0.28。这几个数值供你参考量级，但每个矿的煤质不一样，必须实测。

4.3 峰后软化参数的工程标定与反演

应变软化模型的峰后参数标定相对复杂，业界也没有完全统一的标准。我的做法是这样的：先通过三轴实验的峰后段曲线确定软化段斜率，换算成塑性应变软化模量；然后用FLAC3D或ABAQUS做几个不同软化参数组合的单轴压缩数值实验，和室内单轴实验曲线对比，逐步逼近。实际操作中，软化段的粘聚力折减系数我通常取0.1到0.3，内摩擦角折减幅度则小一些，通常不折减或者从峰值的29°降到残余的26°左右。

这里有个实用技巧：如果你连三轴实验数据都凑不齐，可以参考《煤与岩石物理力学性质测定方法》这类规程里的经验关系，用纵波速度VP估算弹性模量E，公式是E = ρ·Vp²·(1+ν)(1-2ν)/(1-ν)。另一个经验公式是依据单轴抗压强度估算粘聚力，c = σ_c·(1 - sinφ)/(2cosφ)，但这个公式要求你先估一个合理的φ值，一般取25°~35°之间。用估算参数做出来的模型只能用于方案预研，不建议直接用于工程设计。

4.4 参数标定的几个大坑

• 坑一：只用单轴压缩数据标定一切。单轴实验只反映一种应力路径，无法得到围压效应信息，做深部采场模拟必然失真。
• 坑二：忽视尺寸效应。实验室50mm×100mm的煤样强度和现场煤体强度差好几倍，通常现场煤体强度是实验室强度的0.3到0.7倍，需要做尺寸修正。
• 坑三：参数单位搞混。用MPa还是Pa，用米还是毫米，一个疏忽全盘皆输。我的习惯是在建模前把所有单位统一写在一张纸上，挂在屏幕上。

5. 数值实现与仿真实操——从模型搭建到收敛调参

5.1 主流软件怎么选

煤层本构关系最终要落到数值模拟软件里才有工程价值。目前矿业圈主流软件有FLAC3D、ABAQUS、RFPA和PFC。我的使用感受是：FLAC3D在采动应力演化、巷道围岩稳定性这类岩土工程问题上是老牌王者，内置的应变软化模型和FISH语言扩展都很成熟，上手也相对容易；ABAQUS的优势是强大的非线性求解能力和二次开发接口（UMAT/VUMAT），做自定义本构研究的时候优势明显，但三维采矿模型建模和地应力平衡比FLAC3D麻烦；RFPA基于统计损伤理论，能直观模拟煤岩破裂萌生到贯通的全过程，很适合做破坏机理展示和声发射特征研究；PFC是离散元代表，适合研究裂隙扩展机制这类细观问题，但工程尺度应用效率低。

5.2 FLAC3D配置应变软化模型的完整步骤

我用FLAC3D比较多，给新手一个完整的参数配置流程参考。命令大致是这样的：

; 定义煤层材料参数 model mohr-coulomb range group 'coal' property density 1400 bulk 1.67e9 shear 0.94e9 range group 'coal' property cohesion 3.6e6 friction 29 tension 0.8e6 range group 'coal' ; 设置应变软化特性 property table-cohesion 'coal_coh_table' table-friction 'coal_fri_table' range group 'coal'

上面的table需要创建两个表格，分别定义粘聚力和内摩擦角随塑性剪应变的折减关系：

; 定义粘聚力软化表：塑性应变0对应3.6MPa，0.02对应0.5MPa table 'coal_coh_table' 0 3.6e6 0.005 2.4e6 0.01 1.2e6 0.02 0.5e6 ; 定义内摩擦角软化表：塑性应变0对应29度，0.02对应26度 table 'coal_fri_table' 0 29 0.005 28 0.01 27.5 0.02 26

这里要特别说明，软化表第一列是塑性剪应变（无量纲），不是总应变。总应变包含弹性部分，如果拿总应变来控制折减，弹性阶段就已经在掉强度了，模型会“一加载就破坏”，物理上是错的。

跑稳态计算时，FLAC3D的默认求解器基于显式时步迭代，平衡判据是最大不平衡力与典型内力之比小于1e-5。如果模型大、网格多、软化严重，迭代容易震荡。我处理这类问题有几个习惯：一是给煤层网格加密但不过度，软化带区域（比如巷道周边2倍巷径范围内）网格尺寸控制在0.5米左右；二是采用大变形模式（set large on）计算巷道大变形问题；三是如果计算发散，先把软化折减速度放慢（即加大折减对应的应变区间），等模型稳定后再逐步调到目标值。

5.3 ABAQUS用户子程序实现自定义本构的框架

如果要写自定义本构，比如考虑损伤或者应变率效应的煤体本构，那就绕不开ABAQUS用户材料子程序UMAT。UMAT的编程框架我简单说一下：每个增量步开始，ABAQUS会传入当前应变增量Δε和状态变量；你要做的核心工作分三步——根据当前应力状态判断是否屈服、计算塑性流动方向和塑性乘子、更新应力和刚度矩阵DDSDDE。

UMAT的大致骨架是这样：

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT, 1 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS, 2 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT, 3 LAYER,KSPT,KSTEP,KINC) INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS) DIMENSION DSTRAN(NTENS),STRAN(NTENS),TIME(2),PROPS(NPROPS) DIMENSION COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) C PROPS(1)=E, PROPS(2)=NU, PROPS(3)=C0, PROPS(4)=PHI0 C PROPS(5)=C_RES, PROPS(6)=PHI_RES, PROPS(7)=H (软化模量) ... RETURN END

新手写UMAT最头疼的是更新一致切线刚度矩阵DDSDDE，很多教程直接给一个弹性刚度矩阵应付了事，这样做的后果是计算收敛速度慢甚至发散。我的经验是用“数值扰动法”求切线刚度，就是给应力增量加一个微小扰动（比如1e-8），重新计算应力更新，用差分近似偏导数，虽然多耗一点计算量，但通用性和稳定性高很多，非常适合做工程项目的朋友。

5.4 数值模拟的三大常见病

第一个常见病是网格敏感性。软化模型存在应变局部化问题，塑性区宽度可能严重依赖网格尺寸。解决思路有两个：一是给模型引入特征长度（比如Bazant提出的裂纹带模型思想），把软化参数改为网格尺寸的函数；二是用Cosserat连续介质理论改造本构，工程上用得少但理论意义大。实用层面我推荐第一种，用“等效塑性应变乘以网格等效边长”来修正软化曲线，效果立竿见影。

第二个常见病是地应力平衡做不好。采场模型第一步必须把初始地应力场平衡出来，否则后续开挖模拟全是病态结果。FLAC3D的常用做法是先solve elastic得到初始应力场，然后set displacement为零，再切到塑性本构继续计算。这里一个小细节：煤层和岩层弹性模量差异大，地应力平衡时容易出现卸载反弹，建议在模型四周采用应力边界，顶部通过施加重力应力条件来控制。

第三个常见病是塑性流动法则的选取。M-C模型在FLAC3D里的默认流动法则可能不是关联流动，但如果不考虑剪胀角，煤体破坏后的体积膨胀效应会被严重低估。现场煤炮发生时巷道底鼓、两帮剧烈扩容，这部分变形全靠剪胀角控制。我建议剪胀角取内摩擦角的1/4到1/2，具体值用三轴压缩的体积应变曲线来标定。

6. 典型工程应用与常见问题速查

6.1 四种典型的工程应用场景

本构关系选对、参数准了，能解决哪些实际问题？我挑四个常见的应用场景讲。

第一个场景是巷道围岩稳定性分析与支护设计。用应变软化模型模拟煤巷开挖后围岩的塑性区范围和变形量，和现场用多点位移计测得的实际位移做对比验证，可以反过来优化锚杆锚索参数。比如某矿回采巷道用软化模型算出顶板下沉量是26毫米，而用理想弹塑性模型算出来只有9毫米，现场实测数据大约是22毫米，差距一目了然，如果不做软化修正，支护设计明显偏保守。

第二个场景是煤柱设计。煤柱尺寸是压覆资源量和安全之间的平衡点，用带软化模型算出来的煤柱屈服区宽度和应力分布，比解析公式（比如Wilson公式）更贴近复杂地质条件，尤其是断层附近的构造应力区。做这种事前方案比选，数值方法的优势就是可以批量改变参数，形成响应面，供决策用。

第三个场景是冲击地压危险性评价。煤体峰后软化行为与冲击倾向性密切相关，通过数值模拟孕灾过程（应力集中→塑性应变累积→软化加剧→失稳突变），可以划分冲击危险区域并指导卸压钻孔、深孔爆破的布置。这个场景我项目里用过，微震定位数据和模拟的塑性应变高值区拟合度在70%以上。

第四个场景是瓦斯抽采与渗透率演化分析。把本构模型计算的应变场耦合到煤体渗透率模型（比如经典的立方体模型或指数损伤模型）中，模拟抽采过程中煤体变形-渗透率的动态耦合过程，对高位钻孔优化布置有直接指导意义。

6.2 高发问题排查速查表

6.3 本构参数与模型的匹配性检查清单

最后再给大家一个自查清单，建模前逐项打钩，能帮你省掉大量返工时间：

• 是否明确了模拟对象处在地下多少米，对应的原岩应力场大小？
• 煤层是否分层处理？夹矸层、顶底板岩层是否按各自材料参数建模？
• 煤样的实验室参数是否做了尺寸修正、含水率修正？
• 是否根据工程问题选择了合适的本构模型（不要一刀切全用线弹性或理想弹塑性）？
• 软化模型是否定义了粘聚力/内摩擦角随塑性应变的折减表？
• 剪胀角是否结合体积应变曲线做了合理估计？
• 是否做了网格敏感性验证（至少两套不同密度网格对比）？
• 模拟结果是否与现场实测（位移、应力计、微震）做了对比校验？

我自己的经验是，一个真实项目的数值模拟，前期的参数标定和模型验证往往要占用60%的时间，真正跑计算反而很快。很多论文和报告里轻描淡写的一句“参数依据室内实验获得”，背后其实是大量枯燥且关键的标定工作。

最后说一点个人体会。做煤层本构关系这条路，入门容易入深难，难就难在煤这种东西太“不乖”——它既不是标准的弹性体，也不是理想的弹塑性体，而是集软化、损伤、流变、率相关、环境敏感于一身的复杂材料。但恰恰因为复杂，这项工作才有价值。我踩过的最大坑就是迷信某一个“先进模型”而忽视参数标定这个基本功，模型的架子搭得再花哨，参数靠拍脑袋定，结果就是自欺欺人。踏踏实实做实验、认认真真做标定、老老实实做验证，这句话送给每一个正在做数值模拟的同路人。