基于β-屈服函数的炭质泥页岩统一硬化模型

2023-02-24崔凯胡斌崔阿能马利遥刘杨祝鑫

土木与环境工程学报 2023年1期

崔凯,胡斌,崔阿能,马利遥,刘杨,祝鑫

（武汉科技大学资源与环境工程学院；冶金矿产资源高效利用与造块湖北省重点实验室，武汉 430081）

在中国中南、西南地区，如云南、四川、贵州、广西等省存在大范围的含缓倾软弱夹层（岩层和夹层倾角小于25°）的二叠系石灰岩地层[1-3]，该地层拥有丰富的优质石灰石矿产资源，是中国大量基础设施建设的重要建材来源。如图1所示，这些含缓倾软弱夹层通常由炭质泥页岩组成，具有抗剪强度小且容易在外部因素(人为扰动或环境)作用下发生劣化等特点，因此，含缓倾软弱夹层被视作矿山高边坡的薄弱带[4-5]。深入探究炭质泥页岩的力学特性并建立相应的本构关系，对于含缓倾软弱夹层的矿山高边坡的稳定性分析具有重要意义。

图1 矿山高边坡中含炭质泥页岩的软弱夹层Fig.1 Weak intercalation of carbonaceous shale in the high slope of the mine

在复杂加载条件下，炭质泥页岩类似于超固结黏土，会发生弹塑性破坏，其应力—应变曲线具有应变软化的特点，并且伴随着剪胀现象的发生[6-7]。为此，有学者[8]建议将炭质泥页岩等强度较小的软岩视作超固结比很大的黏土，并采用临界状态土力学的理论框架来描述其常规力学特性。宋丽等[9]建立了软岩的三维统一弹黏塑性本构模型，该模型对于日本石川县能登半岛端部的硅藻质泥岩在不同围压和加载速率等条件下的不排水三轴压缩试验结果具有较好的描述效果。熊勇林等[10]同时考虑了围压和温度对软岩残余强度的影响，建立了软岩的热弹黏塑性模型。炭质泥页岩另一个显著的力学特性是具有流变性，在保持荷载强度不变的情况下，增加加载时间将导致额外的体积变形[11]。范庆忠等[12]提出了软岩的非线性蠕变模型，齐明山[13]则系统地研究了大变形软岩流变性态及其在隧道工程结构中的应用。在所有临界状态本构关系中，修正剑桥模型[14]在岩土工程中应用最广泛，然而其只适用于正常固结黏土力学特性的描述，在描述炭质泥页岩等超固结黏土时具有明显的局限性。例如，修正剑桥模型计算得到的超固结黏土的不排水抗剪强度明显偏大，且应力—应变曲线为一条不光滑的曲线。另一方面，修正剑桥模型的剪胀方程在描述某些特定类型土的试验结果时存在一定的偏差，同时，无法考虑超固结比和温度等因素对土的剪胀特性的影响[10,15]。为了克服上述不足，Yamakawa等[16]提出了下加载面的力学概念，在加载的任意过程中，假定土体当前应力点始终位于下加载面上，发生弹塑性体积变形，物理意义明确，能连续平滑地刻画超固结土的弹塑性应力—应变关系。Yao等[17-19]定义了统一硬化参量，建立了适用于超固结重塑土的统一硬化（Unified hardening,UH）模型并取得了成功。统一硬化模型能够较好地描述超固结黏土的应变软化和剪胀特性；同时，在加载过程中，当土体由超固结状态转变为正常固结状态时，统一硬化模型能够计算得到光滑连续的应力—应变曲线。

笔者将炭质泥页岩视作超固结黏土，并且基于β-屈服函数建立了能够描述其应变软化和剪胀等复杂力学行为的修正统一硬化模型。通过将模型计算结果与炭质泥页岩和太古石的排水三轴压缩试验结果进行对比分析，验证了模型的合理性。

1 修正统一硬化模型

假定软岩的体积变形以压为正，其中，有效应力p、剪应力q和剪应力比η分别定义为

式中：σ1、σ2和σ3为主应力。

1.1 β-屈服函数

在三轴压缩条件下，修正剑桥模型的塑性功增量dWin可以表示为

式中：dεpv和dεps分别为塑性体积应变增量和塑性剪切应变增量；M为临界状态应力比。

根据上述能量耗散方程，可以分别得到修正剑桥模型的剪胀方程和塑性势。

式中：d为剪胀因子，可以定义为d=dεpv/dεps；pc为前期固结应力。由于修正剑桥模型采用相关联的流动法则，式（6）同样可以用于描述土的屈服面。

