CN104091009A - 颗粒流与有限差分法耦合计算方法 - Google Patents

Abstract

本发明涉及一种基于颗粒物质力学和有限差分法的耦合计算方法，属于计算岩土力学技术领域。针对现有离散介质与连续介质计算模拟方法的局限性，采用颗粒流/有限差分网格表面任意位置耦合方法实现两种不同介质边界相互作用力的传递，并以牛顿第二定律为基础，分别对颗粒流和有限差分网格进行计算，建立了考虑离散介质-连续介质完全耦合的岩土材料计算模型。本发明能够有效地实现离散介质和连续介质边界耦合计算，可以多尺度分析岩土体的受力特性，同时有力推动了颗粒流和有限差分法在岩土工程技术领域的应用。

Description

颗粒流与有限差分法耦合计算方法

技术领域

本发明涉及一种离散介质与连续介质计算模拟方法，具体涉及一种颗粒流与有限差分法耦合计算方法，属于计算岩土力学技术领域。

技术背景

目前，在岩土力学数值计算领域，常用的方法有有限元法、有限差分法、颗粒流等方法，有限元法和有限差分法将岩土体视为连续介质对其进行计算，只能对岩土体进行宏观分析，不能对土体颗粒的细观特性进行分析；颗粒流法将岩土体视为离散体，针对岩土体的细观特性进行分析，但考虑到空间尺度问题，颗粒半径与工程模型中最小结构的尺寸应该保持一定的比例(内尺度一般取0.01～0.1)，这就造成颗粒流模型需要上百万的三维颗粒，致使进行颗粒流法计算对计算机硬件要求很高，计算成本很大。因此，迫切需要一种切实的计算方法使有限差分法和颗粒流法能够耦合计算，实现岩土体数值模拟计算过程中局部采用颗粒流分析其细观特性，其余部分采用有限差分法进行分析，分析其整体的宏观特性。

发明内容

本发明的目的是针对当前岩土体数值模拟方法主要采用有限元、有限差分法、颗粒流等单一方法进行模拟现状，以及有限元、有限差分法不能分析岩土体细观特性和颗粒流法计算成本大的局限性，提出了一种颗粒流和有限差分法耦合算法，实现两种不同介质边界相互作用力的传递，并以牛顿第二定律为基础，分别对颗粒流和有限差分网格进行计算，建立了考虑离散介质-连续介质完全耦合的岩土材料计算模型。

技术方案：本发明提出的颗粒流与有限差分法耦合计算方法，包括以下步骤：

(1)，模型信息输入及边界节点信息的传输

首先确定问题域，并将问题域划分为离散介质区域和连续介质区域，对于离散介质区域采用颗粒流法建立计算模型，对于连续介质区域采用有限差分法建立计算模型，输入计算问题域的标识、模型参数以及外力信息，并确定离散介质与连续介质的接触边界，将处于接触边界上的有限差分网格节点的坐标、速度初始化并传输到颗粒流计算模块，以此进行接触边界网格划分；

(2)，接触边界网格函数形式

采用自然坐标映射法，运用形函数和节点坐标作为参数表示边界网格，边界网格的空间

形态为多参数曲面：

$x_i^{\left[C\right]}=N_j\×x_i^j$

$N_1=\frac{1}{4}(1+\mathrm{\ξ})(1+\mathrm{\η})$

$N_2=\frac{1}{4}(1-\mathrm{\ξ})(1+\mathrm{\η})---\left(1\right)$

$N_3=\frac{1}{4}(1-\mathrm{\ξ})(1-\mathrm{\η})$

$N_4=\frac{1}{4}(1+\mathrm{\ξ})(1-\mathrm{\η})$

其中，为边界网格上任意一点坐标，为接触界面网格节点坐标，N_j为接触界面形函数，j＝1～4，ξ、η为映射自然坐标。

(3)，接触边界邻近颗粒全域搜索

采用改进空间网格法进行全域搜索，搜索在可能与接触边界接触的所有邻近颗粒；首先对颗粒流模型空间区域进行空间网格划分，经验最佳空间网格尺寸为每个空间网格包罗颗粒数目为4～20个，然后根据空间网格的空间顺序，依次确定各个空间网格所包罗的颗粒，并确定包罗各个有限差分单元接触面的空间网格集合，从而得到每个接触边界网格所有可能与其接触的邻近颗粒；

(4)，接触边界邻近颗粒的区域搜索

