煤岩Holmquist-Johnson-Cook本构模型参数研究

宫能平，张朋朋

(安徽理工大学应用力学研究所，安徽 淮南 232001)

Johnson-Homquist-Cook(HJC)本构模型[1]593是在金属JC本构模型的基础上，考虑了体积应变和静水压的影响，能够很好地描述工程地质材料在高静水压、高应变率和考虑损伤情况下的动态力学性能，在冲击爆破等问题中有着广泛的应用[2]198。数值模拟作为一种常见的研究方法，材料本构模型的选择和参数的确定是研究的重点和难点，将直接影响分析结果的可靠性。所以对材料本构模型的参数进行确定是非常有必要的。文献[3]689通过对HJC模型的分析，确定了抗压强度为8.87MPa、素混凝土HJC模型的参数；文献[4-5]125，145提出了一种通过三轴围压实验获取HJC本构模型极限面参数的方法，分析了混凝土HJC模型参数的敏感性，对工程中混凝土HJC参数的取值给出了建议。文献[6]通过将数值模拟和实验相结合分析比较的方法，得到了C60混凝土HJC模型的相关计算参数，并且分析了HJC本构模型的特性，比较了与金属JC模型的差异。文献[7]14以混凝土的静态抗压强度为基础，利用混凝土计算强度的一些经验公式，同时通过对计算数据的分析，最后给出了一种简易的确定混凝土损伤HJC模型参数的方法。文献[8]将HJC本构模型推广到沥青材料，并通过实验确定了HJC模型所有参数。文献[2]200和文献[9]分别确定了石灰岩和花岗斑岩的极限面参数，并且对所涉及的参数的敏感性进行了分析，给出了一般岩石材料HJC模型参数确定的建议。

但是对煤岩材料HJC模型关键参数的确定却未见报道，在文献[10-11]较多套用已有文献的相似材料的计算参数。本文以淮南潘谢煤田13-1煤层煤岩相关物理、力学实验为基础，通过拟合三轴实验、Hopkinson压杆实验和Hugoniot实验数据，总结出一套获取煤岩HJC模型参数的方法。基于所确定的参数利用有限元分析软件LS-DYNA对煤岩SHPB实验进行了数值模拟验证，并且对关键参数的敏感性进行了分析，对煤岩动载数值计算有参考意义。

1 计算模型及参数分析

1.1 强度模型

HJC强度模型采用特征化的等效应力来描述，表达式为

σ*=[A(1-D)+BP*N](1+Clnε*)

式中：σ*=σ/fc为特征化等效应力，σ、fc分别为实际的等效应力和准静态单轴抗压强度；P*=P/fc为特征化压力，其中P为单元内的静水压力；A、B、N、Smax为极限面参数，在不考虑损伤和率效应时，可通过三轴围压实验获得，C为率效应参数，D为损伤变量，ε*为参考应变率。

1.2 HJC模型损伤的定义

HJC模型的损伤由塑性应变的累积而成，包括体积应变和等效塑性应变，损伤的方程为

1.3 状态方程

HJC模型的状态方程分为弹性阶段、压实变形和压密实后变形三个阶段，用来描述体积应变和静水压之间的关系。其计算方法如下

p=

式中：K为体积模量；u为体积应变；pc和pl分别为煤岩压碎和压实时的压力；uc和ul分别为相应体积应变；K1，K2，K3为压力参数。其状态方程曲线如图1所示。

图1 状态方程曲线

2 煤岩HJC模型参数确定

HJC模型包含极限面参数A、B、N和Smax，损伤参数D1、D2、EFMIIN和率效应系数C，压力参数K1、K2、K3、pc、uc、pl和ul以及基本力学参数fc、G、T和ρ共19个参数。

2.1 基本力学参数的确定

