5.3 基于围岩力学参数概率分布模型的变形敏感性灰关联分析

5.3.1 基于参数概率分布模型的灰关联敏感性分析方法建立

建立基于参数概率分布的灰关联分析模型，首先确定该模型的因素序列与目标序列，根据参数服从的相关概率分布模型将因素序列转化为累积分布概率序列，计算出因素序列与目标序列的关联度。最后，根据关联度大小进行力学参数敏感性评价，具体流程如图5.13所示。

图5.13 基于参数概率分布模型的灰关联模型流程

5.3.1.1 确定敏感因子因素序列矩阵及目标序列矩阵

设有m个影响因素，建立影响因素序列A，并选取各因素序列A所对应的围岩变形量作为目标序列B。

影响因素序列矩阵A：

式中：aij为第i个影响因素的第j个取值（n个应值遍布整个定义域）。

假设m个影响因素服从累积分布函数F1（x）、F2（x）、…、Fm（x），并结合因素序列A建立累积分布概率序列A′，令=F1（A1），则有

式中：a′ij为第i个影响因素在aij点的累计分布概率值。

此时，与因素序列A相对应的目标序列B为

式中：bij为aij所对应的围岩变形量。

5.3.1.2 矩阵的无量纲化

对参数进行敏感性分析需要消除不同参数之间的量纲差异，统一量纲才具有可比性，应用极差变化的方法对数据进行归一化处理，对A′进行数值变换可得

式中：。

同理，将目标序列B进行归一化运算，计算结果为B′。

5.3.1.3 关联度分析及敏感性评价

差异信息构成灰关联差异信息空间矩阵，差异信息的计算公式为

根据差异信息求得差异序列矩阵Δ，然后从中提取最大值和最小值，记为Δmax=max（Δij）、Δmin=min（Δij）。

通过计算比较点和参考点之间的距离分析各因素的差异性和相关性，用关联系数表示因素序列与目标序列之间存在的相关性，表达式为

式中：ζ为分辨系数，ζ∈［0，1］，本文取ζ=0.5，用来增加关联系数具有明显的差异性。

由式（5.42）计算关联度得

分析式（5.43）可知，关联度在区间［0，1］上的取值是变量。在关联度矩阵中，关联度相对越大，其敏感性越大；反之，敏感性越小。

5.3.2 3DEC离散元软件计算节理隧洞围岩变形量

5.3.2.1 3DEC离散元简介

3DEC（Three Dimension Distinct Element Code）是基于离散元基本理论由Cundall和Hart于1985年开始正式开发的用于解决离散介质不连续问题的三维分析软件，其延续了二维的UDEC的核心思想，具备静力学、热力学和动力学三类不同问题的求解模块，同时对流体问题有一定的考虑。另外，还开发提供了FISH高级程式语言，可以自由定义变量和函数，且为命令驱动模式的分析软件，因而为使用者提供了极大的自主空间。该方法计算简单，对计算机内存及硬件要求较低，但需要较长的迭代时间。3DEC离散元数值模拟的基本假设和计算原理主要有以下4个方面：

（1）3DEC离散元程序朴素地看待岩体等介质离散构成特征，将其视为连续性特征（块体等）和非连续性特征（结构面等）这两个基本元素的集合体，并以非常成熟的最基本的牛顿第二定律定义其力与位移变化情况。

（2）3DEC采用凸多面体来描述所在介质中的具有连续性特征的对象的空间形态，并通过凸多面体组合的方式实现凹形的连续对象；另外，非连续特征以若干三角网所组成的曲面进行表示。

（3）连续性特征对象的凸多面体可服从刚体或变形体的受力与变形规律，当为可变形体时，其采用与有限元分析程序一致的快速拉格朗日分析法进行求解。

（4）连续性对象特征之间的相互作用是通过非连续特征来实现的，描述非连续特征对象的受力和变形可遵照接触定律，而模拟凸多面体间在公共接触面上发生滑移或脱离的行为则遵从力学定律。

3DEC采用具有众多优点的显示计算方法（有限差分法），其是一个基于运动方程的解决方案，能够有效地表明不连续系统的潜在破坏模式。计算过程中，每一个时步（Timestep）都会用到本构方程和运动方程，3DEC的计算循环过程如图5.14所示。

图5.14 3DEC计算循环过程

3DEC在求解问题时的一般过程主要包括6个步骤。

第一步：数值计算模型的建立。

通过直接建立基础模型并将其切割或使用PGEN生成或其他三维建模数据导入3DEC中，形成符合实际情况的形体；然后采用JSET等命令根据实际情况生成一系列不连续面。

