一种模拟互层岩体数值模型的仿真方法、系统及处理终端与流程

1.本发明属于互层岩体仿真技术领域，尤其涉及一种模拟互层岩体数值模型的 仿真方法、系统及处理终端。


背景技术：

2.目前：软硬互层状岩体较多存在于隧道、地下洞室等领域，容易造成隧道塌 方、隧道掘进效率低等工程问题。相较于均质岩体软硬互层岩体破坏特征差异显 著，破坏过程较为复杂，更难预测其破坏先兆信息。岩层倾角，岩层层厚，软弱 层含量三类几何特征因子对互层岩体的破幻有较大影响。目前各类研究通过离散 元手段对影响因素作用下互层岩体力学参数变化规律及其细观层面展开深入分 析的研究较少，对多因素共同作用下互层岩体的细观破坏机制进行综合研究的成 果更少，造成现有结果对互层岩体的破坏机制认识不足以服务于实际工程，严重 制约了地下工程的发展。
3.随着我国轨道交通的发展，隧道穿越区域内大量的泥岩、砂岩组成的软硬互 层岩体分布较广，因其力学性质复杂，容易造成地表沉降、隧道塌方等工程稳定 性问题，然而目前关于软硬互层岩体力学特性及破坏机制不够全面、深入。现有 技术对于互层岩体力学特性的研究主要包括室内(物理)试验、数值模拟，室内 试验。室内试验制作试验过程比较复杂，耗时长，代价高，且在定义不同层厚比、 层厚、岩层倾角时会出现较大误差，在试样的养护过程中，养护条件、养护过程 等会出现较大误差，温度，水分等的因素不好控制，在实验时，人为误差，试验 仪器误差等都是较为复杂的影响试验结果的因素。数值模拟主要包括有限元、有 限差分法等，这些方法在模拟过程中，有限元模拟差值是基于网格的，需要人为 做好单元，消耗时间，需要考虑尺寸的敏感性，而离散元是通过单元体间的力学 行为叠加生成宏观响应，不需要对模型预设复杂的本构关系。不需要考虑试验过 程中的试样养护问题、人为误差较小，试验简单方便。能够较为准确的模拟岩石 在加载过程中的局部各向异性行为，较好的模拟裂隙时空演化过程和最终断裂。 对数值模拟有制作试样误差小(如岩层倾角、层厚比、层厚)过程较为简单，主 要是通过代码实现，岩体都是通过代码中的本构关系实现的 通过上述分析，现有技术存在的问题及缺陷为：
4.(1)现有技术不能准确获取互层岩体的相关参数。
5.(2)工程中互层岩体的几何特征因子(指层厚、层厚比、岩层倾角)往往 较复杂，没有有效的方法获取其力学特性。
6.(3)目前对于互层岩体力学特性的研究主要包括室内(物理)试验、数值 模拟，室内试验，室内试验制作试验过程比较复杂，耗时长，代价高，且在定义 不同层厚比、层厚、岩层倾角时会出现较大误差，在试样的养护过程中，养护条 件养护过程等会出现较大误差，在温度，水分等的因素不好控制，在实验时，人 为误差，试验仪器误差等都是较为复杂的影响试验结果的因素。
7.(4)数值模拟主要包括有限元、有限差分法等，这些方法在模拟过程中， 有限元模拟差值是基于网格的，需要人为做好单元，消耗时间，需要考虑尺寸的 敏感性。
8.解决以上问题及缺陷的难度为：
9.离散元是通过单元体间的力学行为叠加生成宏观响应，不需要对模型预设复 杂的本构关系。不需要考虑试验过程中的试样养护问题、人为误差较小，试验简 单方便。能够较为准确的模拟岩石在加载过程中的局部各向异性行为，较好的模 拟裂隙时空演化过程和最终断裂。对数值模拟有制作试样误差小(如岩层倾角、 层厚比、层厚)过程较为简单，主要是通过代码实现，岩体都是通过代码中的本 构关系实现的。
10.解决以上问题及缺陷的意义为：
11.随着我国轨道交通的发展，隧道穿越区域内大量的泥岩、砂岩组成的软硬互 层岩体分布较广，因其力学性质复杂，容易造成地表沉降、隧道塌方等工程稳定 性问题，然而目前关于软硬互层岩体力学特性及破坏机制不够全面、深入。基于 此，为完善互层岩体力学特性理论，揭示其破坏机制。本发明采用颗粒流数值模 拟从细观角度分析了在岩体倾角、层厚比、层厚等因子作用下互层岩体力学特性 响应规律与细观破坏特征。


技术实现要素：

