红砂岩单轴压缩宏细观参数映射关系研究1)

佟 安 张军徽 武 娜

(北方工业大学土木工程学院，北京100144)

近年来，随着我国城市化水平的提高，城市地下空间开发利用得到大幅增长，同时也不可避免地发生着各类地质灾害。研究岩石的破坏机制对于地下工程施工和地质灾害治理等具有极其重要的意义[1]。经典的连续介质方法建立的岩石强度理论，难以解释岩石自身表现出来的非连续性(存在节理裂隙)、大变形以及各向异性等特征，以“均匀连续、小变形”为基本假设，建立的有限元软件不再适用[2]。1971 年，Cundall 最早提出适用于岩石特性的离散元法，1986年，王泳嘉将离散元方法引入国内，并不断发展至今。颗粒流分析程序(particle flow code，PFC)是一种基于离散元方法开发的用于模拟圆形颗粒介质的运动及其相互作用的数值分析软件，可以很好地模拟岩石的宏观力学行为[3]。

目前，PFC存在的主要问题是细观参数难以标定，原因在于极难获得颗粒细观参数与室内试验获取的宏观力学参数之间的关系，两者之间存在诸多难以控制的因素，并且这些因素之间的非线性关系极为明显[4]。多位学者使用“试凑法”解决此问题，但是盲目性较大且耗费大量时间。

反向传播算法(back propagation，BP)神经网络具有优秀的非线性动态处理能力，可实现高度非线性映射，具有较强的学习、存储、计算能力和容错特性，能实现输入样本和输出样本的非线性映射。本文提出以PFC2D 为基础获取宏细观参数样本，利用BP 神经网络建立宏细观参数映射关系的方式代替“试凑法”。

1 PFC2D砂岩单轴压缩试验模拟

PFC2D 是一种常被用于模拟岩石力学行为的二维离散元软件，使用刚性圆盘模拟岩样中的颗粒[5]。在选择接触本构模型时，平行粘结模型是最常用的接触本构模型。但此本构模型存在着明显的问题，其建立的岩石模型不能同时匹配抗压强度(uniaxia compressive strength，UCS)和抗拉强度(tensile strength，TS)，即UCS/TS值(一般为5左右)小于岩石的强度比值(10∼20)[6]，原因在于此模型构造圆形颗粒进行模拟，不能提供足够的抗扭转力，尤其是在模型产生裂缝后。平直节理模型可以解决此类问题，由图1 可以看出，该模型将圆形颗粒构造成多边形颗粒，以此抑制颗粒破坏后的旋转，使得UCS/TS值显著增大。

平直节理模型是一种描述有限尺寸、线弹性且考虑局部损伤的粘结模型，将交界面进行离散化，每个交界面都有粘结和未粘结2种状态，见图1。本文选用平直节理模型用于模拟红砂岩，试件大小选取50 mm×100 mm 的标准试样。为减少计算时间，力学模型尽量简化，由于岩石结构的复杂性，无法做到模拟颗粒与真实颗粒大小一致，当颗粒数超过3000后，模型力学性能基本保持稳定[7]，本文的PFC 模型颗粒总数为3602 个，满足要求。基础参数选取参考文献[8]，具体取值见表1。

图1 平直节理模型

表1 PFC 数值试样的基础细观参数

随着时间步长的增加，上下加载板匀速移动，挤压试件，当加载至应力下降到峰值应力的70%时，停止加载，此时砂岩试件已明显破坏，试验结束。

2 细观参数敏感性分析

本文针对红砂岩抗压强度偏小的特点，选取合适的细观参数范围进行单轴压缩试验的敏感性分析，来确定各个细观参数对试验的影响规律及程度，从而分类筛选出对砂岩单轴压缩应力–应变曲线影响程度较大的主要细观参数，为下一步选取BP 神经网络参数做准备工作。经过大量的敏感性分析，筛选出以下几个主要的参数：颗粒有效模量(E∗)、颗粒刚度比(κ∗)、平直节理有效模量平直节理刚度比抗拉强度(σ)和黏聚力(C)、平直节理摩擦系数(µ)，根据文献[9]，将与值保持一致，取值列于表2中。

表2 敏感性分析参数取值

2.1 弹性模量、刚度比敏感性分析

本文通过控制单一变量实现敏感性分析，保持表2 中各个初始值不变，从左至右依次同时由初始值改变的数值。同理，保持其他数值不变，调整分别分析接触弹性模量、刚度比对应力–应变曲线的影响规律及程度。