试件取自淮南潘谢煤田13-1主采煤层。在顾桥和潘一井田，埋深400～800m为高挥发烟煤，属气煤阶段[12]。其抗压强度fc通过安徽理工大学采矿工程实验室的RMT岩石力学实验系统测得，如图2所示。

图2 煤样RMT静力加载试验

根据材料力学基本假设和静力学公式可得到煤样基本力学参数如表1所示。

表1 煤岩物理和静力学参数

T为最大拉伸强度，由劈裂实验测得。由于缺少实验，采用文献[7]13中的方法。

T=0.62(fc)1/2

所以T=2.61MPa

特征化拉伸强度T*=T/fc=0.148

2.2 极限面参数组的确定

HJC本构模型的极限面参数包括A、B、N和Smax。文献[4]121根据塑性屈服面理论π平面上Mises准则与Mohr-Coulomb准则的关系推导出了无量纲的特征化粘性强度系数

(1)

K为煤岩三轴实验中线性关系的斜率。本文收集整理了淮南地区-780m煤岩的三轴压缩实验[13]数据，通过数据拟合的方法求出K(见图3)。

围压/MPa图3 煤岩轴压跟围压大小的关系

由图3可知，K值为1.52，代入式(1)可得特征化粘性强度系数A=0.41。

HJC本构模型在不考虑损伤和率效应时，其强度表达式为

σ*=A+BP*N

(2)

在三轴围压实验中

(3)

(4)

图4 拟合得到极限面参数

通过处理较高围压下三轴压缩实验数据，拟合得到了煤岩的极限面参数分别为A=0.41、B=1.95、N=0.76。

2.3 率效应参数确定

文献[14]1通过Hopkinson压杆实验，得到了不同应变率下煤岩的强度。但是煤岩强度的提高不仅包含率效应还包括静水压的影响，所以为了得到特征化煤岩强度σ*与应变率的关系，就必须消除静水压的影响。根据文献[3]690的方法，从特征化的静水压力T*=0.148出发，分别画通过数据点的直线，直线的斜率表示了应变率效应的度量。为了消除静水压影响，过对应于恒定的特征化静水压力p*=1/3处做与纵轴平行的直线，与不同斜率的直线相交，确定不同应变率下煤岩的特征化强度。由直线方程拟合数据点得到C=0.001 47，如图5所示。

(a)消除静水压力的影响

(b)拟合得率效应参数图5 获取参数C值的数据处理方法

2.4 压力参数的确定

由2.3节的介绍，pc和uc分别为煤岩压碎时的压力和体积应变，取值可根据文献[1]591得到

pc=fc/3=5.89MPa，uc=pc/K=0.002 86

参数K1、K2、K3、pl和ul可根据冲击压缩经验常数的Hugoniot关系进行拟合选取。冲击压缩经验常数表示的Hugoniot关系式为

p=C2ρ0u(u+1)/[(1-S)u+1]2

(5)

式中：C和S为经验常数，可由文献[15]获得。C=917.4m/s，S=1.446。对比状态方程第三阶段，通过拟合可得K1=1.6GPa，K2=-0.17GPa，K3=5.8GPa。拟合结果如图6所示。同样可获得pl=0.34GPa。

图6 拟合得压力参数

2.5 其余参数的确定

文献[7]13给出了参数D1的计算方法

D1=0.01/(1/6+T*)

其中，T*=T/fc=0.148，所以D1=0.032

其余参数参照文献[16]， 最终确定的煤岩HJC模型参数包括基本力学参数ρ=1 410kg/m3，fc=17.67MPa，G=1.36GPa，T=2.61MPa；极限面参数A=0.41，B=1.95，N=0.76，Smax=7.7；压力参数pc=5.98MPa，uc=0.002 8，K1=1.6GPa，K2=-0.17GPa，K3=5.8GPa，pl=0.34GPa，u1=0.1；损伤参数D1=0.032，D2=1，EFMIN=0.01和率效应参数C=0.001 47。

3 参数验证及敏感性分析

3.1 数值模拟验证

