本发明属于岩石力学与岩石工程领域，具体涉及有限元-离散元耦合数值模拟程序(FDEM)的输入参数标定方法及被标参数的可靠性验证方法。

背景技术：

岩石的破裂问题是研究者及工程师需要重点关注的问题，因为岩石破裂后会直接影响到如隧道围岩及岩石边坡等的稳定性。岩石的破裂经历了从弹性变形到塑性变形，再到岩石破裂后块体的相互接触及块体大运动全过程，一方面岩石材料复杂，天然岩体含有大量微裂隙，且实际工程岩体往往以复合岩体的形式出现；另一方面，岩石的力学性质与其赋存环境、应力历史、应力路径和研究尺度等密切相关；再者，对岩石破裂的研究不仅仅是为了研究其破裂过程及影响因素，更重要的是研究如何促进岩石破裂(如矿石开采)和抑制岩石破裂以便控制其稳定性(如隧道围岩)。

面对如此复杂的材料性质、赋存环境和实际需求，理论解析显得力不从心，而无论室内试验或工程现场试验，只对特定工程有效、不具有预测性，且研究成本高昂、研究周期漫长。为了便于岩石破裂的研究，数值模拟方法应运而生，大致地，可分为连续性方法(如有限元法FEM，有限差分法FDM和边界元法BEM等)、非连续性方法(如颗粒元法PFC，通用离散单元法UDEC和不连续变形法DDA等)和二者的耦合方法(如数值流形法NMM及本说明中的FDEM法等)。其中，FDEM由于其计算效率高、能够重现岩石弹性-塑性-破裂全过程等优点得到了广泛应用，包括室内试验模拟和实际工程的模拟研究。

准确的输入参数是保证数值模拟结果可靠性的前提，遗憾的是，目前尚无高效便捷的FDEM输入参数标定方法。目前普遍采用单轴压缩和巴西劈裂模拟试验不断试算不同的参数组合，以期得到与室内试验结果相吻合的输入参数，需要重复高达上千次的试算方能得到较为理想的输入参数组合。更为严重的是，该输入参数具有严重的网格尺寸依赖性，通过室内试验标定得到的参数直接应用于工程尺度的模拟是存疑的。因此，本发明旨在提出一种快速便捷的FDEM输入参数标定方法，消除参数取值对网格尺寸的依赖性，并提出被标参数的可靠性验证方法。

技术实现要素：

为了解决上述技术问题，本发明提供一种快速的FDEM输入参数标定方法，并消除参数取值对网格尺寸的依赖性，同时提出被标参数可靠性验证方法。

本发明提供的技术方案如下：

一种有限元-离散元耦合数值模拟程序输入参数快速标定方法，包括以下步骤：

(1)基于室内岩石力学实验，获取岩石宏观物理力学参数，包括单轴抗压强度σc，弹性模量E，泊松比ν，抗拉强度σt，粘聚力c，内摩擦角和岩样的破坏模式；

(2)采用不同的节理罚值Pf，进行标准岩样单轴压缩模拟试验，监测峰前应力-应变曲线，得到数值模拟的弹性模量值并与输入值对比，通过不断试算得到节理罚值Pf的取值范围；

(3)采用两个三角形单元的压缩-剪切破坏模拟试验，标定得到法向接触刚度Pn取值；

(4)采用单轴压缩模拟试验，设置抗拉强度σt为较大值，以使岩样仅发生剪切破坏而不发生拉伸破坏，获得II型断裂能GII；

(5)采用标准岩样直接拉伸试验，获得I型断裂能GI；

(6)根据步骤2-5得到的节理罚值Pf、法向接触刚度Pn、II型断裂能GII和I型断裂能GI，再次采用单轴压缩模拟试验，变化不同的切向接触刚度Pt，监测岩样的破坏率(破坏率被定义为破坏节理单元数目与总节理单元数目的比值)，直至获得稳定的破坏率，获得切向接触刚度Pt的取值范围；

(7)根据步骤1-6得到的宏观输入参数和微观输入参数，进行三轴压缩、巴西劈裂和直接剪切模拟试验，验证步骤(2)-步骤(6)所标参数的可靠性。

进一步，所述步骤(1)中岩石力学实验包括单轴压缩、直接拉伸、巴西劈裂和三轴压缩试验。

所述步骤(2)的输入值为室内试验值。

进一步，所述步骤(2)中试算方法如下：采用不同的节理罚值Pf，进行单轴压缩模拟试验，得到对应的峰前弹性模量Eout，直至该比值与输入弹性模量Ein的比值接近于1为止，得到合理的Pf取值范围；变化不同的网格尺寸h和输入弹性模量Ein，验证Pf取值范围的鲁棒性，得到通用的Pf值取值范围。

