弹塑性软化煤柱的承载能力及其有限元模拟

徐晓惠，姚再兴

(1.辽宁工程技术大学矿业学院，辽宁阜新123000;2.煤炭科学研究总院矿山安全技术研究分院，北京100013; 3.北京大学工学院，北京100080)

相应的临界应力与临界应变分别是

在煤炭资源的地下开采中，为承受上覆岩层的压力，通常会形成长条形煤柱 (也叫矿柱)。煤柱的宽度不仅影响煤炭资源的采出率，也直接关系到巷道或采空区的安全。

传统的煤柱设计是以强度理论为依据的，只要认为煤柱某点进入屈服状态，便认为煤柱被破坏。显然，这种设计偏于保守。朱建明等人［1］推导了在SMP(Spatially mobilized plane)准则条件下煤柱的极限强度计算公式，认为煤柱极限强度计算应该考虑中主应力的影响，否则偏低。

突变理论的引入，使煤柱的设计纳入顶底板组成的系统。潘岳等［2］引入突变理论分析非对称开采矿柱的失稳。李江腾等［3］提出了非对称开采时矿柱失稳尖点突变模型。王连国等［4］基于尖点突变模型研究了矿柱失稳机理，提出了矿柱失稳判据。

王学滨［5］对矿柱的受力与稳定性进行了模拟，并考虑了不均匀性及矿柱端部的限制。张晓君［6］的模拟再现了采空区从变形到破坏的全过程。张少杰等［7］采用数值模拟、现场微震和应力监测，研究了工作面回采时，煤体内的应力分布及迁移规律，揭示了工作面的冲击矿压显现特征。

对于长条形的煤柱，在煤柱走向方向上各物理力学参数及受力和变形都没有变化，符合平面应变的条件。因此可以抽象为平面应变问题。假设煤柱材料的破坏满足Drucker－Prager破坏准则［8］。

本文把应力应变的非线性看作由材料的内变量决定，以经典弹塑性理论为基础推导了它们的关系。在积分表达时，采用了数值积分，并用具有软化特点的弹塑性有限元进行了模拟。

1 煤柱模型

1.1 煤柱模型简化

为了计算方便，对煤柱作如下简化:忽略自重及顶底板对煤柱两端的横向限制、煤柱的应力状态及应变状态在整个煤柱内是均匀的、对煤柱的应力及应变分析简化为对煤柱内任一点的应力及应变分析。图1是承受等效压力p的煤柱简化模型。

图1 简化的煤柱模型

1.2 应力特点

按照拉应力为正，且3个主应力由大到小的规定，即σ1≥σ2≥σ3，第一主应力σ1为水平向指向采空区，满足下式:

第二主应力σ2是煤柱的走向方向，为压应力;第三主应力σ3为煤柱横截面的压应力p，即

压应力p以压为正。主应力向量增量记作:

其中，t表示转置。

1.3 应变特点

与3个主应力及其方向分别对应的是煤柱的3个主应变。第一主应变ε1反映煤柱侧向膨胀;第三主应变ε3反映煤柱的压缩;煤柱的第二主应变ε2始终为零。因此，主应变向量增量记作

其中，t表示转置。

2 承载能力分析

2.1 弹性阶段解

2.1.1 弹性本构关系

在足够小的压力下，煤柱可以认为处于弹性受力状态，此时的弹性本构关系如下:

式中，D为弹性矩阵;σ，ε分别是应力向量与应变向量。

考虑到煤柱的受力与变形特点，可以解得弹性状态下的应力与应变关系如下:

式中，E是弹性模量;ν是泊松比。

2.1.2 初始屈服解

当煤柱从弹性状态进入临界塑性状态时，应力满足初始屈服面方程。把式(1)，(2)，(8)代入如下的Drucker－Prager屈服面方程

可求得煤柱进入屈服状态时临界压力ps为

式中，J2是应力偏量的第二不变量;I1是应力的第一不变量，式中的强度参数α，k在初始内变量处取得。I1，J2在主应力空间的表达式如下:

相应的临界应力与临界应变分别是

2.2 塑性软化分析

2.2.1 塑性本构关系

