横观各向同性板岩三轴力学特性及其本构模型

陈晔磊，李地元，王立川 ，李化云，张鹏飞，吴剑

(1. 西华大学 建筑与土木工程学院，四川 成都 610039；2. 中南大学 资源与安全工程学院，湖南 长沙 410083；3. 中南大学 土木工程学院，湖南 长沙 410075；4. 中铁十八局集团有限公司，天津 300222；5. 中铁西南科学研究院有限公司，成都 611731)

岩石的各向异性特征是岩石力学领域的研究重点之一[1]，横观各向同性是各向异性的一种特殊形式，其结构特征、力学参数和应力-应变关系与垂直于各向同性平面的旋转轴对称[2]。板岩具有显著横观各向同性和结构面弱化特征，其破坏模式与层理夹角密切相关[3]。中国西部区域多山，地质条件复杂，广泛分布着大量板岩。随着交通基础设施的大量建设，沿线将不可避免地出现大量穿越板岩地层的隧道等地下工程。因板岩层理发育，层间胶结能力弱，导致工程结构常出现隆起、弯折、断裂等严重病害现象[4]。为此，国内外不少学者对横观各向同性岩体的力学特性[5-6]及其本构模型建立[7]等方面开展了大量研究。岩土体的本构模型对有限元计算有着直接的影响[8-9]，弹塑性本构模型包含了屈服准则、流动法则、硬化定律等理论。Mohr-Coulomb 准则和Drucker-Prager准则在岩土材料的塑性分析中得到了广泛的应用[10-12]，弥补了经典塑性力学仅适用于金属等材料的不足，但这2 种屈服准则主要适用于各向同性材料的分析。在对横观各向同性岩体屈服准则的改进中，层理夹角和软弱层的力学特性是需要被考虑的因素。层理夹角对本构模型的影响不容忽视，与实际地下工程的情况相符。以隧道工程为例，如图1 所示，层理与洞壁切线的夹角随岩石位置的变化而变化(A位置夹角约为45°，B位置夹角约为0°)。在一些本构模型的研究中[12-13]，通常会从微观角度去考虑软弱层的影响，这从一定程度上提高了计算精度，但计算繁琐。当软弱层很小或计算大型开挖模拟时，这种方法是不经济的[14]。横观各向同性材料本构模型参数较多，需多组试验才能确定模型参数[15]，所以参数的确定是值得研究的。

图1 层理夹角与隧道轮廓线的关系Fig. 1 Relationship between bedding angle and tunnel profile

综上所述，本文以板岩为研究对象，通过常规三轴压缩试验的结果，分析其力学特性与层理夹角间的关系。在实验基础上，推导出在整体坐标中的弹性刚度矩阵；采用完全等效模型将层理面效应纳入岩体的连续体描述中，建立适用于板岩的弹塑性本构模型，并求解模型参数。最后通过UMAT 子程序二次开发将其应用于ABAQUS 数值模拟软件中，结合常规三轴试验的结果验证其合理性。

1 板岩三轴压缩试验及力学特性分析

岩石三轴压缩试验主要包括常规三轴压缩试验(σ1>σ2=σ3)和真三轴压缩试验(σ1>σ2>σ3)，本文开展的是圆柱体试样的常规三轴压缩试验。试验所采用的围压分别为5 MPa 和10 MPa(根据板岩取样地所在隧道工程的地应力量级大小确定)，试样层理夹角分别为0°，30°，45°，60°和90°，以此分析围压和层理夹角对板岩力学特性的影响，并为后续本构模型的建立及参数的求解提供可靠的试验数据。

1.1 试验仪器及试样加载

本文采用GCTS RTX-1500 岩石三轴仪，对不同层理夹角的板岩试样(自然风干状态)开展常规三轴压缩力学试验。该试验机轴向最大试验力为1 500 kN，最大围压可达140 MPa。试样取自峨汉高速豹狸岗隧道中的板岩，加工制成不同层理倾角的标准圆柱试样。试样直径为(50±3) mm，高度为100 mm，两端不平整度不超过0.05 mm，端面垂直于试样轴线，最大偏差不超过0.25°。

