[发明专利]颗粒流与有限差分法耦合计算方法

权利要求书

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

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

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

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

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

xi[C]=Nj×xij]]>

N1=14(1+ξ)(1+η)]]>

N2=14(1-ξ)(1+η)---(1)]]>

N3=14(1-ξ)(1-η)]]>

N4=14(1+ξ)(1-η)]]>

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

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

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

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

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

d=|x[D]-x[C]|=(x[D]-x[C])(x[D]-x[C])---(3)]]>

其中，为原点到颗粒球心的向量，为原点到接触界面上颗粒球心投影点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)，颗粒与边界网格的切向接触力

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

边界接触点速度：x.i[C]=ΣNjx.ij,(j=1,2,3,4)---(6)]]>

相对速度：Vi=x.i[D]-x.k[C]+eijkωj[D](xk[C]-xk[D])---(7)]]>

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

切向相对位移：ΔUis=VisΔt---(9)]]>

切向接触力增量：ΔFis=-ksΔUis---(10)]]>

切向接触力累积量方向改变：{Fis}rot.1=Fjs(δij-eijkekmnnm[old]nm){Fis}rpt.2={Fis}rot.1(δij-eijk(ωk[D])Δt)---(11)]]>

切向接触力：Fis={Fis}rot.2+ΔFis---(12)]]>

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

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

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

Fmaxs=μ|Fn|---(13)]]>

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

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

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

Fij=Nj×Fi---(14).]]>