由图2可知，随着弹性模量的增加，单位应变条件下应力增长幅度呈减小趋势，表现出应力–应变曲线初始端越缓、斜率越小即宏观的弹性模量减小，但弹性模量的增加对峰值应力影响不大。同理，由图3可知，随着刚度比的增加，应力–应变曲线的弹性模量和峰值应力均增大。弹性模量与刚度比的变化对试样的应力–应变曲线作用显著。

图2 弹性模量敏感性分析

图3 刚度比敏感性分析

2.2 抗拉强度、黏聚力敏感性分析

从初始值开始依次向右改变表2 中的C取值，其他细观参数取值不变，同理改变σ的取值。分析黏聚力、抗拉强度对应力–应变曲线的影响规律及程度。

由图4 和图5 可知，黏聚力和抗拉强度的敏感性分析结果类似，随着数值的增加，峰值强度显著增大，而弹性模量基本保持不变。黏聚力、抗拉强度对试样的应力–应变曲线作用显著。

图4 黏聚力敏感性分析

图5 抗拉强度敏感性分析

2.3 摩擦系数敏感性分析

从初始值开始依次向右改变表2 中µ的取值，其他细观参数取值不变。分析µ对模型应力–应变曲线的影响规律及程度。

从图6 中可以看出，平直节理摩擦系数对应力–应变曲线的峰值应力和弹性模量影响并不明显，因此，将此变量取固定值，不作为BP 神经网络的输出参数。

图6 摩擦系数敏感性分析

3 建立宏细观参数映射关系

3.1 BP神经网络基本原理

本文采用的BP 神经网络是一种按误差逆向传播算法训练得到的多层前馈网络，在岩土工程领域应用较为广泛[10]。此方法的优势在于能够存储与学习数量庞大的输入输出映射关系，且无需映射关系之间有着明确的数学表达式。BP 神经网络的应用已十分广泛，其原理此处不多赘述。

3.2 BP神经网络建立的总体流程

从PFC2D 模拟红砂岩单压缩试验得到的应力–应变曲线上获取的宏观力学参数与细观参数建立BP 神经网络关系，主要分为6 个工作流程：(1)获取与细观参数相对应的宏观力学参数；(2)宏观力学参数作为输入样本，细观参数作为输出样本，用于训练BP神经网络；(3) 调整BP神经网络程序，使宏细观参数建立精确的映射关系；(4)建立BP 神经网络映射关系后，输入专门用于测试的宏观力学参数，反演相应的细观参数；(5)将反演出的细观参数输入到PFC2D中，进行数值模拟试验，并将数值模拟计算所得到的宏观力学参数与第4 步输入的宏观力学参数进行对比，验证BP 神经网络的精准性；(6)不断重复上述工作流程，直到建立出较为精准高效的BP神经网络。

3.3 BP神经网络的创建

首先，选取建立网络的宏细观参数。采用“试凑法”基于模拟相似应力–应变曲线的原理，试用4 个宏观参数描述应力–应变曲线作为网络的输入样本，分别为：弹性模量、峰值处的应力、应变值以及破坏斜率，如图7。此4 个宏观参数是针对PFC2D 模拟红砂岩单轴压缩试验应力–应变曲线的特点总结，其弹性段和最终破坏段，都具有相当明显的线性增长(下降)的趋势，便于精确获取且利于网络寻找宏细观参数映射关系。使用6 个细观参数作为神经网络的输出样本，分别为：

图7 基于应力–应变曲线选取宏观力学参数

利用PFC2D内置的FISH以及Python语言，二次开发红砂岩单轴压缩程序，使其可以自动在6 组细观参数的取值范围内产生随机数，见表3。并记录相应的宏细观参数，程序自动运行40 次作为网络输入输出的样本库。

表3 PFC 细观参数取值范围

由于用于训练的输入样本和输出样本不同，需要进行归一化处理，便于快速寻找映射关系，归一化公式如式(1)所示。利用此式将训练样本数据，全部归一化至[0,1]之间。

随后创建BP 神经网络需要考虑网络层数及中间层神经元个数，一般情况下具有一个中间层的BP神经网络可以很好地模拟任何有理函数，故本文选取具有一个中间层的三层BP神经网络。4个宏观力学参数作为输入层样本参数，6 个细观参数作为输出层样本参数，其中取相同值，因此，输出层样本参数为4 个。最终，输入节点为4，输出层节点数为4。

中间层神经元个数的选取对整体性能影响较大，数目过少，网络将不能建立复杂的判断界，难以训练出合理的网络，兼容性较差；数目过多会导致训练时间较长，网络的泛化能力下降，计算误差可能变大。因此，选取合适的神经元个数至关重要，但目前并没有确定的科学办法选取神经元个数，一般根据如下经验公式[11]来确定