SHPB实验在安徽理工大学能源与安全学院Φ75霍普金森压杆实验系统上完成，钢质压杆弹性模量为210GPa，密度为7 890kg/m3、泊松比为0.3，输入杆和输出杆长度分别为4 000mm和2 500mm，子弹长度为400mm，试件的长度和直径分别为50mm和75mm。按照实际尺寸建模，依据所确定的煤岩HJC模型参数，对SHPB实验进行数值模拟研究。重构了煤岩的动态应力应变曲线，与实验结果[14]3进行对比如图7所示。

图7 实验和数值模拟的动态应力-应变曲线结果对比

从图7可知，数值模拟结果与实验的结果吻合较好，说明本文所采用的参数确定方法是有效的。煤岩的动态应力应变曲线表现出明显的应变硬化现象，且与无烟煤的实验结果[17]对比得知，淮南13-1煤层煤岩破坏的应变更小，表现出高更明显的脆性，与文献[18]的结论一致。

3.2 参数敏感性分析

在已有的HJC本构模型参数确定与敏感性分析中，都发现极限面参数A、B对结果的影响较大[2，5]197,147，本文保持确定的参数不变，仅通过改变目标参数的方法来研究目标参数对煤岩动态力学性能的影响。

1) 参数A的影响

参数A为归一化的粘性强度，其定义为在一定的静水压力下，无损伤强度(D=0)与完全损伤(D=1)之差。参数A的取值范围约为0.35～0.90。在数值模拟中，保持所确定的参数不变，仅通过改变A的值，可得到此范围内不同A值下的应力-应变曲线，如图8所示。

图8 参数A对煤岩动态应力-应变曲线的影响

从图8可以看出，当A的取值从0.35增加到0.55时，峰值应力从46.66MPa增加到50.33MPa，增加了7.87%，假设在此范围内应力峰值呈线性增加，则A每增加0.1，峰值应力增加3.94%；同样，在A取值在0.55到0.75时，峰值应力从50.33MPa增加到52.0MPa，A每增加0.1，峰值应力增加1.66%；在A取值在0.75到0.9范围内时，峰值应力从52.0MPa增加到52.33MPa，峰值应力几乎没有变化。说明在0.35～0.9范围内，A的增加能够显著提高峰值应力，但是增长率降低。此外，峰值应变也会随着A的增加而变大，变化量较小。

2) 参数B的影响

参数B为归一化压力硬化系数，B值越大，说明静水压力在屈服方程中所占的比例越大[19]。参数B的取值范围约为1.4～2.0。在数值模拟中，保持所确定的参数不变，仅通过改变B的值，可得到此范围内不同B值下的应力-应变曲线，如图9所示。

图9 参数B对煤岩动态应力-应变曲线的影响

从图9可以看出，参数B从1.4增加到1.6时，峰值应力从49.57MPa增加到50.62MPa，增加了2.12%；参数B从1.6增加到1.8时，峰值应力从50.62MPa增加到51.41MPa，增加了1.56%；参数B从1.8增加到2.0时，峰值应力从51.41MPa增加到52.33MPa，增加了1.79%。说明在1.4～2.0范围内，参数B可以提高峰值应力，增长率较均匀。参数B对峰值应变的影响较小。

4 结论

(1)将HJC本构模型推广到煤岩材料，确定了模型参数，给出了煤岩HJC模型各参数的确定方法。基于确定的参数的数值模拟结果与SHPB实验结果吻合的较好，说明本文的参数确定方法是有效的。可为煤岩动载数值计算提供帮助。

(2)对极限面参数A和B对煤岩动态性能的影响规律进行了研究，结果表明，在一定范围内，随着A的增大，峰值应力和峰值应变都有一定程度的增长，但峰值应力的增长率降低；参数B的增加能够提高峰值应力和峰值应变，且B对峰值应力的影响幅度要大于对峰值应变的影响。