在描述超固结黏土的剪胀特性时，修正剑桥模型的剪胀方程具有一定的局限性，即低应力条件下该剪胀方程预测的剪胀因子d明显小于土的实测值。此外，由该剪胀方程得到的屈服函数假定土的物理屈服面具有椭圆形几何形状，与部分土的试验结果不符。

为了克服上述不足，首先将修正剑桥模型的剪胀方程改写为

式中：β为材料参数。令β=2，式（7）将直接退化为修正剑桥模型的剪胀方程。

如图2所示，与修正剑桥模型的剪胀方程相比，该剪胀方程能够较好地描述太古石的试验结果[7]，尤其是在低应力水平条件下。

图2 炭质泥页岩剪胀曲线Fig.2 Dilatancy curve of carbonaceous mud shale

通常可以根据剪胀方程得到塑性势并确定土的塑性应变率，在p-q应力空间，塑性势具有式（8）所示积分表达式[14]。

将式（7）代入式（8）并假定相关联流动法则，可得到式（9）所示的β-屈服函数。

为了便于探究材料参数β对屈服面几何形状的影响，可以将β-屈服函数改写为

根据式（10），当β=2时，可以得到修正剑桥模型的屈服函数。

材料参数β对屈服面几何形状的影响结果图3所示。由图3可以看出，增大参数β将导致计算得到的屈服面出现明显收缩，其弹性区间会随之减小。

图3 参数β对屈服面几何形状的影响Fig.3 The influence of parameter β on the geometry of the yield surface

β-屈服函数关于有效平均主应力p和剪应力q的偏导数分别为

将式（11）、式（12）代入相关联的流动法则中，可以计算得到炭质泥页岩在剪切过程中的塑性应变率。

1.2 相似因子

炭质泥页岩具有相对较好的连续性，其变形特征与典型的弹塑性材料类似，即在围压很小的情况下表现为脆性破坏，而随着围压的增大，将逐渐转变为延性破坏[7]。

在加载的初始阶段，炭质泥页岩的应力—应变曲线具有线性增长的特点，此时试样仅发生弹性变形。进一步加载将导致试样内部的胶结作用逐渐破坏并形成微裂纹，这些微裂纹的发展和贯通将导致试样发生屈服，其应力—应变曲线的非线性变化特点也将逐渐增强，表现为应变硬化。在达到峰值强度后，试样的剪应力会逐渐减小，其应力—应变曲线表现为应变软化。

可以看出，在复杂应力条件下，炭质泥页岩的应力—应变曲线变化特点与超固结黏土十分类似，均表现为应变硬化和应变软化。为此，将炭质泥页岩视作超固结黏土，并且定义其超固结比OCR为抗剪强度fc与围压pic之间的比值，即OCR=fc/pic。

参考姚仰平等[19]提出的适用于超固结黏土的统一硬化模型，假定在任意加载过程中炭质泥页岩当前的应力状态点A（p,q）始终位于下加载面上，其参考应力状态点B（pˉ,qˉ）则位于参考屈服面上，如图4所示。

图4 下加载面和参考屈服面力学示意图Fig.4 Schematic diagram of the mechanics of lower loading surface and reference yield surface

在p-q平面，下加载面和参考屈服面之间具有几何相似性，两者之间的相对位置关系可以采用相似因子R进行描述，R的定义为

根据式（13），当R=1时，下加载面会与参考屈服面重合，此时炭质泥页岩处于正常固结状态。R与OCR具有类似的物理意义，均可以用于描述炭质泥页岩的超固结性。对于三轴压缩试验，剪切伊始，不同围压条件下炭质泥页岩的初始相似因子R0为

采用本文提出的β-屈服函数，参考应力状态点B（）所在的参考屈服面可以表示为

当应力状态点位于参考屈服面上时，炭质泥页岩将处于正常固结状态，因此，类似于修正剑桥模型，可以采用体积硬化法则来描述硬化参量p0-pc的演化规律。

式中：cp=(λ-κ)/(1+e0)，其中，λ、κ分别为炭质泥页岩的压缩模量和回弹模量；e0为炭质泥页岩的初始孔隙比。

将式（16）代入式（15）中，可以得到

整理式（17），可以得到任意加载时刻有效平均主应力的表示式

根据相似因子R的定义，可以进一步得到R的表达式

可以看出，R与应力水平、应力路径和累积塑性体积应变等状态变量有关，式（19）能够综合反映上述影响因素对下加载面大小变化规律的影响。

1.3 统一硬化参量

