魏 炯，高富强

(1.中煤科工开采研究院有限公司，北京 100013；2.煤炭科学研究总院 煤炭资源高效开采与洁净利用国家重点实验室，北京 100013)

岩石力学是一门研究岩石材料力学性能的理论和应用学科，尤其注重于研究岩石的内力、变形、空间和时间之间的相互关系。随着计算机和计算数学的发展，科学计算继理论研究和科学实验之后，成为现代科学研究的第三种方法[1]。在过去数十年中，数值计算岩石力学蓬勃发展，采用数值计算研究岩石在载荷作用下的变形与破裂等力学行为成为岩石力学研究的热点[2-3]。

数值计算可以理解为用计算机来做实验，数值计算结果应与室内物理实验结果存在一一对应的关系。但由于计算能力的限制，一般将三维问题简化二维问题，将准静态问题简化为静力问题。事实上，自然界所有的物质都处于由空间和时间构成的四维时空之中，任何力学过程都是与空间和时间相关的过程。人们已经认识到岩石力学问题是三维的力学问题，将计算岩石力学从二维扩展到了三维，然而准静态加载条件下的时间维度依然为人们所忽略[4-9]。一般情况下，高速加载(地震、爆炸、冲击)条件下采用动力计算，而低速(准静态)加载条件下采用静力计算。

在进行岩石单轴压缩实验时，加载速率是必须输入的边界条件，一般采用载荷控制(kN/min)或位移控制(mm/s)。然而，在有限元数值计算中，一般将其简化为静力问题，按照kN/step或mm/step进行加载，加载速率与时间无关[10-16]。这种处理方式虽然可以实现岩石损伤破裂全过程模拟并得到全应力应变曲线，但不能真实反映物理世界的变化。它与真实物理实验在时间上不存在对应关系。从严格意义上讲，岩石力学实验中没有静力问题，静力问题只是动力问题的简化。即便在准静态加载条件下，也应考虑时间变量。动力问题和静力问题(准静态问题)应具有统一的表达。

数值计算的核心问题之一便是建立反映物理世界本质的数学模型，即建立反映真实世界各物理量之间关系的微分方程及相应的定解条件。没有正确完善的数学模型，数值计算便无从谈起。针对动力学问题，人们已经建立了运动微分方程，可以用来解决地震、爆炸、冲击等问题。对于低速加载的准静态力学问题，运动微分方程及相应的定解条件是否适用，是一个值得讨论的问题。因为动力计算涉及到应力波的传播、透射和反射，这与静力计算完全不同，两者计算结果是否相等鲜有文献提及。本文采用数值计算，探讨准静态加载条件下岩石破裂过程动态分析的可行性，讨论其时间步长的选取原则。

1 准静态加载的动态分析验证

假设对标准方形试件(50 mm×50 mm×100 mm)进行单轴压缩试验，采用轴向载荷控制，加载速率为30 kN/min(0.2 MPa/s)。根据上述条件，建立一个50 mm×100 mm的平面应力模型，弹性模量取10 GPa，泊松比取0.35，密度取2 400 kg/m3。在试件两端设置钢垫板，垫板尺寸为50 mm×10 mm，弹性模量取200 GPa，泊松比取0.33，密度取7 850 kg/m3，模型示意图见图1。采用轴向载荷控制，上边界施加如图2所示的载荷，加载速率为30 kN/min(0.2 MPa/s)。下边界采用固定约束，左右两侧采用自由边界。

图1 岩石单轴压缩试验示意图

图2 边界应力-加载时间曲线

在瞬态求解器下，当t=5 s时，其计算结果如图3所示，x方向最大位移为0.89 μm，最小位移为-0.89 μm；y方向最大位移为0 μm，最小位移为-10.00 μm；x方向最大正应力为1.18 MPa/s，x方向最小正应力为-0.76 MPa/s；y方向最大正应力为-0.80 MPa/s，y方向最小正应力为-3.00 MPa/s。当t=5 s时，对应的上边界应力为1 MPa。在静态分析模式下，其它条件不变，上边界施加1 MPa静应力，其计算结果如图4所示，x方向最大位移为0.89 μm，最小位移为-0.89 μm；y方向最大位移为0 μm，最小位移为-10.00 μm；x方向最大正应力为1.18 MPa/s，x方向最小正应力为-0.76 MPa/s；y方向最大正应力为-0.80 MPa/s，y方向最小正应力为-3.00 MPa/s。

图3 瞬态求解器下计算结果(t =5 s)

图5给出了两种求解器计算的Mises应力对比结果(竖向中心线处)，两者完全重合。以上结果表明在轴向载荷控制加载下，动态计算与静态计算的结果是一样的。可以预见，对于任意时刻都可以得到相应边界条件下与动态计算结果相同的静态计算