与Drucker－Prager屈服面方程相应的屈服面函数记作

材料进入屈服阶段的本构关系满足

式中Dep表示弹塑性矩阵，定义为

Dp为塑性矩阵，定义为

取塑性体应变θp作为内变量κ时，则

参数A定义为

塑性体应变θp的微分表达式为

式中，k'和α'分别是k和α对内变量κ的导数;C是弹性矩阵D的逆矩阵;G是剪切模量 (刚性模量)，G=E/(2(1+ν));K是体积模量，K= E/(3(1－2ν))。

2.2.2 强度参数的变化模型

通常实验室给出的是Mohr－Coulomb形式的强度参数，即黏聚力c和内摩擦角φ。通过这两个参数再等效给出Drucker－Prager的强度参数α，k。这里采用Drucker－Prager圆与Coulomb六边形外顶点重合的等效方法，即:

内摩擦角φ(弧度，对应的角度表示为φ0)保持常数，满足下列关系式:

黏聚力c对塑性体应变θp的依赖关系是

式中，cr表示残余黏聚力;cm为凸出黏聚力，它是最大黏聚力cM与残余黏聚力cr的差;为最黏塑性体应变，塑性体应变达到该值，黏聚力达到最大值cM;ξ称为胖度，胖度越大高斯函数曲线越缓; c0表示初始黏聚力，其变化曲线如图2所示。

图2 黏聚力随塑性体应变变化的高斯曲线

ξ与c0之一，及φ0，c0，cm，θp0由实验给出，ξ与c0的关系是:

式 (24)中的α'为零，k'的表达式如下:

2.2.3 软化限制的分析

在公式 (22)塑性矩阵Dp的定义中，参数A包含的材料弹塑性的重要性质。强度参数按式(27)和式 (28)给出时，k'＞0反映的是材料的强化;k'=0反映的是材料的理想塑性;k'＜0反映的是材料的软化。在材料的弹塑性模拟，特别是软化模拟时，把A控制在合理的范围内非常重要。A随内变量的变化而变化，在时，A取得极大值;在时，A取得极小值。图3是按实例分析中给定参数后A的变化情况。

图3 A随塑性体应变的变化曲线

由于A＞0，要求在A的极小值Amin时也能成立，即

需要注意，式中α是φ的函数。在数值计算中，Amin的合理性关系到计算能否进行下去或是否收敛的重要性。以下分析影响Amin取值的各个因素。

(1)Amin对ξ的导数

说明Amin关于ξ是增函数。

(2)Amin对cm的导数

说明Amin关于cm是减函数。

(3)Amin对φ的导数

式中，cm的数量级是106，K的数量级是109，ξ的数量级是 10－3。若成立，则关于 φ是增函数;若 cm成立，则关于φ是减函数。

(4)Amin对E的导数

说明Amin是对E的增函数。

(5)Amin对ν的导数

2.3 塑性软化的数值解

当煤柱受力进入塑性状态时，本构关系满足公式 (20)，其中的应力向量增量dσ与应变向量增量dε分别由公式(3)和公式(4)确定。由于煤柱的软化，不能事先给出煤柱的最大承压力，可以按控制应变方法进行求解。因此，对于给定的状态，弹塑性矩阵认为是已知的，第三主应变增量也是已知的，则可以求解dε1，dσ2和dp。

由于弹塑性矩阵是内变量的函数，与公式(20)对应的积分表达不易给出。数值积分的计算过程如下:

(1)首先依据公式(10)到公式(18)给出初始屈服状态的解，此时塑性体应变是零。

(2)依据公式 (21)计算弹塑性矩阵。

(3)任意给出第三主应变增量Δε3＜0(其绝对值越小，精度越高);依据与公式(20)相应的增量式Δσ=DepΔε，求出相应的增量Δε1，Δσ2和Δp。

(4)依据与公式 (25)对应的增量式Δθp=［1，1，1］(Δε－CΔσ)确定塑性体应变的增量，进而求出塑性体应变。

(5)如果还需要继续，转到第 (2)步。

采用这种方法，事实上给出一个弹塑性软化分析的精确解。这不仅给出煤柱承载能力，而且为软化弹塑性分析计算机程序给出一个测试算例。