由于当前应力状态点所在的下加载面几何相似于参考屈服面，因此，其屈服函数可以表示为

式中，统一硬化参量H的定义式为

根据Yao等[18]的研究，Ω可以表示为

式中：Mf为潜在峰值强度剪应力比。

在统一硬化模型[17-19]中，Mf被认为是关于相似因子R和临界状态剪应力比M的函数，满足关系式

式中，材料参数k可表示为

统一硬化参量Ω和剪应力比η之间的关系如图5所示，可以看出，η=0时，Ω=1；η=M时，Ω=0；η=Mf时，Ω=-∞。

图5 剪应力比η与统一硬化参数Ω的关系曲线Fig.5 Relationship curve between shear stress ratio η and unified hardening parameter Ω

在模型数值实现过程中，需要确定炭质泥页岩的体变模量K和剪切模量G，参考修正剑桥模型，可以采用式（25）计算。

式中：ν为泊松比。

2 参数试验

在统一硬化模型的基础上，采用了β-屈服函数来描述炭质泥页岩的剪胀特性，由此建立的模型能够合理地描述应变硬化和软化特性，同时考虑了围压对炭质泥页岩力学行为的影响。

模型包含的材料参数为M、e0、ν、λ、κ和β。其中，e0、λ、κ的具体取值可以通过炭质泥页岩的等向压缩曲线得到，临界状态应力比M则可以由炭质泥页岩在p—q平面的临界状态线加以确定。不同类型岩土材料的泊松比ν的取值范围通常为0.2～0.3，可以由特定围压下材料的初始应力—应变曲线得到。

2.1 超固结比对模型计算结果的影响

为了初步验证所建立的适用于炭质泥页岩的统一硬化模型的适用性，并且探究初始超固结比对炭质泥页岩力学特性的影响，开展了围压为200 kPa的排水三轴压缩数值试验，试样的初始超固结比分别为1、2、4和8，材料参数如表1所示。

表1 数值试样的材料参数Table 1 Material parameters of carbonaceous mud shale numerical samples

模型计算结果如图6所示。由图6可以看出，随着初始超固结比逐渐增大，试样的应力—应变曲线表现出明显的应变软化特性，同时，其体积变形将逐渐由剪缩变为剪胀。

图6 初始超固结比对炭质泥页岩排水三轴压缩试验结果的影响Fig.6 Influence of initial over-consolidation ratio on the results of drainage triaxial compression test of carbonaceous shale

2.2 材料参数β对模型计算结果的影响

为了更好地描述炭质泥页岩的剪胀特性，对修正剑桥模型进行了适当的修正，并由此得到β-屈服函数。在上述参数试验的基础上，分别令β为1.5、2、2.5和3，计算得到超固结比为8的试样的三轴压缩排水试验结果，其中，加载围压为200 kPa。

根据图7，增大β的取值会导致试样的峰值强度逐渐增大，其对应的轴向影响则会相应减小。此外，随着β的逐渐增大，试样的最终体积剪胀量也将随之增大。

图7 参数β对炭质泥页岩排水三轴压缩试验结果的影响Fig.7 Effect of parameter β on the results of drainage triaxial compression test of carbonaceous mud shale

3 模型验证

3.1 炭质泥页岩

为了进一步验证模型的适用性，采用该模型对炭质泥页岩[20]的三轴压缩试验结果进行计算，材料参数如表2所示。假定该炭质泥页岩的前期固结应力为50 MPa，当围压分别为100、200、300、400 kPa时，其初始超固结比分别为500、250、165、125。

表2 炭质泥页岩模型参数Table 2 Model parameters of carbonaceous mud shale

在不同围压下，试样的应力—应变曲线试验结果和模型计算结果如图8（a）所示。对比结果表明，该模型能够较好地描述炭质泥页岩的强度和变形特点，尤其是达到峰值强度后的应力—应变曲线，此时炭质泥页岩表现出明显的应变软化特性。不过，在低围压条件下，模型计算得到的峰值强度明显小于试验结果，峰值强度对应的轴向应变也明显偏大，可以调整部分材料参数的取值来进一步提高模型的计算精度。

图8（b）给出了模型计算得到的炭质泥页岩的剪胀曲线，可以看出，在加载的初始阶段，试样会发生剪缩变形，而进一步加载将导致试样发生剪胀破坏。

图8 不同围压条件下炭质泥页岩排水三轴压缩试验结果与模型计算结果对比Fig.8 Comparison of the results of triaxial compression test and model calculation of carbonaceous shale under different confining pressures