针对不同层理夹角的板岩试样进行围压分别为5 MPa 和10 MPa 的常规三轴压缩试验。本次试验围压以0.05 MPa/s的速度加载至预设值，并在试验过程中保持不变。轴压以0.01 MPa/s的速度预加载至2 MPa，待试样上方柱状钢座接触压力机上压头后，以0.05 MPa/s的速度加载至试样破坏。试样破坏后如图2所示。

图2 常规三轴压缩试验破坏试样Fig. 2 Failure specimens for conventional triaxial compression tests

1.2 试验结果

图3为不同层理夹角试样在不同围压下的差应力-应变关系。由图3 可知，围压为5 MPa 和10 MPa时，在试样的差应力-应变关系图中均可观察到明显的峰值跌落现象。但在水平层理试样中，当围压为5 MPa时，由于加载过程中试样压碎破坏的响声明显，表现出较强的脆性特征，为了保护试验仪器，防止试样碎屑掉入压力仓，在未观察到明显的应力跌落现象时便停止了加载，但仍可以看出试样达到应力峰值后下降的趋势。本文主要研究岩石试样破坏前的弹塑性变形阶段，所以并不影响最终分析结果。

如图3 所示，板岩的差应力-应变关系受岩石层理夹角和围压的影响。岩石试样的破坏强度随着围压的增加而增大，初始阶段没有观察到差应力-应变曲线上凹，说明无明显微裂隙闭合阶段，主要发生弹性变形。在达到破坏强度之前，差应力-应变曲线下凹，此时岩石试样发生应变硬化，弹性模量减小，产生部分塑性变形。达到峰值强度之后，岩石试样的微裂隙扩张、贯通，呈现应变软化现象。最后，岩石试样进入残余强度阶段，可认为岩石试样已失去承载能力。

图3 不同层理夹角试样在2 种围压下的差应力-应变关系Fig. 3 Differential stress-strain curves of specimens with different bedding angles under two confining pressures

1.3 力学特性分析

基于板岩常规三轴试验的结果，可以获得其在不同层理夹角下的三轴力学特性，为后续本构模型的构建及弹性参数的求解提供依据。

如图4所示，试样峰值强度与层理夹角呈非线性关系。当围压为5 MPa 时，层理夹角为45°的岩石试样的峰值强度最小，其值为112.4 MPa，层理夹角为0°的岩石试样的峰值强度最大，其值为164.6 MPa。当围压为10 MPa 时与围压为5 MPa 时的规律一致，峰值强度最小为45°层理夹角试样，峰值强度最大为0°层理夹角试样。围压从5 MPa到10 MPa 对0°～90°层理夹角的试样峰值强度的提升幅度分别为59%，67%，70%，67%和69%。试样残余强度与围压呈正相关，随着围压的增加，岩石试样的残余强度也随之增大，岩石的塑性增强。试样残余强度与层理夹角呈非线性关系，60°层理夹角试样的残余强度最小。围压对岩石试样径向应变的影响小于对轴向应变的影响以及对应力的影响。如图3(d)所示，当围压由5 MPa 增加到10 MPa 后，岩石试样弹性阶段的轴向应变增大了17%，径向应变增大了6%。这是因为围压的增加，使岩石内部的晶体凝聚力增大，岩石的破坏表现为滑移破坏，围压又从一定程度上限制了其横向滑移，从而增大了岩石的峰值强度，增强了岩石的塑性。

图4 2种围压下板岩试样峰值强度与层理夹角的关系Fig. 4 Relationship between the peak strength and bedding angle of slate specimens under two confining pressures

2 弹塑性本构模型

2.1 弹性本构模型

目前，岩土材料的弹性阶段主要采用在局部坐标系中的广义胡克定律来描述。该定律认为在三维空间中，应力与应变呈线性关系，其表达式为[16]：

式中：E为弹性模量，ν为泊松比，G=E/2(1+ν)，εi(i=x,y, z)为主应变分量，σi(i=x,y, z)为主应力分量。

2.1.1 局部坐标弹性本构模型

如图5所示，当岩体位于局部坐标系中，即不考虑坐标轴旋转时，板岩的弹性本构模型为：

图5 局部坐标系中的板岩试样Fig. 5 Slate samples in local coordinate system

局部坐标中的板岩弹性柔度矩阵[Ce]为：

