[发明专利]基于GPU矩阵的离散元流固耦合数值模拟方法及系统

权利要求书

1.一种基于GPU矩阵的离散元流固耦合数值模拟方法，其特征在于，包括以下步骤：

步骤1，建立孔隙网络模型，包括建立网格编号矩阵：建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒；在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C；利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN；然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙；

步骤2，建立流体单元运移模型：

步骤2.1，在初始模型中设定初始流体单元状态以及边界条件；

步骤2.2，通过水压与温度以及密度关系式计算各流体单元的密度ρ_i，i表示流体单元的编号；

步骤2.3，计算各流体单元的质量M_i＝ρ_iV_i，V_i为孔隙的面积；

步骤2.4，计算模型中的孔喉渗透系数：孔喉渗透系数与通道长度以及宽度有关，孔喉渗透系数：

其中，k(d)表示孔喉渗透系数，k为设定的数值，L表示两颗粒圆心之间的距离，R₁、R₂表示两颗粒的半径，minR表示所有颗粒的最小半径，l表示孔喉通道长度；

步骤2.5，计算单个流体单元的渗流量：在离散元计算的一个时间步dt内，某个流体单元pore_i与周围流体单元pore_n发生流体运移，流体运移服从达西定律，定义流体流入为正，单个孔喉由于压力梯度导致体积渗流量：

q_i＝K(d)*A*J*dT

J＝(P_n-P_i)/min(R_a，R_b)

其中，P_n和P_i分别是流体单元pore_n和pore_i的压力，R_a和R_b分别是组成该孔喉的两颗粒a和b的半径，q_i表示渗流量，A表示孔喉通道宽度，J表示压力梯度系数，dT表示时间步；

步骤2.6，计算总运移量：pore_i在一个时间步内与周围流体单元进行的运移总量为：

其中，各流体单元的密度V_i为孔隙的体积，M_i为孔隙内各流体单元的质量；

步骤2.7，计算运移后流体单元质量：pore_i在第一个时间步后的总质量：

M′_i＝M_i+Q_n

其中，M′_i表示流体单元发生流体运移后的质量，M_i表示运移前的质量，Q_n表示总运移量；

步骤2.8，更新流体单元状态：一个时间步dt后流体单元的密度流体单元的压力通过函数关系式：P′_i＝f(ρ′_i，T_i)进行更新，其中，ρ′_i表示流体单元的密度，M′_i表示流体单元的质量；

步骤3，流体固体耦合作用：孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用；

步骤3.1，固体颗粒重新平衡：由于流体的运移，各流体单元内的水压发生改变，从而各流体单元对固体颗粒之间的力也发生改变，颗粒受流体单元作用，在一个时间步内的位移距离为：

其中，x_i为在一个时间步内的位移距离，P_i为流体单元的压力，l_i为流体单元与中心颗粒的接触长度，F_i为颗粒间的接触作用力，m₀为固体颗粒的质量；

步骤3.2，孔隙体积更新：由于固体颗粒的位置改变，所以流体单元的体积发生了改变，V_i′为固体颗粒改变之后孔隙的体积；

步骤3.3，流体单元状态更新：计算孔隙体积变化后的流体单元密度，在一个时间步后其密度V′₁为固体颗粒改变之后孔隙的体积；压力为P′_i＝f(ρ″_i，T1)，其中，ρ″_i表示流体单元的密度，P″_i表示流体单元的压力；

步骤3.4，重复步骤1到步骤3.3，即可实现二维流固耦合数值模拟。