(19)中华人民共和国国家知识产权局
(12)发明专利申请
(10)申请公布号 (43)申请公布日 (21)申请号 202011186326.X
(22)申请日 2020.10.30
(71)申请人 武汉大学
地址 430072 湖北省武汉市武昌区珞珈山
武汉大学
(72)发明人 邓鹏海 刘泉声 吴坚
(74)专利代理机构 武汉科皓知识产权代理事务
所(特殊普通合伙) 42222
代理人 李炜
(51)Int.Cl.
G01N 3/62(2006.01)
(54)发明名称
(57)摘要
本发明公开了一种有限元‑离散元耦合数值
模拟程序(FDEM)输入参数快速标定方法，步骤如
下：(1)采用单轴压缩、巴西劈裂和三轴压缩室内
压缩模拟试验峰前阶段曲线标定节理罚值P f ；
(3)基于算法改进和单个节理单元压剪破坏标定
法向接触刚度P n ；(4)采用单轴压缩模拟试验全
过程曲线标定II型断裂能G II ；(5)采用直接拉伸
试验标定I型断裂能G I ；(6)再次采用单轴压缩模
拟试验标定切向接触刚度P t ；(7)采用三轴压缩、
巴西劈裂和直接剪切模拟试验验算所标定参数
的可靠性。本发明提出的标定方法效率高，消除
了参数取值对网格尺寸的依赖性，能够适用于工
程尺度的模拟。权利要求书2页 说明书6页 附图4页CN 112362520 A 2021.02.12
C N 112362520
A
1.一种有限元-离散元耦合数值模拟程序输入参数快速标定方法，其特征在于，包括以下步骤：
(1)基于室内岩石力学实验，获取岩石宏观物理力学参数，包括单轴抗压强度σc，弹性模量E，泊松比ν，抗拉强度σt，粘聚力c，内摩擦角和岩样的破坏模式；
(2)采用不同的节理罚值P f，进行标准岩样单轴压缩模拟试验，监测峰前应力-应变曲线，得到数值模拟的弹性模量值并与输入值对比，通过不断试算得到节理罚值P f的取值范围；
(3)采用两个三角形单元的压缩-剪切破坏模拟试验，标定得到法向接触刚度P n取值；
(4)采用单轴压缩模拟试验，设置抗拉强度σt为较大值，以使岩样仅发生剪切破坏，获得II型断裂能G II；
(5)采用标准岩样直接拉伸试验，获得I型断裂能G I；
(6)根据步骤2-5得到的节理罚值P f、法向接触刚度P n、II型断裂能G II和I型断裂能G I，再次采用单轴压缩模拟试验，变化不同的切向接触刚度P t，监测岩样的破坏率(破坏率被定义为破坏节理单元数目与总节理单元数目的比值)，直至获得稳定的破坏率，获得切向接触刚度P t的取值范围；
(7)根据步骤1-6得到的宏观输入参数和微观输入参数，进行三轴压缩、巴西劈裂和直接剪切模拟试验，验证步骤(2)-步骤(6)所标参数的可靠性。
2.根据权利要求1所述的方法，其特征在于：所述步骤(1)中岩石力学实验包括单轴压缩、直接拉伸、巴西劈裂和三轴压缩试验。
3.根据权利要求1所述的方法，其特征在于：所述步骤(2)的输入值为室内试验值；步骤(2)中试算方法如下：采用不同的节理罚值P f，进行单轴压缩模拟试验，得到对应的峰前弹性模量E out，直至该比值与输入弹性模量E in的比值接近于1为止，得到合理的P f取值范围；变化不同的网格尺寸h和输入弹性模量E in，验证P f取值范围的鲁棒性，得到通用的P f值取值范围。
4.根据权利要求1所述的方法，其特征在于：所述步骤(3)中，两个三角形单元之间能够由粘接状态平稳过渡到接触状态。
5.根据权利要求1所述的方法，其特征在于：所述步骤(3)中法向接触刚度P n计算公式，如下：
式中，P n(i-j)为三角形单元i与三角形单元j间的法向接触刚度，h i、h j分别为三角形单元i、j的高，l i、l j分别为三角形单元i、j的长，P b为基本刚度，采用下式计算：
P b＝αP f
式中，α为系数。
6.根据权利要求5所述的方法，其特征在于：所述步骤(3)中，压缩-剪切破坏模拟试验中，通过改变不同的α值，直至压剪破坏瞬间所有的节点力平稳地由粘接状态平稳过渡到接触状态，消除节理单元破坏瞬间节点力的异常震荡，以获得最佳α取值。
7.根据权利要求1所述的方法，其特征在于：所述步骤(4)中II型断裂能G II的获得方法
为：单轴压缩模拟试验后，宏观输入参数取为试验值，节理罚值P f和法向接触刚度P n分别根据步骤2和步骤3的方法确定，切向接触刚度P t设为较大值(如P t＝10E)，变化II型断裂能
G II，直至获取与室内试验相吻合的模拟结果。
8.根据权利要求1所述的方法，其特征在于：所述步骤(4)中，修正了FDEM中现有的节理单元本构模型，将现有的应力-位移本构模型修正为应力-应变本构模型，消除了参数取值对网格尺寸的依赖性，如下：
式中，o r，s r分别为节理单元的极限拉伸和剪切位移，f s为节理单元抗剪强度。通过上式表明，节理单元的破裂由应变控制，而不再是位移控制，消除了G I和G II取值对网格尺寸的依赖性。
9.根据权利要求1所述的方法，其特征在于：所述步骤(6)中，岩样的破坏率达到稳定值时，认为切向接触刚度P t是合理的，如下：
式中，η为岩样的破坏率，n failed，n total分别为破坏节理单元数目和节理单元总数目。
10.根据权利要1所述的方法，其特征在于：综合采用了三轴压缩、巴西劈裂和直接剪切模拟试验验证所标参数的可靠性，因为上述三种试验的模拟结果可直接与理论结果对比，可便捷、准确判定模拟结果的正确性。
一种有限元-离散元耦合数值模拟程序(FDEM)输入参数快速
标定方法
技术领域
[0001]本发明属于岩石力学与岩石工程领域，具体涉及有限元-离散元耦合数值模拟程序(FDEM)的输入参数标定方法及被标参数的可靠性验证方法。
背景技术
[0002]岩石的破裂问题是研究者及工程师需要重点关注的问题，因为岩石破裂后会直接影响到如隧道围岩及岩石边坡等的稳定性。岩石的破裂经历了从弹性变形到塑性变形，再到岩石破裂后块体的相互接触及块体大运动全过程，一方面岩石材料复杂，天然岩体含有大量微裂隙，且实际工程岩体往往以复合岩体的形式出现；另一方面，岩石的力学性质与其赋存环境、应力历史、应力路径和研究尺度等密切相关；再者，对岩石破裂的研究不仅仅是为了研究其破裂过程及影响因素，更重要的是研究如何促进岩石
破裂(如矿石开采)和抑制岩石破裂以便控制其稳定性(如隧道围岩)。
[0003]面对如此复杂的材料性质、赋存环境和实际需求，理论解析显得力不从心，而无论室内试验或工程现场试验，只对特定工程有效、不具有预测性，且研究成本高昂、研究周期漫长。为了便于岩石破裂的研究，数值模拟方法应运而生，大致地，可分为连续性方法(如有限元法FEM，有限差分法FDM和边界元法BEM等)、非连续性方法(如颗粒元法PFC，通用离散单元法UDEC和不连续变形法DDA等)和二者的耦合方法(如数值流形法NMM及本说明中的FDEM 法等)。其中，FDEM由于其计算效率高、能够重现岩石弹性-塑性-破裂全过程等优点得到了广泛应用，包括室内试验模拟和实际工程的模拟研究。
[0004]准确的输入参数是保证数值模拟结果可靠性的前提，遗憾的是，目前尚无高效便捷的FDEM输入参数标定方法。目前普遍采用单轴压缩和巴西劈裂模拟试验不断试算不同的参数组合，以期得到与室内试验结果相吻合的输入参数，需要重复高达上千次的试算方能得到较为理想的输入参数组合。更为严重的是，该输入参数具有严重的网格尺寸依赖性，通过室内试验标定得到的参数直接应用于工程尺度的模拟是存疑的。因此，本发明旨在提出一种快速便捷的FDEM输入参数标定方法，消除参数取值对网格尺寸的依赖性，并提出被标参数的可靠性验证方法。
发明内容
[0005]为了解决上述技术问题，本发明提供一种快速的FDEM输入参数标定方法，并消除参数取值对网
格尺寸的依赖性，同时提出被标参数可靠性验证方法。
[0006]本发明提供的技术方案如下：
[0007]一种有限元-离散元耦合数值模拟程序输入参数快速标定方法，包括以下步骤：[0008](1)基于室内岩石力学实验，获取岩石宏观物理力学参数，包括单轴抗压强度σc，弹性模量E，泊松比ν，抗拉强度σt，粘聚力c，内摩擦角和岩样的破坏模式；
[0009](2)采用不同的节理罚值P f，进行标准岩样单轴压缩模拟试验，监测峰前应力-应
变曲线，得到数值模拟的弹性模量值并与输入值对比，通过不断试算得到节理罚值P f的取值范围；
[0010](3)采用两个三角形单元的压缩-剪切破坏模拟试验，标定得到法向接触刚度P n取值；
[0011](4)采用单轴压缩模拟试验，设置抗拉强度σt为较大值，以使岩样仅发生剪切破坏而不发生拉伸破坏，获得II型断裂能G II；
[0012](5)采用标准岩样直接拉伸试验，获得I型断裂能G I；
[0013](6)根据步骤2-5得到的节理罚值P f、法向接触刚度P n、II型断裂能G II和I型断裂能G I，再次采用
单轴压缩模拟试验，变化不同的切向接触刚度P t，监测岩样的破坏率(破坏率被定义为破坏节理单元数目与总节理单元数目的比值)，直至获得稳定的破坏率，获得切向接触刚度P t的取值范围；
[0014](7)根据步骤1-6得到的宏观输入参数和微观输入参数，进行三轴压缩、巴西劈裂和直接剪切模拟试验，验证步骤(2)-步骤(6)所标参数的可靠性。
[0015]进一步，所述步骤(1)中岩石力学实验包括单轴压缩、直接拉伸、巴西劈裂和三轴压缩试验。
[0016]所述步骤(2)的输入值为室内试验值。
[0017]进一步，所述步骤(2)中试算方法如下：采用不同的节理罚值P f，进行单轴压缩模拟试验，得到对应的峰前弹性模量E out，直至该比值与输入弹性模量E in的比值接近于1为止，得到合理的P f取值范围；变化不同的网格尺寸h和输入弹性模量E in，验证P f取值范围的鲁棒性，得到通用的P f值取值范围。
[0018]进一步，所述步骤(3)中，两个三角形单元之间能够由粘接状态平稳过渡到接触状态。
[0019]进一步，所述步骤(3)中法向接触刚度P n计算公式，如下：
[0020]
[0021]式中，P n(i-j)为三角形单元i与三角形单元j间的法向接触刚度，h i、h j分别为三角形单元i、j的高，l i、l j分别为三角形单元i、j的长，P b为基本刚度，采用下式计算：
[0022]P b＝αP f
[0023]式中，α为系数。
[0024]进一步，所述步骤(4)中II型断裂能G II的获得方法为：单轴压缩模拟试验后，宏观输入参数取为试验值(弹性模量E，泊松比ν，粘聚力c，内摩擦角)，节理罚值P f和法向接触刚度P n分别根据步骤2和步骤3的方法确定，切向接触刚度P t设为较大值(如P t＝10E)，变化II型断裂能G II，直至获取与室内试验相吻合的模拟结果。
[0025]进一步，所述步骤(4)中，修正了FDEM中现有的节理单元本构模型，将现有的应力-位移本构模型修正为应力-应变本构模型，消除了参数取值对网格尺寸的依赖性，如下：
[0026]
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议