首先对全域搜索过程得到的潜在接触颗粒与边界网格进行空间几何运算，剔除不符合空间接触条件的颗粒，然后采用牛顿迭代法计算颗粒中心到边界网格的投影点自然坐标，从而得到颗粒中心到边界网格的最小距离，并将其与颗粒半径进行比较，若其最小距离小于颗粒半径，则判定颗粒与边界网格接触；

$d=|x^{\left[D\right]}-x^{\left[C\right]}|=\sqrt{(x^{\left[D\right]}-x^{\left[C\right]})(x^{\left[D\right]}-x^{\left[C\right]})}---\left(3\right)$

其中，为原点到颗粒球心的向量，为原点到接触界面上颗粒球心投影点C的向量，d为颗粒球心D到边界网格的投影点C的距离，x^[D]为颗粒球心的坐标，x^[C]为球心D到边界网格的投影点C的坐标；

(5)，颗粒与边界网格的法向接触力

颗粒采用线弹性接触模型，计算颗粒与边界网格法向接触力，首先通过颗粒半径减去球心到边界的最小距离求得颗粒与边界法线方向上的最大重叠量，并求出法向单位向量n，根据颗粒的法向刚度求得颗粒所受到的法向接触力的各向分力；

法向单位向量：

法向接触力各方向分力：F_i ⁿ＝k_nUⁿn_i                        (5)

其中，n为颗粒与边界接触平面的法向单位向量，x^[D]、y^[D]、z^[D]为颗粒球心在坐标各轴的分量，x^[C]、y^[C]、z^[C]为颗粒与边界接触点C在坐标各轴的分量，为沿坐标各轴方向的单位向量，为法向接触力各方向分力，k_n为颗粒法向接触刚度，Uⁿ为颗粒与接触边界重叠量，n_i为法向单位向量各方向分量；

(6)，颗粒与边界网格的切向接触力

切向接触力采用增量叠加的方式计算，在单个时间步内，增量值采用颗粒与边界接触点垂直于法线平面的相对位移与颗粒切向刚度相乘计算，切向接触力方向的改变主要通过新旧两个接触平面的公共线方向、新接触平面的法向方向，两个方面来实现；

边界接触点速度： ${\stackrel{.}{x}}_i^{\left[C\right]}=\mathrm{\Σ}{N}_j{\stackrel{.}{x}}_i^j,(j=\mathrm{1,2,3,4})---\left(6\right)$

相对速度： $V_i={\stackrel{.}{x}}_i^{\left[D\right]}-{\stackrel{.}{x}}_k^{\left[C\right]}+e_{\mathrm{ijk}}{\mathrm{\ω}}_j^{\left[D\right]}(x_k^{\left[C\right]}-x_k^{\left[D\right]})---\left(7\right)$

切向相对速度：V_i ^s＝V_i-V_i ⁿ＝V_i-V_jn_jn_i   (8)

切向相对位移： ${\mathrm{\ΔU}}_i^s=V_i^s\mathrm{\Δt}---\left(9\right)$

切向接触力增量： ${\mathrm{\ΔF}}_i^s=-k_s{\mathrm{\ΔU}}_i^s---\left(10\right)$