12.针对现有技术存在的问题，本发明提供了一种模拟互层岩体数值模型的仿真 方法、系统及处理终端。
13.本发明是这样实现的，一种模拟互层岩体数值模型的仿真方法，包括：
14.利用distribute命令生成球体，并设置密度；
15.迭代计算使颗粒达到平衡；
16.生成不同岩层倾角、层厚比、层厚的互层岩体数值模型。
17.进一步，所述仿真方法进一步包括：
18.利用ball distribute命令生成球体，同时设置球体孔隙率、半径、球体生 成范围，并设置球体属性密度、阻尼；如：ball distribute porosity 0.01radius 3.0e
‑
4 6.0e
‑
4box 0 0.05 0 0.1ball attribute density 2180damp 0.7， 表示孔隙比0.01半径范围3.0e
‑
4至6.0e
‑
4范围内，在盒子5
×
10cm范围内生 成，密度2180，阻尼0.7。ball distribute命令按孔隙比分配可重叠的球到一 个区域，在迭代计算颗粒运动达到设定的孔隙比时颗粒停止运动。
19.进一步，强风化泥岩层厚a与强风化砂岩层厚b的比值m，层厚比的比值， 强风化泥岩层厚a表示为：强风化砂岩层厚b表示为： 20.进一步，所述利用distribute命令生成球体，并设置密度前，需进行： 设置整个模型生成的空间范围；通过domain命令，如domain extent
‑
0.1 0.2
ꢀ‑
0.1 0.2condition destroy为10cm
×
20cm的范围。设置空间范围内实体的 力学参数，包括粘结模量、刚度比、阻尼；通过cmat命令设置如：cmat defaultmodel linear method deform emod 1.0e9 kratio 0.0cmat default propertydp_nratio 0.5，表示范围内实体的变形模量为1gpa，
刚度比为0，阻尼比为 0.5。
21.进一步，所述生成不同岩层倾角、层厚比、层厚的互层岩体数值模型后，需 进行：
22.对颗粒进行分组，使互层岩体数值模型的不同图层中充满颗粒；
23.删除左右两侧墙体；
24.对不同组的颗粒赋予不同的平行粘结参数包括：平行粘结模量、刚度比、平 行粘结刚度比、颗粒摩擦系数、法向粘结强度、切向粘结强度、平行粘结摩擦角 等参数，用以模拟不同层的岩体属性。
25.本步骤的实现通过软件自带的parallel_bonded文件，对pb_deform emod、 krat等参数进行分组，对应不同图层的颗粒，并对pb_deform emod、krat等参 数赋值。
26.同时在openfoamt定义网格，确定网格坐标划分网格数量及比例。设置压强， 入口压强设置为zerogradient，出口压强设置为固定值，设置流体速度，入口速 度y方向为1，其余方向为0，出口速度设置为zerogradient。同时设置流体粘 度。
27.通过离散元软件读取网格数据节点及边界条件流体压强、速度。并设置运行 时间间隔，使数据在离散元中运行一定时间后，把数据交换到openfoam软件， 在openfoam中运行一定时间把数据交换到离散元，以此类推，不断进行交换， 使互层岩体与流体耦合。
28.进一步，所述对不同组的颗粒赋予不同的平行粘结参数后，还需进行：
29.打开应力、应变、裂纹记录及能量开关，记录裂纹数量及发展趋势、摩擦能、 应变能、动能变化趋势，记录应力、应变变化趋势；
30.裂纹数量及趋势的方法：通过离散元自带的fracture功能，不断对颗粒间 的断裂进行记录和更新，并通过history命令记录裂纹的数量及趋势。
31.摩擦能、应变能、动能变化趋势获取方法：通过set energy on函数调用， history记录个能量的变化趋势。
32.上下墙体以一定的速度分别向下、向上移动，应力为先上升后下降的趋势， 当下降应力值小于峰值应力的70％时(实现方法：通过math.abs函数记录试样的 实时应力，把应力下降那一刻的应力记为峰值应力，当实时应力小于峰值应力乘 以0.7时，试样停止加压)；墙体停止移动，试验结束，并测量单轴压缩试验的 弹性模量、及泊松比。
33.本发明的另一目的在于提供一种模拟互层岩体数值模型的仿真系统，包括：
34.互层岩体数值模型参数生成模块，用于设置整个模型生成的空间范围，还用 于设置空间范围内实体的力学参数，包括粘结模量、刚度比、阻尼；还用于墙体 的生成；
35.互层岩体数值模型生成模块，利用distribute命令生成球体，并设置密度； 迭代计算使颗粒达到平衡；生成不同岩层倾角、层厚比、层厚的互层岩体数值模 型；
36.不同层的岩体属性模拟模块，用于对颗粒进行分组，使互层岩体数值模型的 不同图层中充满颗粒；删除左右两侧墙体；对不同组的颗粒赋予不同的平行粘结 参数包括：平行粘结模量、刚度比、平行粘结刚度比、颗粒摩擦系数、法向粘结 强度、切向粘结强度、平行粘结摩擦角等参数，用以模拟不同层的岩体属性；
37.测量模块，用于记录裂纹数量及发展趋势、摩擦能、应变能、动能变化趋势， 记录应力、应变变化趋势；还用于上下墙体以一定的速度分别向下、向上移动， 应力为先上升后下降的趋势，当下降应力值小于峰值应力的70％时；墙体停止移 动，并测量单轴压缩试验的弹性模量、及泊松比。
38.本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和 处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时， 使得所述处理器执行所述方法。
39.本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序， 所述计算机程序被处理器执行时，使得所述处理器执行所述方法。
40.本发明的另一目的在于提供一种信息数据处理终端，所述信息数据处理终端 用于实现所述方法。
41.结合上述的所有技术方案，本发明所具备的优点及积极效果为：
42.本发明互层岩体的定义及仿真方法在软件中即可实现。本发明采用颗粒流数 值模拟从细观角度分析了在岩体倾角、层厚比、层厚等因子作用下互层岩体力学 特性响应规律与细观破坏特征。
附图说明
43.为了更清楚地说明本技术实施例的技术方案，下面将对本技术实施例中所需 要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本技术的一 些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以 根据这些附图获得其他的附图。
44.图1是本发明实施例提供的模拟互层岩体数值模型的仿真方法流程图。
45.图2是本发明实施例提供的互层岩体实物图。
46.图3是本发明实施例提供的裂纹图岩层倾角15
°
、层厚16.6mm、层厚 比为1:1仿真效果图。
47.图4是本发明实施例提供的(裂纹发展曲线岩层倾角15
°
、层厚16.6mm、 层厚比为1:1仿真效果图。
48.图5是本发明实施例提供的(能量发展曲线岩层倾角15
°
、层厚16.6mm、 层厚比为1:1仿真效果图。
49.图6是本发明实施例提供的(能量发展曲线岩层倾角15
°
、层厚 16.6mm、层厚比为1:1仿真效果图。
50.图7是本发明实施例提供的(岩体组构，岩层倾角15
°
、层厚16.6mm、 层厚比为1:1仿真效果图。
具体实施方式
51.为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对 本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释 本发明，并不用于限定本发明。
52.针对现有技术存在的问题，本发明提供了一种模拟互层岩体数值模型的仿真 方法，下面结合附图对本发明作详细的描述。
53.如图1所示，本发明提供的模拟互层岩体数值模型的仿真方法包括：
54.s101，设置整个模型生成的空间范围，10cmx20cm。
55.s102，设置空间范围内实体的力学参数，包括粘结模量、刚度比、阻尼。
56.s103，墙体的生成5cmx10cm，并设置墙体id号，及墙体生成位置(此时墙 体是固定的)。
57.s104，利用ball distribute命令生成球体，并设置球体的半径大小、孔隙率、 生成范围，并对求球体赋予属性包括密度、阻尼。
58.s105，迭代计算使颗粒达到孔隙率，试样平衡。
59.s106，导入可以生成不同岩层倾角、层厚比、层厚的cad图。
60.s107，对颗粒进行分组，使cad的不同图层中充满颗粒。
61.s108，删除左右两侧墙体。
62.s109，对不同组的颗粒赋予不同的平行粘结参数包括：平行粘结模量、刚度 比、平行粘结刚度比、颗粒摩擦系数、法向粘结强度、切向粘结强度、平行粘结 摩擦角等参数，用以模拟不同层的岩体属性。
63.s110，打开应力、应变、裂纹记录及能量开关，记录裂纹数量及发展趋势、 摩擦能、应变能、动能变化趋势，记录应力、应变变化趋势(均为软件自带函数)。
64.本步骤的应力计算方法:利用pfc内嵌的fish语言记录上下加压板上接触点 垂直力的总和除以试样的宽度，得到轴向应力。试样宽度的获得：遍历每个球体 的信息包括位置、半径等信息，通过x方向位置坐标最小的球体减去其半径得到 x方向最小值，通过x方向位置坐标最大的球体加上其球体半径得到x方向最大 值，如loop foreach bp ball.list；
65.local ball_xmin＝ball.pos.x(bp)
‑
ball.radius(bp)
66.xmin＝math.min(xmin,ball_xmin)
67.local ball_xmax＝ball.pos.x(bp) ball.radius(bp)
68.xmax＝math.max(xmax,ball_xmax)
69.end_loop
70.def axial_stress_wall
71.local force1＝
‑
wall.force.contact(wp_top,vertical_direction)
72.local force2＝wall.force.contact(wp_bottom,vertical_direction)
73.axial_stress_wall＝0.5*(force1 force2)/cross_sectional_area
74.end然后通过x方向最大值减去最小值得到试样的宽度如local diameter_＝ xmax
‑
xmin。轴向应力的获得：通过wall.force.contact内置函数获得作用在 上下加圧板上的所有接触力的总和，并除以2获得平均轴向应力。
75.应变计算方法：通过wall.disp函数获得上下加压板的实时位移，通过 wall.find获取上下加压板的信息，并由上部加圧板的位置减去下部加圧板的位 置获取试样初始高度。实时位移除以试样初始高度获得应变，如global wp_top＝ wall.find(1)；assume wall 1is the top wall
76.global wp_bottom＝wall.find(2)；assume wall 2is the bottom wall global vertical_direction＝global.dim
77.global sample_height＝wall.pos(wp_top,vertical_direction)
‑ꢀ
wall.pos(wp_bottom,vertical_direction)
78.def axial_strain_wall
79.；
80.；compute axial strain(positive tension)using walls
81.；
82.；assumes global variable sample_width has been set
83.；
84.axial_strain_wall＝ 2.0*wall.disp(wp_top,vertical_direction)/sample_height
85.end。
86.s111，上下墙体以一定的速度分别向下、向上移动，应力为先上升后下降的 趋势，当下降应力值小于峰值应力的70％时。墙体停止移动，试验结束，并测量 单轴压缩试验的弹性模量、及泊松比。
87.s112，对破坏试样进行组构分析，具体实现方法为：首先把圆分为72份，遍 历球与球之间的组构数，获得破坏岩体的组构分布
88.在本发明中，关于互层岩体的定义如下：
89.互层岩体为两种或两种以上属性的岩体组成的岩体。
90.步骤s109中，同时在openfoamt定义网格，确定网格坐标划分网格数量及 比例。设置压强，入口压强设置为zerogradient，出口压强设置为固定值，设置 流体速度，入口速度y方向为1，其余方向为0，出口速度设置为zerogradient。 同时设置流体粘度。
91.通过离散元软件读取网格数据节点及边界条件流体压强、速度。并设置运行 时间间隔，使数据在离散元中运行一定时间后，把数据交换到openfoam软件， 在openfoam中运行一定时间把数据交换到离散元，以此类推，不断进行交换， 使互层岩体与流体耦合。
92.如图2所述岩层倾角15
°
、层厚16.6mm、层厚比为1:1(设定1、3、5为强 风化泥岩，2、4、6为强风化砂岩)。
93.岩层倾角：岩层与水平面的夹角。
94.层厚：每一层的厚度，设定试样高度h为10cm，层数n为6层层厚为16.6mm， 一共4层层厚为25mm。
95.层厚比：强风化泥岩层厚a与强风化砂岩层厚b的比值m，如层厚比为0.5， 则为的比值，那么强风化泥岩层厚a可以表示为：强风化砂 岩层厚b表示为：
96.本发明提供的方法进行了试验验证。具体结果如表1
‑
3所示。
97.固定层厚为16.6mm、岩层倾角为0
°
，层厚比对互层岩体力学特性影响如下 表(图中0代表纯砂岩，∞代表纯泥岩)结果如表1
[0098][0099]
层厚比为0.8、层厚为16.6mm的互层岩体，岩层倾角对互层岩体力学特性的 无影响结果如表2：
[0100][0101]
层厚比为0.8、岩层倾角为50
°
、不同层厚对互层岩体力学特性的影响结果如
[0102]
表3：
[0103][0104]
在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上；术 语“上”、“下”、“左”、“右”、“内”、“外”、“前端”、“后端”、
ꢀ“
头部”、“尾部”等指示的方位或位置关系为基于附图所示的方位或位置关系， 仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须 具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。 此外，术语“第一”、“第二”、“第三”等仅用于描述目的，而不能理解为指 示或暗示相对重要性。
[0105]
应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来 实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适 当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术 人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器 控制代码中来实现，例如在诸如磁盘、cd或dvd
‑
rom的载体介质、诸如只读 存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供 了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、 诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设 备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软 件实现，也可以由上述硬件电路和软件的结合例如固件来实现。
[0106]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此， 任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神 和原则之
内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围 之内。