[发明专利]一种基于DEM-CFD耦合计算的滑坡-堵江-涌浪灾害链模拟方法

权利要求书

1.一种基于DEM-CFD耦合计算的滑坡-堵江-涌浪灾害链模拟方法，其特征在于，包括以下步骤：

第一步，利用数字高程模型对滑坡-堵江-涌浪灾害链区域进行三维建模，并通过笛卡尔坐标系对三维模型计算域进行网格剖分，设置边界命名并导出网格文件；

第二步，确定滑坡体模拟参数，其中，滑坡体模拟参数包括材料参数及接触参数；

第三步，将第一步导出的网格文件导入DEM模块并建立滑坡体源区，并根据所确定滑坡体模拟参数，在滑坡体源区中填充与实际方量相当的颗粒以模拟滑坡体；

第四步，将第一步导出的网格文件导入CFD模块，并进行流体材料参数设定及流体域流体填充；所述流体材料包括水及空气，所述流体材料参数包括流体密度及流体粘度；所述流体域流体填充根据河流实际水面高程确定，其中水面高程以下流体材料由水填充，水面高程以上流体材料由空气填充；

第五步，根据河流实际入流流量确定CFD模块中的流体域边界条件，包括上流入流边界条件，下游压力出流边界条件及顶部开放边界条件；

第六步，根据高程方向确定重力方向，并设置重力加速度；

第七步，设置模拟时长，分别设置所述DEM模块及CFD模块时间步长，并确定耦合频率；所述耦合频率为CFD模块的时间步长与DEM模块时间步长的比值，其数值应在1-100以内；所述模拟时长应根据实际案例需求确定；

第八步，根据第二至第七部设置初始化DEM模块及CFD模块，并通过双向耦合计算开始模拟滑坡-堵江-涌浪灾害链演化过程；在双向耦合计算中，在DEM模块计算颗粒作用力与位移；在CFD模块中计算流场演化，并在CFD模块中引入VOF模型，该模型允许引入空气相及水相的多相流体，进一步通过追踪空气-水相交界面来模拟河流演化；耦合作用模块分别读取DEM模块的颗粒位置和速度数据，及CFD模块的流体压力场及速度场数据，进一步通过虚拟球模型计算局部孔隙率，在虚拟球模型中，通过基于真实DEM颗粒构建虚拟球，从而在CFD网格中确定孔隙率计算的范围，进一步计算耦合作用力并分别传递回DEM模块及CFD模块。

2.根据权利要求1所述的一种基于DEM-CFD耦合计算的滑坡-堵江-涌浪灾害链模拟方法，其特征在于，第八步中，所述双向耦合计算具体包括：

8.1)在DEM模块中通过Hertz-Mindlin接触模型计算颗粒间及颗粒与边界间接触力；

8.2)DEM模块通过牛顿第二定律求解式(1)及式(2)的运动方程：

式中，m_i为颗粒质量，v_i为颗粒速度，ω_i为颗粒角速度，F_ij和M_ij分别为颗粒的接触力及扭矩，F_f为流体对颗粒作用力，g为重力加速度；

8.3)在DEM模块中更新求解运动方程得到的颗粒位置及速度信息，重复上述步骤至单次耦合计算内DEM循环次数达到耦合频率次数；

8.4)颗粒位置及速度信息导入耦合作用模块，通过虚拟球模型计算局部孔隙率，并更新动量源项；进一步的，所述虚拟球模型计算局部孔隙率包括下述步骤：

根据真实DEM颗粒的位置构建相同圆心，且直径为4倍颗粒直径的虚拟球，在虚拟球覆盖范围内的CFD网格确定真实DEM颗粒的孔隙率计算关联网格；随后将颗粒体积均匀划分至虚拟球所确定孔隙率计算关联网格；同一CFD网格允许关联至不同DEM颗粒，即存在重叠孔隙率计算关联网格，因此某一网格内局部孔隙率可由式(3)计算：

式(1)中，ε为局部孔隙率，V_p为颗粒i体积，V_rc为颗粒i的孔隙率计算关联网格总体积；

8.5)耦合作用模块将局部孔隙率及动量源项导入至CFD模块；

8.6)在CFD模块中首先通过考虑颗粒相体积分数的VOF模型求解基于流体速度的输运方程，追踪河流自由液面演化过程实现滑坡涌浪传播的模拟：

式(4)中，α₁为水相体积分数，为u流体速度，u_r为相对速度，也称压缩速度；

8.7)在CFD模块中通过PISO算法，求解流体控制方程，所述流体控制方程包括动量方程(5)及连续性方程(6)，根据所述流体控制方程的求解结果更新流体域速度场及压力场，直至收敛；

8.8)将流体域速度场、压力场及流体材料参数传输至耦合作用模块，随后在耦合作用模块中计算耦合作用力并传递至DEM模块；重新返回至8.1步开始下一个计算循环，直至达到模拟时长。