切向接触力累积量方向改变： $\left\{\begin{array}{c}{\left\{F_i^s\right\}}_{\mathrm{rot}.1}=F_j^s({\mathrm{\δ}}_{\mathrm{ij}}-e_{\mathrm{ijk}}{e}_{\mathrm{kmn}}{n}_m^{\left[\mathrm{old}\right]}{n}_m)\\ {\left\{F_i^s\right\}}_{\mathrm{rpt}.2}={\left\{F_i^s\right\}}_{\mathrm{rot}.1}({\mathrm{\δ}}_{\mathrm{ij}}-e_{\mathrm{ijk}}\left({\mathrm{\ω}}_k^{\left[D\right]}\right)\mathrm{\Δt})\end{array}\right.---\left(11\right)$

切向接触力： $F_i^s={\left\{F_i^s\right\}}_{\mathrm{rot}.2}+{\mathrm{\ΔF}}_i^s---\left(12\right)$

其中，为接触边界网格各节点速度，为颗粒球心速度，为颗粒转动角速度，Δt为计算时步长，k_s为颗粒切向接触刚度，n^[old]为上计算步颗粒与边界接触面的法向向量；

(7)，颗粒与边界网格的滑移模型

在颗粒与边界的接触设置滑动模型，防止颗粒与边界之间的切向接触力无限增大，其切向接触力最大值为法向接触力与摩擦系数的乘积；

$F_{\max }^s=\mathrm{\μ}\left|F^n\right|---\left(13\right)$

其中，μ为摩擦系数，Fⁿ为法向接触力；

(8)，颗粒流与有限差分法信息的传递

颗粒与接触边界网格所受到的接触力为颗粒对接触边界的反作用力，可通过形函数等效作用在接触边界节点所在的有限差分网格节点上，如此循环往复，直至模型达到稳定状态；

$F_i^j=N_j\×F_i---\left(14\right).$

与现有技术相比，本发明的优势如下：

(1)当前在计算岩土力学技术领域，尚无系统的岩土材料采用颗粒流-有限差分法耦合计算模拟的方法，通过本发明的计算方法，能够有效地将两种数值方法进行耦合计算，在提高计算效率的同时，可以精确地分析局部颗粒流模拟区域的细观特性及岩土体整体的宏观特性，为有限差分法和颗粒流法在岩土力学中的应用提供了有力的技术支撑。

(2)本发明首次采用改进空间网格法进行颗粒流-有限差分法接触边界的全域搜索。改进空间网格搜索方法实现了颗粒在不同网格中的移动，消除了颗粒粒径对网格尺寸的影响，网格划分尺寸不受颗粒大小限制，有效地解决了大尺寸粒径颗粒的检索问题，很大程度上提高了计算效率。

(3)对于颗粒中心与接触边界的最小距离计算，属于空间点到多参数曲面的距离问题，本发明采用牛顿迭代法求解其自然映射坐标表示的微分方程组获得其接触点，进而求解其最小距离，该算法计算效率高，计算结果精确。

(4)目前现有的耦合计算技术，一般对耦合计算区域划定范围，往往局限于小变形问题分析，本发明采用形函数和节点坐标描述接触边界曲面，在计算过程中随着节点坐标的变化，边界曲面的空间位置也在不断变化，并且不存在范围问题，因此，本发明不仅适用于小变形问题，在极大变形分析上同样具有较高的计算精度和计算效率。

(5)目前现有的技术对于离散介质-连续介质耦合的方法，一般采用过渡区进行耦合或仅考虑法向接触力，往往会导致耦合计算精度的下降。本发明打破了这一局限性，对于颗粒与接触边界的接触力计算，本发明将颗粒与接触边界的接触力分为法向接触力和切向接触力，法向接触力通过颗粒与接触边界的最大重叠量计算，切向接触力采用累加计算，并且设置了滑移接触模型，使耦合计算模型更符合岩土力学原理，计算结果更加精确。

附图说明：

图1为本发明的流程图。

图2为与边界网格接触颗粒的搜索流程图。

图3为边坡案例计算模型初始状态图。

图4为模型计算30000步时边坡滑坡运动的计算模拟结果。

图5为模型计算60000步时边坡滑坡运动的计算模拟结果。

图6为模型计算100000步时边坡滑坡运动的计算模拟结果。

图7为模型计算30000步时边坡土体位移矢量和基覆面剪切速率图。

图8为模型计算60000步时边坡土体位移矢量和基覆面剪切速率图。

图9为模型计算100000步时边坡土体位移矢量和基覆面剪切速率图。

具体实施方式

本实施例仅为介绍本发明的使用方式和有效性，其具体的数据及内容不在本发明保护范围之内。

实施例1：岩体基覆面上边坡滑坡运动过程模拟

计算模型的基岩倾角为10°，岩体为花岗岩；边坡堆积体材料为砂土，边坡坡脚为30°，边坡高度为1.15m，对边坡土体强度进行折减，折减系数为1.67，边坡失稳滑坡。

应用本发明的计算方法，对边坡滑坡过程进行计算。具体计算过程如图1和图2所示：

(1)首先确定边坡土体与基岩范围，并将问题域划分为离散介质区域和连续介质区域，对于边坡土体部分采用颗粒流法建立计算模型，对于基岩部分采用有限差分法建立计算模型(如图3所示)，输入计算问题域的标识、模型物理、力学参数(如表1所示)，并确定边坡土体与基岩的接触边界，将处于接触边界上的有限差分网格节点的坐标、速度初始化并传输到颗粒流计算模块，以此进行接触边界网格划分；

表1边坡模型材料计算参数

(2)采用自然坐标映射法，运用形函数和节点坐标作为参数表示边界网格，边界网格的空间形态为多参数曲面(公式1)；

(3)采用改进空间网格法进行全域搜索，根据空间网格的空间顺序，依次确定各个空间网格所包罗的颗粒，并确定包罗各个有限差分单元接触面的空间网格集合，从而得到每个接触边界网格所有可能与其接触的邻近颗粒；

(4)首先对全域搜索过程得到的潜在接触颗粒与边界网格进行空间几何运算，剔除不符合空间接触条件的颗粒，然后采用公式2、3计算颗粒中心到边界网格的最小距离，并将其与颗粒半径进行比较，若其最小距离小于颗粒半径，则判定颗粒与边界网格接触；

(5)颗粒采用线弹性接触模型，通过公式4、5计算颗粒与边界网格法向接触力；

(6)通过公式(6～12)采用增量叠加的方式计算颗粒与边界网格的切向接触力；

(7)在颗粒与边界的接触设置滑动模型，防止颗粒与边界之间的切向接触力无限增大，其切向接触力最大值为法向接触力与摩擦系数的乘积，对步骤(6)所得的切向接触力进行修正。

(8)计算颗粒与接触边界网格所受到的接触力，颗粒对接触边界的反作用力，可通过形函数等效作用(公式14)在接触边界节点所在的有限差分网格节点上，如此循环往复，直至模型达到稳定状态(如图4～6)

(9)对模拟计算数据进行后处理，分析边坡滑移位移趋势以及对基岩的影响(图7～9)，结果证明边坡滑移位移趋势与基岩剪切速率趋势吻合，并且结果与有限元软件刚度折减方法得到的边坡滑移位移趋势与基岩剪切速率相吻合，证明本发明对于颗粒流与有限差分法耦合计算的精度满足计算要求。

Claims (1)

1.颗粒流与有限差分法耦合计算方法，其特征在于具体步骤如下：

(1)，模型信息输入及边界节点信息的传输

首先确定问题域，并将问题域划分为离散介质区域和连续介质区域，对于离散介质区域采用颗粒流法建立计算模型，对于连续介质区域采用有限差分法建立计算模型，输入计算问题域的标识、模型参数以及外力信息，并确定离散介质与连续介质的接触边界，将处于接触边界上的有限差分网格节点的坐标、速度初始化并传输到颗粒流计算模块，以此进行接触边界网格划分；

(2)，接触边界网格函数形式

采用自然坐标映射法，运用形函数和节点坐标作为参数表示边界网格，边界网格的空间形态为多参数曲面：

$x_i^{\left[C\right]}=N_j\×x_i^j$

$N_1=\frac{1}{4}(1+\mathrm{\ξ})(1+\mathrm{\η})$

$N_2=\frac{1}{4}(1-\mathrm{\ξ})(1+\mathrm{\η})---\left(1\right)$

$N_3=\frac{1}{4}(1-\mathrm{\ξ})(1-\mathrm{\η})$

$N_4=\frac{1}{4}(1+\mathrm{\ξ})(1-\mathrm{\η})$

其中，为边界网格上任意一点坐标，为接触界面网格节点坐标，N_j为接触界面形函数，j＝1～4，ξ、η为映射自然坐标；

(3)，接触边界邻近颗粒全域搜索

采用改进空间网格法进行全域搜索，搜索在可能与接触边界接触的所有邻近颗粒；首先对颗粒流模型空间区域进行空间网格划分，经验最佳空间网格尺寸为每个空间网格包罗颗粒数目为4～20个，然后根据空间网网格的空间顺序，依次确定各个空间网格所包罗的颗粒，并确定包罗各个有限差分单元接触面的空间网格集合，从而得到每个接触边界网格所有可能与其接触的邻近颗粒；

(4)，接触边界邻近颗粒的区域搜索

首先对全域搜索过程得到的潜在接触颗粒与边界网格进行空间几何运算，剔除不符合空间接触条件的颗粒，然后采用牛顿迭代法计算颗粒中心到边界网格的投影点自然坐标，从而得到颗粒中心到边界网格的最小距离，并将其与颗粒半径进行比较，若其最小距离小于颗粒半径，则判定颗粒与边界网格接触；

$d=|x^{\left[D\right]}-x^{\left[C\right]}|=\sqrt{(x^{\left[D\right]}-x^{\left[C\right]})(x^{\left[D\right]}-x^{\left[C\right]})}---\left(3\right)$

其中，为原点到颗粒球心的向量，为原点到接触界面上颗粒球心投影点C的向量，d为颗粒球心D到边界网格的投影点C的距离，x^[D]为颗粒球心的坐标，x^[C]为球心D到边界网格的投影点C的坐标；

(5)，颗粒与边界网格的法向接触力

颗粒采用线弹性接触模型，计算颗粒与边界网格法向接触力，首先通过颗粒半径减去球心到边界的最小距离求得颗粒与边界法线方向上的最大重叠量，并求出法向单位向量n，根据颗粒的法向刚度求得颗粒所受到的法向接触力的各向分力；

法向单位向量

法向接触力各方向分力：F_i ⁿ＝k_nUⁿn_i                          (5)

其中，n为颗粒与边界接触平面的法向单位向量，x^[D]、y^[D]、z^[D]为颗粒球心在坐标各轴的分量，x^[C]、y^[C]、z^[C]为颗粒与边界接触点C在坐标各轴的分量，为沿坐标各轴方向的单位向量，为法向接触力各方向分力，k_n为颗粒法向接触刚度，Uⁿ为颗粒与接触边界重叠量，n_i为法向单位向量各方向分量；

(6)，颗粒与边界网格的切向接触力

切向接触力采用增量叠加的方式计算，在单个时间步内，增量值采用颗粒与边界接触点垂直于法线平面的相对位移与颗粒切向刚度相乘计算，切向接触力方向的改变主要通过新旧两个接触平面的公共线方向、新接触平面的法向方向，两个方面来实现；

边界接触点速度： ${\stackrel{.}{x}}_i^{\left[C\right]}=\mathrm{\Σ}{N}_j{\stackrel{.}{x}}_i^j,(j=\mathrm{1,2,3,4})---\left(6\right)$

相对速度： $V_i={\stackrel{.}{x}}_i^{\left[D\right]}-{\stackrel{.}{x}}_k^{\left[C\right]}+e_{\mathrm{ijk}}{\mathrm{\ω}}_j^{\left[D\right]}(x_k^{\left[C\right]}-x_k^{\left[D\right]})---\left(7\right)$

切向相对速度：V_i ^s＝V_i-V_i ⁿ＝V_i-V_jn_jn_i   (8)

切向相对位移： ${\mathrm{\ΔU}}_i^s=V_i^s\mathrm{\Δt}---\left(9\right)$

切向接触力增量： ${\mathrm{\ΔF}}_i^s=-k_s{\mathrm{\ΔU}}_i^s---\left(10\right)$

切向接触力累积量方向改变： $\left\{\begin{array}{c}{\left\{F_i^s\right\}}_{\mathrm{rot}.1}=F_j^s({\mathrm{\δ}}_{\mathrm{ij}}-e_{\mathrm{ijk}}{e}_{\mathrm{kmn}}{n}_m^{\left[\mathrm{old}\right]}{n}_m)\\ {\left\{F_i^s\right\}}_{\mathrm{rpt}.2}={\left\{F_i^s\right\}}_{\mathrm{rot}.1}({\mathrm{\δ}}_{\mathrm{ij}}-e_{\mathrm{ijk}}\left({\mathrm{\ω}}_k^{\left[D\right]}\right)\mathrm{\Δt})\end{array}\right.---\left(11\right)$

切向接触力： $F_i^s={\left\{F_i^s\right\}}_{\mathrm{rot}.2}+{\mathrm{\ΔF}}_i^s---\left(12\right)$

其中，为接触边界网格各节点速度，为颗粒球心速度，为颗粒转动角速度，Δt为计算时步长，k_s为颗粒切向接触刚度，n^[old]为上计算步颗粒与边界接触面的法向向量；

(7)，颗粒与边界网格的滑移模型

在颗粒与边界的接触设置滑动模型，防止颗粒与边界之间的切向接触力无限增大，其切向接触力最大值为法向接触力与摩擦系数的乘积；

$F_{\max }^s=\mathrm{\μ}\left|F^n\right|---\left(13\right)$

其中，μ为摩擦系数，Fⁿ为法向接触力；

(8)，颗粒流与有限差分法信息的传递

颗粒与接触边界网格所受到的接触力为颗粒对接触边界的反作用力，可通过形函数等效作用在接触边界节点所在的有限差分网格节点上，如此循环往复，直至模型达到稳定状态；

$F_i^j=N_j\×F_i---\left(14\right).$