式中：Ex=Ey，为xy层面内的弹性模量；νxy为xy层面内的泊松比；Ez为垂直于层面的弹性模量；νyz=νxz为垂直于层面的泊松比；Gyz=Gzx为垂直于层面的剪切模量。

2.1.2 整体坐标弹性本构模型

在实际工程中，最小主应力与岩体层理面的夹角会因开挖后的应力重分布而发生改变。如图6所示，假定层理面始终平行于局部坐标系中的xy平面，最小主应力始终平行于整体坐标系中的x＇y＇平面。当局部坐标系绕x轴旋转时，最小主应力与岩体层理面的夹角θ发生改变，该夹角简称为层理夹角。

图6 局部坐标与整体坐标的转换Fig. 6 Conversion of local coordinates to global coordinates

韩昌瑞等[17]基于正交各向异性材料的弹性理论，通过坐标轴旋转，实现了广义胡克定律从局部坐标系到整体坐标系的转换，从而体现了层理夹角因开挖后应力重分布而发生的改变。将广义胡克定律在整体坐标系中表示为：

整体坐标中的板岩弹性柔度矩阵[Ce]'为：

式中：Sij(i,j=1~6)为每个单位应力分量所产生的应变分量。

当层理夹角θ=0 或90°时，[Ce]'可以退化为局部坐标系中的弹性柔度矩阵[Ce]。通过对整体坐标系中的弹性柔度矩阵[Ce]'求逆，可得整体坐标系中的弹性刚度矩阵[De]'：

式中：Kij(i,j=1~6)为每个单位应变分量所对应的应力分量。

2.2 塑性本构模型

岩土材料应力-应变关系的塑性部分通常需要屈服准则、流动法则、硬化定律等理论来描述。但目前常用的屈服准则只适用于各向同性材料，关联流动法则一般适用于金属类材料，不适用于岩土材料。本文考虑层理夹角的影响，对Mohr-Coulomb准则进行了推广，并采用非关联流动法则描述塑性流动，从而推导出合理的板岩塑性本构模型。

2.2.1 屈服准则

以Mohr-Coulomb 定律为基础的屈服准则被广泛采用，但是该准则的屈服面存在尖角，计算时会出现奇点问题，且常用于各向同性材料中，不能反映材料方向性。HILL[16]基于米塞斯准则提出了能反映材料方向性，适用于岩土材料的屈服准则：

式中：F，G，H，L，M，N为模型材料参数，用以描述材料的方向性，其数值通过常规三轴压缩试验的结果求解。

对于Mohr-Coulomb 这类在π平面上为非圆形的屈服函数，可用广义Von·mises 圆曲线去逼近，从而解决了计算时的奇点问题，表达式如下所示[18]：

式中：φ为内摩擦角，c为黏聚力。

在式(8)中的应力张量第二不变量J2引入模型材料参数，使其能反映材料方向性。并且因为板岩的各向异性对称于z轴，模型材料参数F=G=M=L，H=N，则J2可被推广为J＇2：

将式(13)代入式(8)可得：

式中：θσ为应力洛德角μσ为洛德参数，I1为应力张量第一不变量，I1=σx+σy+σz。

2.2.2 流动法则

流动法则是描述塑性流动的规律。Von·mises假设存在某种塑性势函数Q(σ)，其塑性流动方向与塑性势函数Q(σ)的梯度或外法线方向保持一致。基于塑性位势理论可得：

在传统塑性力学中，常采用关联流动法则。该法则将塑性势函数看作屈服函数，即Q(σ)=f(σ)。大量的土工试验表明[19]：相关联的流动法则不考虑应力主轴的旋转，若Mohr-Coulomb 模型采用关联流动法则，将会导致出现远大于实际的剪胀变形。综上所述，关联流动法则不适用于岩土类材料，根据广义塑性理论的规定选取非关联流动法则，Q(σ)≠f(σ)。依据文献[19]中的方法，将Q(σ)假设为：

式中：β为剪胀角，根据文献[20]得β=φ/2，R为塑性拟合参数。

2.2.3 硬化定律

硬化定律是描述应力增量对塑性应变影响的准则，它决定了屈服面的变化情况。本文采用塑性主应变为硬化参量H[19]：