其中，p为中间层神经元个数，n为输入层节点数。本文根据式(2)选神经元个数为5个左右，在此范围上下浮动5 个神经元，选出期望值与反演输出值均方差最小的神经元个数，最终确定神经元个数为9，如图8所示。

图8 神经元个数选取对比

最后选取建立网络的必要函数，网络输入层与隐藏层及隐藏层与输出层之间的传递函数采用logsig 函数和purelin 函数；考虑到网络的规模和学习时间，选用trainlm 函数，trainlm 函数使用Levenberg–marquardt算法，训练网络。

3.4 BP神经网络的训练与反演

利用之前获取的40 组宏细观参数作为训练样本集，训练次数为2000，训练的目标误差为0.008，学习率设为0.05。将训练好的网络用sim 函数进行反演模拟，并不断调试。

3.5 BP神经网络反演能力验证

再次使用3.3节中的方法，获取3组不同于训练样本集的宏细观参数来验证BP 神经网络的反演能力，如表4所示。

表4 BP 神经网络反演值与期望值对比

表4 中反演输出值与期望值之间总体误差保持在目标误差0.01 以内。并将第1 组反演输出结果输入至PFC2D中进行模拟，对比结果如图9所示，可以看出两曲线高度相似，表明创建的BP 神经网络具有较强的反演能力，且精度较高。

图9 模拟曲线与原模拟曲线对比

4 BP神经网络的室内试验验证

红砂岩试样取自武汉一采石场，严格按照规范[12]要求加工成50 mm×100 mm的长方体，并确保试样没有较大裂缝。本试验选用吉林省金力试验技术有限公司研制的YAW-2000 微机控制电液伺服压力试验机，加载方式为位移控制，加载速率为0.005 m/s,以上条件均满足PFC2D的参数选取。试验获取多组应力–应变曲线，来验证BP 神经网络的可行性。使用本文方法获取试验曲线的4 个宏观参数作为输入样本，输入到上述创建成功的BP 神经网络中，网络的反演输出得到的细观参数如表5 所示。再将得到的细观参数输入至PFC2D 单轴压缩程序中，模拟得到的应力–应变曲线与试验曲线进行对比。从图10中可以看出，PFC模拟曲线与试验曲线十分接近(由于砂岩试件内部结构复杂，分布并不均匀，试验曲线弹性段未能体现出线性增长)，进一步从数据上看，如表6 所示，输出的宏观参数值与室内单轴压缩试验所获取的宏观参数值极为接近，同时如图11所示两破坏形式也很相近。由此证明，本文所创建的BP神经网络系统是合理可行的。

表5 BP 神经网络反演的PFC 细观参数值

图10 PFC 模拟曲线与试验曲线对比

表6 PFC 模拟获取的宏观参数与单轴压缩试验参数对比

图11 PFC 试样与试件破坏形式对比

5 UCS/TS值验证

图12 为上述神经网络模拟室内试验最终识别出的红砂岩单轴压缩曲线，其UCS值为34.45 MPa。同样将网络输出的细观参数输入至PFC2D 红砂岩拉伸模拟试验中，获得的单轴拉伸应力–应变曲线如图13 所示，其TS 值为3.18 MPa。UCS/TS 值约为11，在合理范围内，表明平直节理模型可以正常模拟出岩石UCS/TS 的真实值，弥补了平行粘结模型在此方面的缺陷。

图12 PFC 单轴压缩曲线

图13 PFC 拉伸曲线

6 结论与展望

本文基于BP 神经网络并二次开发PFC2D 岩石单轴压缩程序来解决参数难以标定的问题，在BP 神经网络中输入应力–应变曲线上的宏观力学参数，反演出细观力学参数，再输入至PFC2D 中进行模拟，得到与试验曲线高度吻合的模拟曲线，得出以下结论：

(1)通过上述对比可以看出，BP 神经网络具有非线性映射能力、泛化能力等优势，可以精确地建立宏细观参数映射关系，反演得到的应力–应变曲线与室内试验曲线高度吻合，且本文创建的网络仅需要40组样本数据，较“试凑法”节省大量时间，表明BP神经网络可以很好地解决细观参数标定难题。

(2)本文选取的平直节理模型有效地解决平行粘结模型关于UCS/TS 过小问题，展现出平直节理模型模拟岩石的优势。

(3) PFC2D 与BP 神经网络的联合使用效果显著，但是过程中依旧存在一些不确定因素，比如神经元个数的选取依旧难以寻找最优解，目前只有经验公式可供选择，有待解决。

(4) PFC2D 模拟红砂岩的单轴压缩试验所得到的应力–应变曲线与室内试验获得的曲线仍存在差别，弹性段为直线，过于理想，程序有待再度开发。