进一步，所述步骤(3)中，两个三角形单元之间能够由粘接状态平稳过渡到接触状态。

进一步，所述步骤(3)中法向接触刚度Pn计算公式，如下：

式中，Pn(i-j)为三角形单元i与三角形单元j间的法向接触刚度，hi、hj分别为三角形单元i、j的高，li、lj分别为三角形单元i、j的长，Pb为基本刚度，采用下式计算：

Pb＝αPf

式中，α为系数。

进一步，所述步骤(4)中II型断裂能GII的获得方法为：单轴压缩模拟试验后，宏观输入参数取为试验值(弹性模量E，泊松比ν，粘聚力c，内摩擦角)，节理罚值Pf和法向接触刚度Pn分别根据步骤2和步骤3的方法确定，切向接触刚度Pt设为较大值(如Pt＝10E)，变化II型断裂能GII，直至获取与室内试验相吻合的模拟结果。

进一步，所述步骤(4)中，修正了FDEM中现有的节理单元本构模型，将现有的应力-位移本构模型修正为应力-应变本构模型，消除了参数取值对网格尺寸的依赖性，如下：

式中，or，sr分别为节理单元的极限拉伸和剪切位移，fs为节理单元抗剪强度。通过上式表明，节理单元的破裂由应变控制，而不再是位移控制，消除了GI和GII取值对网格尺寸的依赖性。

进一步，所述步骤(6)中，岩样的破坏率达到稳定值时，认为切向接触刚度Pt是合理的，如下：

式中，η为岩样的破坏率，nfailed，ntotal分别为破坏节理单元数目和节理单元总数目。

进一步，综合采用了三轴压缩、巴西劈裂和直接剪切模拟试验验证所标参数的可靠性，因为上述三种试验的模拟结果可直接与理论结果对比，可便捷、准确判定模拟结果的正确性。

与现有技术相比，本发明提供了一种快速的FDEM输入参数标定方法及被标参数可靠性验证方法，并消除了参数取值对网格尺寸的依赖性，实现了输入参数的快速获取，提高了参数取值的准确性，为数值模拟结果的可靠性提供了基本保障，确保了所有参数取值来源有据可依。

除此以外，本发明还具有以下有益效果：

(1)所有的宏观输入参数可采纳为室内试验值，减少了被标参数数目；

(2)基于节理单元本构模型修正，消除了参数取值对网格尺寸的依赖性；

(3)所标定的节理罚值Pf、法向接触刚度Pn和切向接触刚度Pt取值范围具有强的适用性，对任意材料参数均是适用的；

(4)对I型和II型断裂能的标定效率高，约需要20次的标定试验即能得到较为理想的输入值；

(5)综合采用三轴压缩、巴西劈裂和直接剪切模拟试验验证了被标参数的可靠性。

附图说明

图1(a)为室内试验单轴压缩破坏的真实岩样；图1(b)为室内试验巴西劈裂破坏的真实岩样；

图2(a)为本发明涉及的单轴压缩模拟试验数值模型；图2(b)为本发明涉及的直接拉伸模拟试验数值模型；

图3为本发明涉及的三角形法向接触刚度计算示意图；

图4为本发明涉及的三角形法向接触刚度计算公式中系数α标定方法示意图；

图5(a)为本发明涉及的被标参数可靠性验证三轴压缩模拟数值模型；图5(b)为本发明涉及的被标参数可靠性验证巴西劈裂模拟数值模型；图5(c)为本发明涉及的被标参数可靠性验证直接剪切模拟数值模型；

图6为本发明提出的FDEM输入参数标定流程。

具体实施方式

下面结合具体实施例对本发明的内容进一步说明，本发明的内容完全不限于此。

实施例

以下结合附图对本发明涉及的有限元-离散元耦合数值模拟程序(FDEM)输入参数快速标定具体实施方案进行详细说明。FDEM输入参数的具体标定流程如下：

步骤1.将野外获得的岩心加工成标准岩样，根据国际岩石力学与工程学会(ISRM)推荐的方法进行室内单轴压缩(直径50mm×高度100mm)、直接拉伸(直径50mm×高度100mm)或巴西劈裂(直径50mm×高度25mm)和三轴压缩(直径50mm×高度100mm)试验，获得真实岩样的宏观物理力学参数，包括单轴抗压强度σc，弹性模量E，泊松比ν，抗拉强度σt，粘聚力c，内摩擦角和岩样的破坏模式。以某次室内试验为例，得到的岩石宏观力学参数如下：σc＝23.6MPa,E＝12.5GPa,ν＝0.25,σt＝2MPa,c＝7MPa,岩样破坏模式如图2所示。