2.3 参数求解

2.3.1 弹性参数

基于围压为5 MPa时的常规三轴试验，可得到Ez，Ex，Ey，νyz，νxz，νxy等6 个基本弹性参数，为本构模型的数值实现提供弹性参数依据。

WANG 等[15]通过水平层理试样求解Ez，νyz和νxz：

式中：dσi为主应力增量，dεi为主应变增量。

对于竖直层理试样，围压不变，dσx=dσy=0，可得如下各弹性参数：

式中：dεr为径向应变，dεr=(dεx+dεy)/2。

基于5 MPa围压作用下水平层理和竖直层理试样的常规三轴试验结果，板岩弹性参数的计算结果如表1所示。

表1 板岩弹性参数Table 1 Elastic parameters of slate

2.3.2 塑性参数

板岩弹塑性本构模型的屈服函数中的F和H为模型材料拟合参数，通过分离差应力-应变关系中的塑性应变增量求得其具体数值，为本构模型的数值实现提供塑性参数依据。

式(15)中的塑性因子dλ可以表示为：

水平层理试样的常规三轴试验中围压保持不变，σx=σy，由势函数得：

径向塑性应变增量：

轴向塑性应变增量：

其中塑性应变增量可以通过试验分离。结合常规三轴试验，联立式(24)和式(25)得到模型材料拟合参数F=4.67，H=1。

3 本构模型的数值实现

本节基于ABAQUS 6.14，对UMAT 子程序进行二次开发。用户材料子程序UMAT 是ABAQUS用于定义特殊材料属性的二次开发接口[21]。

3.1 计算模型与参数

建立标准圆柱试样模型，高100 mm，直径为50 mm。模型再现常规三轴试验全过程。在模拟过程中，采用应力加载的方式。模型单元类型为C3D8I，边界条件采取固定底端竖直方向(Z方向)位移，限制模型顶端和底端的转动，模型的计算参数如表2所示。

表2 ABAQUS模型计算参数Table 2 Calculation parameters of ABAQUS model

3.2 计算结果及模型验证

在围压为5 MPa的条件下，分别对不同层理夹角的试样进行了数值模拟计算，其计算结果与板岩常规三轴试验对比，以验证板岩本构模型的合理性。本文以不同层理夹角试样的轴向应变为例，数值模拟与常规三轴试验的对比结果如图7所示。

图7 模型计算结果验证Fig. 7 Validation of model calculation results

图7表明：板岩弹塑性本构模型的计算结果与常规三轴试验的结果基本吻合，差应力-应变的变化趋势一致，当达到板岩的峰值强度时，应变略微偏小。除90°层理夹角试样以外，其余层理夹角试样均在20～30 MPa 范围内开始出现拐点，这是因为试样开始出现塑性变形逐渐屈服，与试验结果一致。数值模拟结果再次体现了板岩的材料方向性对板岩的变形和破坏强度都产生了一定的影响。除此之外，在数值模拟的过程中发现层理夹角为30°，45°和60°的试样模型计算迭代步数明显多于水平层理和竖直层理，导致数值计算的速度较慢，其原因是当层理夹角为30°，45°和60°时，UMAT 子程序计算循环次数增加，后期研究可继续优化。

4 结论

1) 开展了不同层理夹角、不同围压条件下的板岩常规三轴压缩试验，基于试验结果分析了板岩的力学特性。板岩试样的峰值强度、残余强度都随着围压的增加而增大，但试验范围内的围压对其径向变形的影响有限，板岩峰值强度和残余强度均与层理夹角呈非线性关系，45°层理夹角试样强度最小，0°层理夹角试样强度最大。

2) 建立了板岩弹塑性本构模型，模型弹性部分采用了整体坐标系下的广义胡克定律，并推导出了整体坐标系中的弹性刚度矩阵；采用改进的Mohr-Coulomb 屈服准则和势函数、非关联流动法则和应变硬化准则描述其塑性部分。基于常规三轴试验结果，求解得到了试样的弹、塑性参数。

3) 通过对UMAT 子程序的二次开发，完成了板岩弹塑性本构模型在ABAQUS 中的数值实现，数值计算结果与常规三轴试验结果的对比表明板岩弹塑性本构模型是合理的。