图5 竖向中心线处Mises应力对比(t =5 s)

结果。因此，在准静态加载条件下选择动态求解器计算是可行的。应力波的传播、透射和反射对计算结果没有影响，动力问题和静力问题(准静态问题)具有统一性。

2 动态分析的时间步长讨论

在真实单轴压缩试验过程中，岩石是连续体，且处于连续实时变化之中。但计算机只能处理离散数据，在进行数值计算时，必须对求解域进行分割，离散为若干单元。涉及动态计算时，对时间维度也必须进行离散，设置一定的时间步长。动态分析过程中，时间步长的选取是十分重要的问题。Tedesco等[17]认为对于动力问题，在数值积分中，为了保证积分精度，最大时间步长必须有所限制。最大时间步长的数值与应力波速度和有限单元的尺寸有关。最大时间步长的选择要确保在此时间内应力波能传过单元内高斯积分点之间的距离，该时间步长定义为：

(1)

式中:Lc为单元在波传播方向的长度，Cd为纵波波速。当Δt≤1/3(Δt)max时，可以保证计算精度。

对于高速加载(高应变率)力学过程，这种处理方式是合理的。但对于准静态加载(低应变率)力学过程，这种处理方式有待商榷。一般岩石单轴压缩试验持续时间为10 min～20 min。如果采用上述时间步长，电脑内存的消耗是巨大的，而且计算时间是漫长的。以本文单轴压缩数值计算为例，计算机处理器为英特尔酷睿i5处理器，内存16 G，三角形单元数目20 992个，自由度数目84 802，终止计算时间5 s。图6给出了不同时间步长下计算耗时。当时间步长为1 s时，计算耗时为12.33 s；当时间步长为1×10-1s时，计算耗时为14.06 s；当时间步长为1×10-2s时，计算耗时为27.96 s；当时间步长为1×10-3s时，计算耗时为163.03 s；当时间步长为1×10-4s时，计算耗时为2 424 s；当时间步长为1×10-5s时，求解运行24 h后，计算机内存溢出，求解失败。

图6 不同时间步长下计算耗时

可见，时间步长对计算效率的影响是巨大的。所以应该选择既能保证计算精度，又使计算量尽量小的时间步长。时间步长选取不但与单元尺寸和波速有关，而且与加载速率相关。不同加载速率或应变率(爆炸、冲击、普通压力机加载、蠕变)应该选择不同的时间步长。图7给出了单轴压缩数值计算中不同时间步长下竖向中心线的Mises应力。从图中可以看出，在时间步长分别取1 s、0.1 s、0.01 s、0.001 s、0.000 1 s时，得到的计算结果相同。这说明针对具体问题，必须具体分析，时间步长不能一概而论。在准静态加载条件下，选择较大的时间步长依然可以保证计算精度。

图7 不同时间步长下竖向中心线处Mises应力(t =5 s)

3 准静态加载的高速计算与观测

高速摄像机能将肉眼无法分辨的高速世界赋予新生命，显示肉眼看不见的瞬间动作。在岩石动力学试验中，有学者成功采用高速摄影观测到岩石的破裂过程[18]。但对于准静态加载力学过程，想要在室内实验室成功实现高速摄影是困难的。岩石的波速一般在2 mm/s～5 mm/s，为了观察其传播过程，数据采集频率必须在微秒级。否则，就无法保证有效捕捉到裂纹扩展过程。当摄像机运作如此快速时，大功率的照明设备是必需的，而且摄影成员必须使用烤箱手套来操作聚焦透镜。

如果摄像机以1帧/μs的速度记录图片，假设一张图片是500 KB，乘以一百万就是476.84 GB，而且这只是一秒的数据，而岩石试验一般持续时间为10 min～20 min。因此，它需要有一个惊人的储存器，而且每张图片像素不能太高。可以想象，如此高性能的相机，技术要求颇高，而且一定价格不菲。最重要的是，摄像机只能观测岩石表面结构，无法观测应力和位移的变化。

在数值模型中，边界条件类似于力学试验机，时间步长类似于摄像机。这为我们研究准静态加载下岩石破裂问题提供了一个很好的手段，可以得到任意时刻岩石内的应力和应变分布。图8给出了在0.2 MPa/s的加载速率下，不同时刻竖向中心线处米塞斯应力，t表示时间步长。

图8 不同时刻竖向中心线处Mises应力

图8(a)的时间步长为1 μs，在1 μs～5 μs时刻，岩石内部应力有明显的动态效应，岩石内部应力随着时间的推移而增加；图8(b)的时间步长为10 μs，也可以观测到岩石内部的动态效应。这种动态效应在应力云图中展示地更加清楚(见图9)，应力从上边界逐渐向下传递。但图8(c)—图8(d)中岩石的动力效应消失，这是由于加载速率很低，岩石的波速约为3 313 m/s，应力波从上边界传播到下边界约需要36 μs，这段时间只能产生7.2 Pa的应力增量，相比于已有的内部应力，其数值增量太小。因此，准静态问题可以简化为静态问题。但值得注意的是，惯性效应是否可以忽略，不但与加载速率相关，而且与研究对象的尺寸有关。如果研究对象尺寸无限大，即使在低加载速率下，惯性效应也不能忽略。

