连续-非连续数值模型连接方法与流程

1.本技术涉及岩土工程技术领域，特别是一种连续-非连续数值模型连接方法。


背景技术：

2.在岩土工程领域，主要采用数值模拟方法研究地震波的传播及动力作用下的破坏过程。连续数值模拟方法如有限单元法、有限差分法具有明确的力学意义，应用最为广泛，其波场连续、效果最好，但在模拟岩土破坏方面效果较差；而颗粒离散单元法等非连续数值模拟方法虽然可以模拟大变形引起的破坏，但其位移场、应力场不连续。
3.鉴于连续数值模拟方法与非连续数值模拟方法各有优缺点，因此将非大变形区域处理为连续数值模型、大变形破坏区域处理为颗粒离散元等非连续模型成为岩土地震等动力破坏研究的重要手段。
4.当采用连续-非连续数值耦合方法开展数值模拟时，最常规的方法是在连续部分与非连续部分设置墙，令离散的球与墙作用，而构成墙的顶点与连续部分的单元节点相匹配一致运动，实现连续-非连续数值模型力-位移的传递。
5.但是在连续-非连续数值模型之间起中介作用的墙不满足牛顿第二定律，因此应力波在连续-非连续数值模型的界面上无法反映正常的应力波透射条件。导致采用现有模型计算后所得结果不符合地震动力载荷和地震力学机理。


技术实现要素：

