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

摘要

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

主权项

颗粒流与有限差分法耦合计算方法，其特征在于具体步骤如下：(1)，模型信息输入及边界节点信息的传输首先确定问题域，并将问题域划分为离散介质区域和连续介质区域，对于离散介质区域采用颗粒流法建立计算模型，对于连续介质区域采用有限差分法建立计算模型，输入计算问题域的标识、模型参数以及外力信息，并确定离散介质与连续介质的接触边界，将处于接触边界上的有限差分网格节点的坐标、速度初始化并传输到颗粒流计算模块，以此进行接触边界网格划分；(2)，接触边界网格函数形式采用自然坐标映射法，运用形函数和节点坐标作为参数表示边界网格，边界网格的空间形态为多参数曲面：$x_i^{\left[C\right]}=N_j\times x_i^j$$N_1=\frac{1}{4}(1+\xi )(1+\eta )$$N_2=\frac{1}{4}(1-\xi )(1+\eta )---\left(1\right)$$N_3=\frac{1}{4}(1-\xi )(1-\eta )$$N_4=\frac{1}{4}(1+\xi )(1-\eta )$其中，为边界网格上任意一点坐标，为接触界面网格节点坐标，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]}=\Sigma 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}}{\omega }_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{\Delta U}}_i^s=V_i^s\mathrm{\Delta t}---\left(9\right)$切向接触力增量：${\mathrm{\Delta F}}_i^s=-k_s{\mathrm{\Delta U}}_i^s---\left(10\right)$切向接触力累积量方向改变：$\left\{\begin{array}{c}{\left\{F_i^s\right\}}_{\mathrm{rot}.1}=F_j^s({\delta }_{\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}({\delta }_{\mathrm{ij}}-e_{\mathrm{ijk}}\left({\omega }_k^{\left[D\right]}\right)\mathrm{\Delta t})\end{array}\right.---\left(11\right)$切向接触力：$F_i^s={\left\{F_i^s\right\}}_{\mathrm{rot}.2}+{\mathrm{\Delta F}}_i^s---\left(12\right)$其中，为接触边界网格各节点速度，为颗粒球心速度，为颗粒转动角速度，Δt为计算时步长，k_s为颗粒切向接触刚度，n^[old]为上计算步颗粒与边界接触面的法向向量；(7)，颗粒与边界网格的滑移模型在颗粒与边界的接触设置滑动模型，防止颗粒与边界之间的切向接触力无限增大，其切向接触力最大值为法向接触力与摩擦系数的乘积；$F_{\max }^s=\mu \left|F^n\right|---\left(13\right)$其中，μ为摩擦系数，Fⁿ为法向接触力；(8)，颗粒流与有限差分法信息的传递颗粒与接触边界网格所受到的接触力为颗粒对接触边界的反作用力，可通过形函数等效作用在接触边界节点所在的有限差分网格节点上，如此循环往复，直至模型达到稳定状态；$F_i^j=N_j\times F_i---\left(14\right).$