[发明专利]一种基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法

权利要求书

1.一种基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法，其特征在于，该方法包括以下步骤：

步骤1、将流场划分为欧拉网格，并设定欧拉网格中的初始参数，得到LBM流场；根据所要模拟的堆积体块体的级配，编程生成代表堆积体块体几何形状的块体，并随机选择产生位置，块体边界区域内的网格设为拉格朗日网格，定义块体的接触参数，得到DDA块体；进而得到LBM-DDA流固耦合模型；

步骤2、在LBM-DDA流固耦合模型中加入插值格式的边界条件方程，将DDA块体边界对欧拉网格的作用离散到整个LBM流场中，得到DDA块体对LBM流场的单向作用；

步骤3、完成DDA块体对LBM流场的单向作用计算后，进行LBM流场的欧拉网格之间的碰撞和迁移，完成LBM流场的计算；

步骤4、根据LBM流场，应用基于伽利略不变的动量交换法计算LBM流场作用于DDA块体的流体曳力，再将流体曳力作为载荷施加到相应的DDA块体上；

步骤5、DDA块体的接触搜索：判定DDA块体间的相对状态；相对状态包括接触状态和张开状态，接触状态包括锁定状态和滑动状态；

若判定结果为接触状态，则设置接触弹簧形成的接触子矩阵，并对DDA块体间的接触面上的粘聚力进行计算；若判定结果为张开状态，则不设置接触子矩阵，DDA块体间不存在接触面上的粘聚力；在粘聚力的计算中，如果DDA块体间出现过滑动状态或张开状态，则此DDA块体间的接触面上的粘聚力消失，在后续的时间步中此DDA块体间永远不存在粘聚力；

步骤6、DDA块体的动力学计算：根据步骤5的判定结果，当DDA块体间的相对状态为接触状态时，计算DDA块体的单元刚度矩阵和外力矩阵，得到DDA块体产生相对滑动或碰撞时的接触力；然后将该接触力与步骤4得到的流体曳力进行合成，得到合力以及合力矩；再根据牛顿第二定律求出DDA块体的加速度，再对时间积分进而求得DDA块体的速度和位移，从而不断累加，实现DDA块体的大位移计算；

步骤7、首先，判断是否到达最大计算时间，若没有达到则返回步骤2，若达到，首先通过输出DDA块体的坐标，获取渗流侵蚀破坏后的LBM-DDA流固耦合模型的孔隙率以及模型受到渗流侵蚀破坏后损失的DDA块体的质量和直径分布；然后，通过LBM流场得到该模型在受到渗流侵蚀破坏后的渗透系数，进而得到渗流侵蚀破坏后的LBM-DDA流固耦合模型的计算结果。

2.根据权利要求1所述的基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法，其特征在于，步骤1中，欧拉网格的初始参数包括分布密度函数和速度边界条件。

3.根据权利要求2所述的基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法，其特征在于，步骤1中，欧拉网格中的分布密度函数由平衡分布密度函数确定，平衡分布密度函数如式(1)所示：

式(1)中，f_i^eq为平衡分布密度函数；ρ为欧拉网格的密度，f_i为分布密度函数；ω_i为权重；e_i为离散速度矢量；v＝(v_x,v_y)代表出入口流速。

4.根据权利要求1所述的基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法，其特征在于，步骤2中，插值格式的边界条件方程如式(2)所示：

式(2)中，f_i'和f′_inv(i)为完成碰撞之后的分布密度函数，下标inv(i)和i代表其速度方向相反；q为某个欧拉网格到DDA块体边界的直线距离；δf_i'为DDA块体表面运动引起的速度矫正：

式(3)中，c_s为网格速度，u_ω为DDA块体边界处的速度。

5.根据权利要求1所述的基于LBM-DDA耦合的堆积体渗流侵蚀破坏模拟计算方法，其特征在于，步骤3中，碰撞和迁移过程是：对每个欧拉网格中的分布密度函数在每个时间步长内通过碰撞进行重新分配，再将重新分配后的分布密度函数沿离散速度方向迁移到最近的欧拉网格上。