3 煤柱承载能力实例分析

3.1 输入参数

残余黏聚力与凸出黏聚力的比cr/cm=0.35;峰值黏聚力内变量θp0=0.01;内摩擦角φ0=15°;胖度ξ=0.01;凸出黏聚力cm=1MPa。

对于如上给出的参数，A随内变量的变化见图3。可见，在这一算例中，A的最小值也是大于0的。

黏聚力随内变量的变化 (Gauss曲线)见图4。黏聚力先增加后减小，对应材料的强化和软化。

图4 黏聚力随塑性体应变的变化曲线

强度参数k随内变量的变化见图5。强度参数k先增加后减小，对应材料的强化和软化。

图5 k随塑性体应变的变化曲线

3.2 应变强化及软化曲线

按数值积分给出煤柱压力随压应变变化的曲线，见图6。

图6 煤柱压力随压应变的变化曲线

可见，煤柱在2.2MPa的压力下开始进入屈服状态;在约5.0MPa的压力下达到最大承载能力;煤柱的残余承载能力约1.4MPa。

4 弹塑性软化煤柱压缩有限元模拟

4.1 软化材料弹塑性有限元原理

与传统有限元程序相比，这一程序引入了弧长法，从而在延拓算法中不需要预先指定压力的增加或是减小，这就有可能使煤柱压力达到最大值时，煤柱的压缩继续增加，而煤柱的压力减小，煤柱进入不稳定平衡状态。

4.2 计算模型及控制参数

根据对称性，取煤柱的1/4建立有限元几何模型，如图7。左侧约束水平向位移;下侧约束竖向位移。上侧加均布压力。划分三角形网格，网络共计2690个，节点共计1416个。模拟得出的煤柱压力与煤柱压应变关系见图8。

图7 煤柱模拟的有限元模型

图8 软化有限元模拟的煤柱压力随压应变的变化曲线

可见，煤柱在2.5MPa的压力下开始进入屈服状态;在约4.8MPa的压力下达到最大承载能力;煤柱的残余承载能力约2.2MPa。

5 结论与展望

(1)在压力作用下，具有软化特性的煤柱开始表现出弹性性质，随后进入塑性状态。在强化阶段，承载的压力仍能增加;达到峰值承载能力后，如果要使煤柱保持平衡，必须减小压力，煤柱进入软化阶段，煤柱处于不稳定的平衡状态。

(2)煤柱按弹性极限或残余承载力设计是非常保守的，在经济上不合理，在能取得较准确的物理力学参数的基础上按峰值承载力设计可充分发挥煤柱的承载能力，但要取足够大的安全系数。

(3)对有软化特点的平面应变问题，对边受均布载荷，可以通过数值积分，得到足够精确的解，可以用作软化弹塑性分析计算机程序的测试算例。

［1］朱建明，彭新坡，姚仰平，等.SMP准则在计算煤柱极限强度中的应用［J］.岩土力学，2010，31(9):2987－2990.

［2］潘 岳，张 勇，吴敏应，等.非对称开采矿柱失稳的突变理论分析［J］.岩石力学与工程学报，2006，25(S2).

［3］李江腾，曹 平.非对称开采时矿柱失稳的尖点突变模型［J］.应用数学和力学，2005，26(8):1003－1008.

［4］王连国，缪协兴.基于尖点突变模型的矿柱失稳机理研究［J］.采矿与安全工程学报，2006，23(2):137－140.

［5］王学滨.具有初始随机材料缺陷的矿柱渐进破坏模拟［J］.中国矿业大学学报，2008，37(2):196－200.

［6］张晓君.矿柱及围岩对采空区破坏影响的数值模拟研究［J］.采矿与安全工程学报，2006，23(1):123－126.

［7］张少杰，王金安，魏现昊.煤柱过渡区的应力迁移与冲击矿压显现特征［J］.中国矿业，2013，22(5):73－78.

［8］Drucker，D.C.，Prager，W..Soil mechanics and plastic analysis for limit design［J］.Quarterly of Applied Mathematics，1952，10 (2):157－165.