图9 不同时刻Mises应力云图

4 时间步长对岩石损伤破裂的影响

需要注意的是，岩石在加载初期是连续的，但在中后期会发生损伤、破裂直至失稳破坏，它是一个从连续到非连续的过程。如果我们用高速摄影拍摄岩石单轴压缩试验过程，无论我们采用1帧/s，1帧/μs，还是1帧/ns的拍摄速度，试验的最终结果是一样的，不会因拍摄速度的不同而变化。时间步长类似于拍摄速度，如果计算程序可以真实地再现岩石损伤与破裂全过程，那么计算结果也应该与时间步长无关。但事实上是有关的，数值计算永远都不可能做到真实，这是由计算机的特性决定的。因为现实世界的空间和时间都是连续的，而计算机只能处理离散数据，它所构建的虚拟世界必然是非连续的。任何模拟计算都只能做到“形似”，无法做到“神似”。

作者所在的研究团队一直致力于岩石损伤与破裂过程的数值计算。通过考虑岩石的非均匀性，将复杂的宏观非线性问题转化成简单的细观线性问题。通过引入数学连续物理不连续概念，将复杂的非连续介质力学问题转化成简单的连续介质力学问题。本节应用此方法研究时间步长对单轴压缩岩石损伤与破裂的影响。计算过程中按照位移加载，加载速率为0.002 mm/s，时间步长分别取0.1 s、1 s、10 s、20 s、50 s、100 s、200 s、400 s。图10给出了其应变随时间变化曲线。从图中可以看出，岩石应变随着时间匀速增加。这也从另一方面说明了动态计算的有效性。

图10 单轴压缩岩石应变随时间变化曲线

图11和图12分别给出了单轴压缩岩石应力随时间变化和应力应变曲线。因为应变随着时间匀速变化，所以两者变化规律一致。当时间步长小于10 s时，时间步长对应力曲线影响较小，计算结果不随着时间步长变化而变化。当时间步长大于10 s时，时间步长对应力曲线影响显著，岩石的峰值应力和弹性模量都随着时间步长的增加而减小。但在同一时刻，不同时间步长下的计算结果基本重合。这是因为时间步长取值过大，从而跳过了关键数据节点，造成曲线失真。图13给出了单轴压缩岩石损伤随时间变化曲线。从图中也得出类似结论，当时间步长小于20 s时，不同时间步长下的损伤曲线基本重合。当时间步长大于20 s时，不同时间步长下的损伤曲线差别显著。

图11 单轴压缩岩石应力随时间变化曲线

图12 单轴压缩岩石应力应变曲线

图13 单轴压缩岩石损伤随时间变化曲线

因此，时间步长必须小于某一阈值。本文模型峰值应力的对应时刻tp约为280 s，10/280≈0.04，建议选取0.03tp作为时间步长。图14给出了不同时间步长下岩石破裂模式。从图中可以看出，岩石损伤与破裂受时间步长的影响较大。当时间步长小于20 s时，岩石总是沿一条剪切破裂带破坏，其破裂模式基本不受时间步长的影响。当时间步长大于20 s时，岩石破裂模式开始发生变化，岩石呈共轭剪切破坏。时间步长为50 s和100 s时，形成两条主破裂带。时间步长为200 s、300 s和400 s时，形成多条主破裂带。从岩石最终破裂模式也可以看出，时间步长取0.04tp比较合理。

图14 不同时间步长下岩石破裂模式

5 结 论

本文利用数值计算，证明了准静态加载条件下的岩石破裂过程动态分析的可靠性，主要结论如下：

(1) 在准静态加载条件下选择动态求解器计算是正确可行的。应力波的传播、透射和反射，对计算结果没有影响，动力问题和静力问题(准静态问题)具有统一性。惯性效应是否可以忽略，不但与加载速率相关，而且与研究对象的尺寸有关。

(2) 在动力计算中，时间步长选取不但与单元尺寸和波速有关，而且与加载速率相关。不同加载速率或应变率(爆炸、冲击、普通压力机加载、蠕变)应该选择不同的时间步长。在准静态加载条件下，选择较大的时间步长依然可以保证计算精度。

(3) 时间步长对岩石的应力应变曲线和损伤与破裂模式影响较大。时间步长取值过大将造成曲线失真和破裂模式的不确定性。当时间步长小于0.03tp时，岩石应力应变曲线和损伤与破裂模式趋于稳定。