步骤2.将输入弹性模量E和泊松比ν设定为室内试验值，即：Ein＝E，νin＝ν，根据图2(a)所示的模型，进行单轴压缩模拟试验，仅观测峰前应力-应变曲线，得到输出弹性模量Eout和输出泊松比νout，输出泊松比被定义为横向应变与轴向应变之比的绝对值，变化不同的节理罚值Pf，直至且为止，确定该Pf值为合适的节理罚值。变化图2(a)的网格尺寸和输入弹性模量，验证上述Pf值与网格尺寸的无关性及对任意输入弹模的可适用性。经过多次试算表明，Pf＝20-30E即可满足要求。

步骤3.根据图3采用下式计算三角形单元间的法向接触刚度Pn，以三角形i与j间的接触为例，其法向接触刚度Pn(i-j)为

式中，hi，hj分别为三角形单元i、j的高，li、lj分别为三角形单元i、j的长，Pb为基本刚度，采用下式计算：

Pb＝αPf

式中，α为系数。

采用图4所示的模型对任意两个三角形单元进行压剪模拟试验，实时监测6个节点的节点力，采用二分法获得合适的α值，直至在节理单元发生破坏瞬间，6个节点的节点力均实现了由粘结至接触的平稳过渡。变化不同的网格尺寸、三角形单元形状、加载速率、计算时步、弹性模量、II型断裂能，确定α取值不受上述参数的影响。经过试验表明，α＝0.1448能够满足要求。

步骤4.采用图2(a)的模型，进行单轴压缩模拟试验，将输入参数中的抗拉强度σt设定为极大值，使得岩样不发生拉伸破裂，其余宏观输入参数直接采用室内试验值(弹性模量E，泊松比ν，粘聚力c，内摩擦角)，节理罚值Pf和法向接触刚度Pn分别根据步骤2和步骤3确定，切向接触刚度Pt设定为较大值(如Pt＝10E)，不断变化II型断裂能GII，直至模拟得到的单轴抗压强度和岩样破坏模拟与室内试验结果接近，则确定为合适的GII取值。变化网格尺寸，采用上述获得的GII值再次进行单轴压缩模拟试验，验证该值与网格尺寸的无关性。

变化网格尺寸，确定GII取值不存在网格尺寸的依赖性，因为该依赖性可通过下式消除：

式中，or，sr分别为节理单元的极限拉伸和剪切位移，fs为节理单元抗剪强度。通过上式表明，节理单元的破裂由应变控制，而不再是位移控制，消除了GI和GII取值对网格尺寸的依赖性。经过多次试算表明，在本次模拟中，GII＝30000J/m2能够满足要求。

步骤5.采用图2(b)的直接拉伸模拟试验标定I型断裂能GI，抗拉强度取为室内试验值，不断变化GI，当模拟得到的岩样破裂模式和拉伸应力-应变曲线与室内试验结果一致时，则确定该GI值是合理的。采用不同的网格尺寸，验证GI取值与网格尺寸的无关性。经过多次试算表明，在本次模拟中，GI＝10000J/m2能够满足要求。

步骤6.再次采用图2(a)的单轴压缩模拟试验，将节理罚值Pf，法向接触刚度Pn，II型断裂能GII和I型断裂能GI分别设定为步骤2-步骤5得到的数值，进行切向接触刚度Pt的标定模拟试验。采用不同的Pt值，获得破裂后的岩石样品，定义破裂率＝破裂的节理单元数目/总的节理单元数目，当破裂率得到稳定值，即Pt的改变不再对模拟结果有影响，则视为得到了合适的Pt时。再次变化不同的网格尺寸，验证Pt取值与网格尺寸的无关性。

岩样的破坏率达到稳定值时，认为切向接触刚度Pt是合理的，如下：

式中，η为岩样的破坏率，nfailed，ntotal分别为破坏节理单元数目和节理单元总数目。经过多次试算表明，Pt≥5E能够满足要求。

步骤7.采用图5的三轴压缩、巴西劈裂和直接剪切模型，验证步骤2-步骤6标定的微观输入参数的可靠性，即所有的输入参数可取为表1中的值。

表1 FDEM所有输入参数来源依据

因此，根据上述步骤，本发明提供了如图6所示的FDEM输入参数标定流程。根据本发明提出的标定流程，所有的宏观输入参数可取为室内试验值，避免了二次标定，减少了被标参数的数目，提高了参数标定效率，实现了所有参数取值来源有据可依。

以上实施例仅是对本发明技术方案所做的举例说明。本发明所涉及的有限元-离散元耦合数值模拟程序(FDEM)输入参数标定方法并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。