第二步：设置材料属性。

赋予块体和结构面基本的物理力学参数，块体在应力较低时可考虑成变形体，并可采用GEN EDGE NUM命令为块体自动划分网格。

第三步：设置边界条件和初始条件。

上述基本模型建立完成后，可通过BOUNDARY命令施加力学边界条件，采用INSITU命令施加初始原岩应力条件（注意速度边界条件只能施加在可变形体上）。

第四步：调整至初始平衡。

程序在正式求解前要达到应力平衡，采用差分格式解出一个微小时段的速度和位移，逐渐通过STEP命令增加时段，达到平衡状态。是否达到平衡可以通过对不平衡力的监测得到结果，当不平衡力远远小于所施加外部载荷时，我们就认为模型已经达到平衡。

第五步：进行问题计算求解。

3DEC允许模型在求解过程的任意点处，状态发生变化，实际求解时可根据需要改变荷载、应力和材料属性等以获得目标结果。

第六步：计算结果的存储及调用。

3DEC中采用SAVE和RESTORE命令可以有效地保存和调用计算中任意时刻的结果文件，有助于节省运算时间。

5.3.2.2 3DEC计算模型的建立

依据工程需要，采用美国ITASCA公司开发的3DEC三维离散元计算软件建立新疆布仑口-公格尔水电站地下洞室3号标段地质模型，进行围岩力学参数拱顶位移敏感性研究。三维模型尺寸为60m×60m×60m，城门洞隧洞宽4.6m，高5.3m，划分单元66037个，选取埋深800m。边界条件为：①模型左右边界和下边界施加位移约束；②模型上部边界施加21.6MPa的垂直压力。屈服准则采用摩尔-库仑准则。模型中考虑3条遍布整个模型的节理组J1、J2和J3，采用摩尔-库仑滑动屈服准则，依据勘查设计资料得节理参数见表5.6，地下洞室三维模型如图5.15所示。

5.3.2.3 影响参数及考核指标的确定

本文借助三维离散元软件对地下洞室围岩稳定性参数进行拱顶位移敏感性分析，在进行最不利断面选择时，发现Y=15.0断面处拱顶位置处于结构面附近且围岩下沉量较其他断面大，故将Y=15.0m断面为研究对象，拱顶（0，15，2.3）处位移作为考核指标，具体位置见图5.16。选择岩体密度ρ、弹性模量E、泊松比μ、黏聚力c、内摩擦角φ及节理内摩擦角φj等6个参数为拱顶位移敏感性影响参数，由于缺乏详细的勘察资料，在对布仑口-公格尔水电站某标段围岩地质勘查设计资料和相关试验数据进行整合分析的基础上，与相关类似工程进行类比统计，综合分析得该标段围岩力学参数基准值及分布概率模型，具体见表5.7，6个参数的统计频率直方图及概率密度曲线见图5.17。

图5.15 围岩节理岩体模型

图5.16 Y=15.0m断面示意图（单位：m）

由图5.17可知，各个参数统计直方图与正态分布累积概率密度曲线图吻合度较好，故选择正态分布形式满足要求。

5.3.3 基于参数概率分布模型的变形敏感性灰关联分析计算

将6个参数分别在定义区间范围内变化，其他参数固定在基准值，求得拱顶下沉量以及累积分布概率，计算结果见表5.8～表5.13。

图5.17 参数统计频率直方图及概率密度曲线

由表5.8～表5.13可知，围岩力学参数累积分布函数矩阵为

拱顶位移目标因素矩阵为

根据式（5.40）对影响因素累积分布函数矩阵A′和拱顶位移量目标因素矩阵B进行无量纲化处理得

根据式（5.42）求得关联系数矩阵为

根据式（5.43）求得关联度为

为了对比验证基于影响因素概率分布模型的灰关联变形敏感性的合理性，对不考虑影响因素概率分布模型的常规敏感性灰关联分析，将二者计算得到的参数敏感性结果列于表5.14。

从表（5.14）可知：①基于影响因素概率分布模型的灰关联变形敏感性程度大小排序为：密度ρ、弹性模量E、黏聚力c、泊松比μ、内摩擦角φ及节理内摩擦角φj；②常规灰关联变形敏感性程度大小排序为：密度ρ、黏聚力c、泊松比μ、弹性模量E、内摩擦角φ、节理内摩擦角φj。二者计算结果出现不一致情况，而考虑影响参数的概率分布更加符合实际地质参数的空间变异性，减小由于参数选择不当带来的计算误差。因此，本节方法能够更加准确、合理地对参数进行客观评价。