6.本技术针对上述技术问题提供了一种连续-非连续数值模型连接方法，改善地震作用下应力波的模拟效果，提高模拟的真实性与参考价值。
7.本技术提供了一种连续-非连续数值模型连接方法，包括以下步骤：
8.步骤s1：将建好的连续数值模型导入flac3d6.0平台，遍历节点找出模型范围，并将模型左、右、前、后、底边界设置为粘滞边界以吸收能量；
9.步骤s2：设置非连续数值模型对应的分组，将分组边界放大以保证放大后范围比实际区域超出至少2个平均颗粒尺寸为准；
10.步骤s3：在模型区域内生成特定颗粒尺寸的圆球，采用颗粒伺服法均匀化处理颗粒体系，使颗粒体系趋于平衡；
11.步骤s3包括以下步骤：
12.步骤s31：采用步骤s2定义的非连续区域边界几何定义为几何组，名称可自指定，此处以“exterior”称之，同时将其转化为墙，组名也为“exterior”；
13.步骤s32：利用flac3d平台中随机投放圆球的功能，在步骤s2生成的几何组“exterior”控制下生成一系列圆球(半径r)，r值处于圆球最小半径rmin最大半径rmax之间，呈均匀随机分布，所有颗粒半径的平均值约等于最小与最大半径的平均值，记作rave，其中最大颗粒半径与最小颗粒半径比值以1.0～3.0之间为宜；
14.步骤s33：将模型区域内的所有圆球之间的接触设置为线性接触模型，设置法向刚
度参数kn＝1e8n/m，切向刚度参数ks＝1e8n/m，令颗粒体系以时间步长1.0运行一定迭代步，如2万步，令颗粒体系基本平衡；
15.步骤s34：遍历所得几何组“exterior”的顶点，计算出容纳几何范围的最小长方体，由xmin0，xmax0，ymin0，ymax0，zmin0，zmax0限定，在该长方体范围内，设定测量圆半径取平均颗粒半径的5倍，生成一系列规则的测量圆；
16.步骤s35：通过测量圆控制应力大小，通过缩小、放大所有颗粒半径，令测量圆的平均应力值为100kpa使得颗粒体系趋于平衡，平衡后模型即为非连续模型部分控制域内的非连续颗粒体系；
17.步骤s4：在模型内对与实体单元重叠的重叠圆球，采用单元内插法计算重叠圆球的位移、速度、加速度，通过重叠圆球颗粒之间的传力将荷载传递到离散模型范围内。
18.步骤s5：标定搭接区域细观力学参数，按照牛顿第二定律将离散模型和连续数值模型中的状态相互传递，逐步计算地震波传播过程，记录节点和离散球的位移、速度、加速度峰值，最终得反应谱数据。
19.步骤s5包括以下步骤：
20.步骤s51：对颗粒体系细观参数标定时，应满足：模型搭接区颗粒参数体现出的宏观性质与连续模型宏观性质一致或接近一致，以保证模型搭接区应力波的稳定；
21.步骤s52：设置颗粒体系计算时间步为一较小值dt，如1.0e-7(足够小以保证应力波稳定),在已知任一时刻t模型各节点的速度时，通过插值法计算重叠区域内颗粒所处位置的速度，将所得速度赋值给处于相应位置的球i，从而实现应力波产生的波动传递到不连续区域的颗粒体系中；
22.步骤s53：在t+dt时刻颗粒体系所受的力传递到球i后，该力作为集中力传递到连续数值模型中参与计算，重复步骤s52～s53，直至应力波在该边界上连续光滑的传递，从而实现应力波的透射；
23.步骤s54：在每一个时间步上均对比节点和离散球的位移、速度、加速度峰值，并记录所得最大值，得到反应谱数据。
24.优选地，步骤s1包括以下步骤：
25.步骤s11：将非连续域作为连续模型处理，将建立的连续模型导入flac3d6.0及以上版本平台。模型由各节点坐标和单元构成，其中节点坐标为(x，y，z)，单元为(每个单元含8个节点以及单元分组名称)；
26.步骤s12：原则上模型在水平面上投影为规则矩形，故设定模型范围在x向最小值为x
min
(初值100000)，最大值为x
max
(初值为-10000)；模型范围在：y向最小值为y
min
(初值100000)，最大值为ymax(初值为-100000)；模型在z向最小值z
min
(初值100000)，遍历模型所有节点并更新x
min
、x
max
、y
min
、y
max
、z
min
。
27.优选地，更新步骤为：若节点x的坐标小于x
min
，则用x替换x
min
；若节点坐标x值大于x
max
，则用x替换x
max
；若节点坐标y值小于y
min
，则用y替换y
min
；若节点坐标y大于y
max
，则用y替换y
max
；若节点坐标z小于z
min
，则用z替换z
min
。
28.优选地，步骤s12还包括：设置边界容差error＝0.05，设置模型左边界x坐标处于[xmin-error，xmin+error]，右边界x坐标处于[xmax-error，xmax+error]，前边界y坐标处于[ymin-error，ymin+error]，y坐标处于[ymax-error，ymax+error],底边界z坐标处于
[zmin-error，zmin+error]；并在模型的法线方向和水平方向上分别设置独立的阻尼器以吸收模型内部振动产生的入射波；
[0029]
阻尼器提供的法向粘滞力tn和切向粘滞力ts如式(1)所示：
[0030]
tn＝-ρc
p
vn[0031]
ts＝-ρc
svs
式(1)
[0032]
式中：vn,vs分别为模型边界上法向和切向的速度分量，ρ为介质密度，c
p
,cs分别为模型的p波和s波的波速。
[0033]
优选地，步骤s2中的分组边界放大操作包括以下步骤：
[0034]
步骤s21：针对存在破坏风险的连续模型分组，将其外边界几何输出为非连续模型边界，该边界由一系列三角形集合成，估计连续模型平均颗粒半径rave0的值；
[0035]
步骤s22：对非连续模型边界几何进行法向判断，判断条件及对应操作：
[0036]
1)如果模型边界的法向量为正，则将边界下半部分的位置下移至少2倍平均颗粒半径；
[0037]
2)如果模型上边界为边坡面，则上边界不可更改；
[0038]
3)如果模型上边界不是模型外表面而是中间层，上边界也可向上平移至少2倍平均颗粒半径。
[0039]
优选地，步骤s4包括以下步骤：
[0040]
步骤s41：遍历所有实体单元zone和颗粒ball，根据实体单元和颗粒的几何位置将模型分为：连续模型区、离散模型区、模型搭接区，模型中岩桥区的颗粒即为边界控制颗粒，该区域为实体单元和颗粒重叠的重叠区域；
[0041]
步骤s42：由于一个实体网格单元可能与一个或者多个颗粒重叠，重叠区域位移的传递按下式计算：
[0042][0043]
式中：在离散区域中α等于1，β等于0；在连续区域中，α等于0，β等于1；对于ball-zone耦合方法，α和β在重叠区域中从0到1线性变化，αj和βi分别是离散颗粒和连续单元的系数，取值范围为0.01到0.99；mj为离散颗粒的质量，dj为离散颗粒的位移，f
jtot
是作用在离散颗粒上的外力，λj为离散颗粒的拉格朗日乘数，nj为离散颗粒的数目；mi为连续单元的质量，ui为连续单元的位移，f
itot
是作用在连续单元上的外力，λk为连续单元节点上的拉格朗日乘数，ni为连续单元的数目；k为运动矩阵，k
jk
是由经典插值函数定义的单元运动矩阵，uk是包围离散颗粒j的连续单元i的八个节点(k)的位移函数。
[0044]
本技术能产生的有益效果包括：
[0045]
1)本技术所提供的连续-非连续数值模型连接方法，该方法能有效实现连续数值模型中应力波圆滑地传递到非连续模型区域，使得地震动力荷载的传递更符合实际，得到的结果更符合力学机理。
[0046]
2)本技术所提供的连续-非连续数值模型连接方法，通过在不同时步内，反复利用
连续数值模型与离散数值模型连接处边界控制颗粒，使得不同区域运动和受力相互传递，实现连续数值模型与非连续数值模型之间的应力波按其传播特点进行传播，获得更加真实、可靠地模拟结果，提高模拟结果的真实有效性。
附图说明
[0047]
图1为本技术提供的连续-非连续数值模型连接方法流程示意图；
[0048]
图2是本技术实施例中连续-非连续耦合的边坡数值模型；
[0049]
图3是本技术实施例中单个单元节点及其退化形式；
[0050]
图4是本技术实施例中连续-非连续数值模型的边界；
[0051]
图5是本技术实施例中扩大的非连续数值模型边界；
[0052]
图6是本技术实施例中非连续颗粒状态监控设置的测量球；
[0053]
图7是本技术实施例中非连续数值模型平衡状态的监控曲线；
[0054]
图8是本技术实施例中边界搭接球位移插值处理；其中a)为各区搭接示意图；b)为单元内任一离散颗粒j所处位置示意图；
[0055]
图9是本技术实施例中耦合边界地震响应计算结果；
具体实施方式
[0056]
为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。通常在此处附图中描述和示出的本发明实施方式的组件可以以各种不同的配置来布置和设计。
[0057]
因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
[0058]
本实施方案中所用控制器为现有结构，且控制电路通过本领域的技术人员简单的编程即可实现，属于本领域的公知常识，仅对其进行使用，不进行改造，故不再详细描述控制方式和电路连接。
[0059]
本技术中未详述的且并不用于解决本技术技术问题的技术手段，均按本领域公知常识进行设置，且多种公知常识设置方式均可实现。
[0060]
参见图1，本技术提供的连续-非连续数值模型连接方法，包括以下步骤：
[0061]
步骤s1：将建好的连续数值模型导入flac3d6.0平台，遍历节点找出模型范围，并将模型左、右、前、后、底边界设置为粘滞边界以吸收能量；
[0062]
步骤s2：设置非连续数值模型对应的分组，将分组边界放大以保证放大后范围比实际区域超出至少2个平均颗粒尺寸为准；
[0063]
步骤s3：在模型区域内生成特定颗粒尺寸的圆球，采用颗粒伺服法均匀化处理颗粒体系，使颗粒体系趋于平衡；
[0064]
步骤s3包括以下步骤：
[0065]
步骤s31：采用步骤s2定义的非连续区域边界几何定义为几何组，名称可自指定，
此处以“exterior”称之，同时将其转化为墙，组名也为“exterior”；
[0066]
步骤s32：利用flac3d平台中随机投放圆球的功能，在步骤s2生成的几何组“exterior”控制下生成一系列圆球(半径r)，r值处于圆球最小半径rmin最大半径rmax之间，呈均匀随机分布，所有颗粒半径的平均值约等于最小与最大半径的平均值，记作rave，其中最大颗粒半径与最小颗粒半径比值以1.0～3.0之间为宜；
[0067]
步骤s33：将模型区域内的所有圆球之间的接触设置为线性接触模型，设置法向刚度参数kn＝1e8n/m，切向刚度参数ks＝1e8n/m，令颗粒体系以时间步长1.0运行一定迭代步，如2万步，令颗粒体系基本平衡；
[0068]
步骤s34：遍历所得几何组“exterior”的顶点，计算出容纳几何范围的最小长方体，由xmin0，xmax0，ymin0，ymax0，zmin0，zmax0限定，在该长方体范围内，设定测量圆半径取平均颗粒半径的5倍，生成一系列规则的测量圆；
[0069]
步骤s35：通过测量圆控制应力大小，通过缩小、放大所有颗粒半径，令测量圆的平均应力值为100kpa使得颗粒体系趋于平衡，平衡后模型即为非连续模型部分控制域内的非连续颗粒体系；
[0070]
步骤s4：在模型内对与实体单元重叠的重叠圆球，采用单元内插法计算重叠圆球的位移、速度、加速度，通过重叠圆球颗粒之间的传力将荷载传递到离散模型范围内。
[0071]
步骤s5：标定搭接区域细观力学参数，按照牛顿第二定律将离散模型和连续数值模型中的状态相互传递，逐步计算地震波传播过程，记录节点和离散球的位移、速度、加速度峰值，最终得反应谱数据。
[0072]
步骤s5包括以下步骤：
[0073]
步骤s51：对颗粒体系细观参数标定时，应满足：模型搭接区颗粒参数体现出的宏观性质与连续模型宏观性质一致或接近一致，以保证模型搭接区应力波的稳定；
[0074]
步骤s52：在已知任一时刻t模型各节点的速度时，通过插值法计算重叠区域内颗粒所处位置的速度，将所得速度赋值给处于相应位置的球i，从而实现应力波产生的波动传递到不连续区域的颗粒体系中；
[0075]
步骤s53：在t+dt时刻颗粒体系所受的力传递到球i后，该力作为集中力传递到连续数值模型中参与计算，重复步骤s52～s53，直至应力波在该边界上连续光滑的传递，从而实现应力波的透射；
[0076]
步骤s54：在每一个时间步上均对比节点和离散球的位移、速度、加速度峰值，并记录所得最大值，得到反应谱数据。
[0077]
采用本技术提供方法能通过在模型区域内生成颗粒尺寸r的圆球，并对所得圆球，并对圆球进行上述操作，从而以颗粒的形式实现通过重叠圆球颗粒之间的传力将荷载传递到离散模型范围内，从而有效改善模型的力传导模拟准确性。
[0078]
优选地，步骤s1包括以下步骤：
[0079]
步骤s11：将非连续域作为连续模型处理，将建立的连续模型导入flac3d6.0及以上版本平台。模型由各节点坐标和单元构成，其中节点坐标为(x，y，z)，单元为(每个单元含8个节点以及单元分组名称)；
[0080]
步骤s12：原则上模型在水平面上投影为规则矩形，故设定模型范围在x向最小值为x
min
(初值100000)，最大值为x
max
(初值为-10000)；模型范围在：y向最小值为y
min
(初值
100000)，最大值为ymax(初值为-100000)；模型在z向最小值z
min
(初值100000)，遍历模型所有节点并更新x
min
、x
max
、y
min
、ymax、z
min
。
[0081]
优选地，更新步骤为：若节点x的坐标小于x
min
，则用x替换x
min
；若节点坐标x值大于x
max
，则用x替换x
max
；若节点坐标y值小于y
min
，则用y替换y
min
；若节点坐标y大于y
max
，则用y替换y
max
；若节点坐标z小于z
min
，则用z替换z
min
。
[0082]
优选地，步骤s12还包括：设置边界容差error＝0.05，设置模型左边界x坐标处于[xmin-error，xmin+error]，右边界x坐标处于[xmax-error，xmax+error]，前边界y坐标处于[ymin-error，ymin+error]，y坐标处于[ymax-error，ymax+error],底边界z坐标处于[zmin-error，zmin+error]；并在模型的法线方向和水平方向上分别设置独立的阻尼器以吸收模型内部振动产生的入射波；
[0083]
阻尼器提供的法向粘滞力tn和切向粘滞力ts如式(1)所示：
[0084]
tn＝-ρc
p
vn[0085]
ts＝-ρc
svs
式(1)
[0086]
式中：vn,vs分别为模型边界上法向和切向的速度分量，ρ为介质密度，c
p
,cs分别为模型的p波和s波的波速。
[0087]
优选地，步骤s2中的分组边界放大操作包括以下步骤：
[0088]
步骤s21：针对存在破坏风险的连续模型分组，将其外边界几何输出为非连续模型边界，该边界由一系列三角形集合成，估计连续模型平均颗粒半径rave0的值；
[0089]
步骤s22：对非连续模型边界几何进行法向判断，判断条件及对应操作：
[0090]
1)如果模型边界的法向量为正，则将边界下半部分的位置下移至少2倍平均颗粒半径；
[0091]
2)如果模型上边界为边坡面，则上边界不可更改；
[0092]
3)如果模型上边界不是模型外表面而是中间层，上边界也可向上平移至少2倍平均颗粒半径。
[0093]
优选地，步骤s4包括以下步骤：
[0094]
步骤s41：遍历所有实体单元zone和颗粒ball，根据实体单元和颗粒的几何位置将模型分为：连续模型区、离散模型区、模型搭接区，模型中岩桥区的颗粒即为边界控制颗粒，该区域为实体单元和颗粒重叠的重叠区域；
[0095]
步骤s42：由于一个实体网格单元可能与一个或者多个颗粒重叠，重叠区域位移的传递按下式计算：
[0096][0097]
式中：在离散区域中α等于1，β等于0；在连续区域中，α等于0，β等于1；对于ball-zone耦合方法，α和β在重叠区域中从0到1线性变化，αj和βi分别是离散颗粒和连续单元的系数，取值范围为0.01到0.99；mj为离散颗粒的质量，dj为离散颗粒的位移，f
jtot
是作用在离散颗粒上的外力，λj为离散颗粒的拉格朗日乘数，nj为离散颗粒的数目；mi为连续单元的质量，
ui为连续单元的位移，f
itot
是作用在连续单元上的外力，λk为连续单元节点上的拉格朗日乘数，ni为连续单元的数目；k为运动矩阵，k
jk
是由经典插值函数定义的单元运动矩阵，uk是包围离散颗粒j的连续单元i的八个节点(k)的位移函数。
[0098]
实施例
[0099]
参见图2～9，某水电工程堆积体边坡，规模接近5000万方，水平向尺寸近千米，规模巨大，在地震作用下应力波传播复杂，需要重点关注地震作用下的动力响应与堆积体的变形破坏机制。但由于模型较大，对于堆积体采用不连续数值方法(pfc3d)模拟,其它岩土体采用连续数值模拟方法(flac3d6.0)实现。但如何处理连续模型与非连续模型间的应力波传递对计算结果非常重要。利用本发明一种连续-非连续数值模型连接处透射边界处理方法进行实现，步骤如下：
[0100]
(1)基于flac3d6.0平台，该平台可与pfc3d相耦合，故可同时实现连续数值模型与非连续数值模型计算。首先将建好的连续数值模型(包括滑坡体需考虑为非连续的部分)导入，如图2所示。模型共分为4个组，自模型底部向上组“1”为微风化岩体、组“2”为风化岩体、组“3”位堆积体与风化岩体间的滑带土、组“4”为滑坡体。
[0101]
模型共由276827个节点和531815个单元构成。其中每个单元都可由图3a)所示的标准六面体单元或尤其退化而成的三棱柱单元(图3b))、四面体单元(图3c))。
[0102]
遍历所有的节点坐标，找出模型的范围，左边界面xmin＝-1400；右边界面xmax＝960；前边界面ymin＝-500；后边界面ymax＝2450；底边界面zmax＝1600.将左、右、前、后、底边界设置为粘滞边界，粘滞力施加公式如下并在每个时步更新:
[0103]
tn＝-ρc
p
vn[0104]
ts＝-ρc
svs
[0105]
(2)针对存在破坏风险的连续模型滑坡体组“4”，借助软件平台的几何操作功能，将组“4”外边界几何输出为非连续模型边界，该边界由一系列三角形集合而成，如图4所示。
[0106]
根据经验设置滑坡体采用半径1.5～3.0m的圆球，估计平均颗粒半径rave的值为2.25m。
[0107]
对非连续模型边界几何进行法向判断，如果法向量为正，则表明边界下半部分，将其位置下移2倍～10倍平均颗粒半径，实际取20m。上边界为边坡面，则上边界不可更改，故上半边界面平移距离0.0m。
[0108]
(2)遍历滑坡体模型范围几何节点，分别找出其x轴最小值为-862.0m，x轴最大值452.0m，y轴向最小值为-262.0m，y轴向最大值为153.0m，z轴向最小值2132m，z轴最大值3248m。利用三个坐标轴向(x轴、y轴、z轴)约束的长方体区域与非连续模型边界几何共同约束，在离散模型边界控制区内生成半径1.5m～3m高斯分布的圆球，共125378个。
[0109]
非连续域边界放大后的几何、连续的单元(zone)、非连续的球(ball)构成的数值模型如图5所示。
[0110]
在上述长方体区域生成半径为10m的测量球，如图6a)，删除孔隙率大于0.5或者测量球中心不再非连续模型范围内的测量球，剩余测量球用于颗粒体系监控，如图6b)。
[0111]
基于离散元计算平台pfc3d，设置圆球颗粒体系间监测有效模量5e8mpa，刚度比2.0，令颗粒体系平衡，利用颗粒伺服方法令颗粒体系满足均匀化并所有测量圆的平均应力为0.5mpa。颗粒体系状态逐步逼近伺服应力如图7所示。
[0112]
(3)遍历所有的实体单元zone和颗粒ball，根据其几何位置将模型分为三部分，连续模型区、离散模型区、模型搭接区，如图8a)所示。其中模型搭接区为边界控制颗粒，该区域是一块重叠区域，一个实体网格单元可能与一个或者多个颗粒重叠。
[0113]
具体操作如图8b)所示，不失一般性设一个单元由1，2，3，4，5，6，7，8个节点构成，对应这8个节点的位移分别是对应这8个节点的速度则位于该单元内的圆球(编号j)位移可以由下式插值得到。
[0114]
重叠区域位移的传递计算如式所示：
[0115][0116]
(5)对于颗粒体系细观参数标定时，应注意模型搭接区颗粒参数体现出的宏观性质应与连续模型的宏观性质一致或接近一致，以保证模型搭接区应力波的稳定。
[0117]
在已知任何时刻t节点的速度时，通过插值计算重叠区域内颗粒所处位置的速度，将该速度赋值给相应位置的球j，进而传递到不连续区域的颗粒体系中。
[0118]
在t+dt时刻颗粒体系所受的力传递到i后，该力作为集中力传递到连续数值模型中参与计算。如此反复，应力波在该边界上可以连续光滑的传递。
[0119]
赋予参数后，施加典型地震波，利用以上原理得到的典型时刻应力波传播如图9所示，可见连续-非连续数值模型边界上应力可圆滑传递，说明采用本技术提供方法所得模型所得数据能有效改善地震作用下应力波的模拟效果，该模型拟合所得数据能更符合地震动力荷载的传递的实际，所得结果更符合力学机理。
[0120]
在每一个时刻，对比该时刻连续网格的位移、速度、加速度与该节点历史最大位移、速度、加速度之对比，如果当前位移、速度、加速度分别大于其历史最大之，则用当前值替换历史最大值。计算终了时各节点的最大位移、速度、加速度即为位移、速度、加速度反应谱值。
[0121]
尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。