一种基于CFD-DEM耦合的颗粒运动模拟方法与流程

本发明属于流体力学的，具体涉及一种基于cfd-dem耦合的颗粒运动模拟方法。

背景技术：

1、颗粒物的输运现象广泛存在于各类流化过程、电池热失控喷发、气力输送、污染物扩散控制等诸多应用领域。如何通过数值模拟获取颗粒物行为，对主要工业设备进行设计优化，提高设备运行效率，增强安全性，降低研发成本具有重要现实意义。

2、目前，将cfd计算流体力学方法与dem离散元方法耦合模拟含颗粒的气固两相流过程是研究流场中颗粒物行为的重要技术手段。例如，采用有限体积法fvm计算流场的cfd-dem耦合方法。然而，该现有技术存在几何精度上的限制：首先，对于流场几何要求较高，要求几何形状尽可能规则，对于实际几何需要做出一定简化以便划分计算网格，容易损失一些流场关键信息。其次，精度难以提高到二阶计算精度以上。

技术实现思路

1、本发明的目的在于提供一种基于cfd-dem耦合的颗粒运动模拟方法，旨在解决上述的问题。本发明基于pjfnk方法求解流场的cfd-dem耦合，可以有效提高cfd-dem耦合计算精度，节约计算内存，并具备优良的扩展性与适应性，适用于电池喷热失控喷发颗粒物在复杂的包内环境中运动特性模拟。

2、本发明主要通过以下技术方案实现：

3、一种基于cfd-dem耦合的颗粒运动模拟方法，包括以下步骤：

4、步骤s1：获取流场和颗粒物的初始信息，所述流场的初始信息包括流体的初始温度、速度和物性条件，所述颗粒物的初始信息包括颗粒物的初始位置、温度、物性参数和形貌参数；对于流体而言，建立考虑颗粒相体积分数和颗粒物对流体反作用力影响的流体动力学方程组；对于颗粒而言，建立描述颗粒受力与运动，并考虑流体对颗粒作用力的颗粒动力学方程组；基于流场与颗粒物的初始信息，对所述流体动力学方程组和所述颗粒动力学方程组进行初始化；

5、步骤s2：对流场划分有限元计算网格，基于有限元方法离散流体动力学方程组，获取以网格节点参数值为求解量的非线性方程组；

6、步骤s3：采用pjfnk方法求解非线性方程组，对当前时间步流场的运动状态进行求解，获取网格节点上流场的速度、温度、密度；

7、步骤s4：计算每一个有限元网格节点与颗粒中心位置的距离，并通过核函数计算颗粒与流体之间相互作用力的权重系数，再计算颗粒物与流体的相互作用力；将颗粒物与流体的相互作用力和权重系数结合，反馈到dem与cfd的动量源项计算中，实现dem与cfd的耦合；然后，计算每一个网格节点上所施加的作用力的总和，并加载在颗粒上，作为颗粒的所受流场作用力的合力；

8、步骤 s5：基于离散元方法求解颗粒的受力、速度、位置、温度，然后基于解得的颗粒的受力、速度、位置、温度对步骤s1中的颗粒动力学方程组进行求解，根据颗粒动力学方程组的求解结果更新颗粒的状态。

9、进一步地，所述流体的物性条件包括流体的密度、粘性、导热系数。所述颗粒物的物性参数包括密度、比热容、弹性模量、泊松比。所述颗粒物的形貌参数为：若是球形颗粒，则形貌参数是半径；若是椭球型颗粒，则形貌参数是长短轴长；若是柱形颗粒，则形貌参数是颗粒的长宽高。

10、为了更好地实现本发明，进一步地，所述步骤s1中，所述流体动力学方程组的计算公式如下：

11、质量连续方程：

12、（1）

13、含颗粒作用力源项的动量守恒方程：

14、（2）

15、含颗粒传热量的能量守恒方程：

16、（3）

17、其中：

18、 ε g为流体的局部空隙率；

19、 ρ g为流体的密度；

20、 u g为流体的表观速度；

21、 τ g为应力张量；

22、 s为流体与颗粒相互作用力；

23、 c pg为流体的比热容；

24、 k g为流体的导热系数；

25、 a为流体与颗粒的接触面积；

26、 t p为颗粒的温度；

27、 t g为流体的局部温度；

28、 h pg为颗粒与流体之间的传热系数；

29、为梯度算子，对变量进行求梯度运算；

30、 t为时间。

31、为了更好地实现本发明，进一步地，所述步骤s1中，所述颗粒动力学方程组的计算公式如下：

32、受力平衡方程：

33、（4）

34、力矩平衡方程：

35、（5）

36、其中， m为颗粒的质量；

37、 v为颗粒速度；

38、 g为重力加速度；

39、 t为时间；

40、 f j为第 j个接触颗粒施加的相互作用力；

41、 i为颗粒的转动惯量；

42、 ω为颗粒的角速度；

43、 t tj为第 j个接触颗粒所作用的切向力矩；

44、 t rj为第 j个接触颗粒所作用的摩擦力矩；

45、 s'为颗粒所受流体作用力。

46、为了更好地实现本发明，进一步地，所述步骤s2中，采用间断有限元方法对流体动力学方程组进行离散，得到离散化的非线性方程组，包括以下步骤：

47、步骤a1：将流场控制方程改写为矢量形式：

48、（6）

49、其中， u为含孔隙率的流体密度矢量：

50、（7）

51、其中， u g、 v g、 w g分别为流体表观速度在空间坐标系中的三个分量；

52、其中， f1( u)、 f2( u)、 f3( u)分别为流体动力学方程组中的对流项：

53、，，（8）

54、其中， p为流体的压力，

55、 g1( u)、 g2( u)、 g3( u)分别为流体动力学方程组中的粘性通量：

56、，，（9）

57、其中， τ ij为粘性应力张量，下标 i和 j分别表示力的作用方向，通过 i和 j两个下标表示了空间中粘性力的作用面和方向；

58、下标 i表示粘性力垂直的方向轴， i=1表示粘性力垂直 x轴， i=2表示粘性力垂直 y轴， i=3表示粘性力垂直 z轴；

59、下标 j表示粘性力沿着的方向， j=1表示粘性力沿着 x轴方向， j=2表示粘性力沿着 y轴方向， j=3表示粘性力沿着 z轴方向；

60、 τ11表示粘性力垂直 x轴，沿 x轴方向；

61、 τ12表示粘性力垂直 x轴，沿 y轴方向；

62、 τ13表示粘性力垂直 x轴，沿 z轴方向；

63、 τ21表示粘性力垂直 y轴，沿 x轴方向；

64、 τ22表示粘性力垂直 y轴，沿 y轴方向；

65、 τ23表示粘性力垂直 y轴，沿 z轴方向；

66、 τ31表示粘性力垂直 z轴，沿 x轴方向；

67、 τ32表示粘性力垂直 z轴，沿 y轴方向；

68、 τ33表示粘性力垂直 z轴，沿 z轴方向；

69、 sr( u)是体积力与热源项，其中包含流体与颗粒的相互作用力：

70、（10）

71、其中， s x、 s y、 s z分别为x,y,z三个方向上流体的体积力源项；

72、步骤a2：将矢量形式采用张量指标记法进行简写：

73、（11）

74、其中， f i( u)为对流项，

75、 g k( u)为扩散项，

76、u为变量向量，

77、 x i和 x k均表示为坐标分量，公式中采用张量写法；

78、通过划分有限元网格，对流场几何进行空间离散，在网格中选取高阶单元及其形函数，在单元上乘以试函数 w test，然后积分可以得到：

79、（12）

80、（13）

81、其中， u h为有限元方法中的形函数；

82、 a i为有限元节点上的待求解物理量的值；

83、 ω为有限单元的体积域；

84、 n为高斯积分点个数；

85、 n为有限元节点个数；

86、σ为有限元单元表面的积分单元；

87、将 h( f( u), g( u))记为 h( u)，且 h( u)为对流项与粘度项散度之和；

88、步骤a3：对公式（13）中的单元边界上的通量项采用br2格式进行离散：

89、

90、其中：

91、为界面左侧的插值函数；

92、为界面右侧的插值函数；

93、步骤a4：对于单元上的积分计算采用高斯积分法求解，获得局部单元矩阵，然后再组装获得总刚度阵，与源项构成关于网格节点上变量的方程组：

94、（14）

95、其中，

96、 v h为有限元为物理场的有限元近似解；

97、 w i高斯积分点权重，

98、 x为网格节点空间坐标。

99、为了更好地实现本发明，进一步地，所述步骤s3中，采用经过预处理的pjfnk方法对离散化的方程组进行求解，其中预处理的方法为单矩阵预处理方法smp、对角元预处理方法bdp、有限差分预处理方法fdp、基于物理模型的预处理方法pbp中的任意一种。

100、为了更好地实现本发明，进一步地，所述步骤s3中，采用经过对角元预处理方法bdp预处理的pjfnk方法，对离散化的方程组进行求解；包括以下步骤：

101、步骤b1：采用差分格式构造jacobian近似向量：

102、离散后的方程组可以统一写为 f( u)=0，其中 u是解向量，是每一个节点上待求解的物理值；

103、针对非线性方程组，首先使用牛顿法对方程组进行线性化处理：

104、（15）

105、其中：

106、

107、

108、其中， j( u k)为非线性方程的jacobian矩阵， f( u)为关于 u的离散后的线性方程组；

109、 δ k+1为第 k+1次迭代得到的近似解 u k+1与第 k次迭代所得到的解 u k的差值；

110、步骤b2：采用gmres算法求解离散后的线性方程组，更新解向量 u k：

111、构造jacobian矩阵与矢量的乘积，以构造基向量，采用差分近似求解得到：

112、（16）

113、其中，为扰动参数，

114、 δ k为求解线性方程组时的残差，根据计算获得的 δ k更新 u k+1，其中 u k+1= u k+ δ k，更新为下一次迭代的解向量的值；

115、判断 δ k的数值是否低于预设数值，若不满足，则进入步骤b2，否则进入步骤s4；

116、若完成收敛，则解向量 u k的值即为各个节点上待求物理解的值，便可获得当前时间步下的流场状态。

117、为了更好地实现本发明，进一步地，所述步骤s4中，将单个颗粒的范围扩大，构成虚拟流体区域，在区域内，计算每一个有限元网格节点与颗粒中心位置的距离，将距离带入核函数，计算得到颗粒与流体之间相互作用力的权重系数；再根据颗粒的速度、形状与邻近颗粒的流场信息计算颗粒与流体的相互作用力。

118、为了更好地实现本发明，进一步地，所述步骤s4中，根据颗粒中心位置与邻近网格节点的相对距离，并通过高斯函数计算颗粒与流体之间相互作用力的权重系数：

119、（17）

120、其中： r i为颗粒中心位置，

121、 r j为收颗粒影响的流体域中的网格节点，

122、 k为颗粒半径放大倍数，

123、 d i为颗粒直径，

124、d()为高斯核函数。

125、为了更好地实现本发明，进一步地，所述步骤s4中，颗粒的所受流场作用力的合力计算如下：

126、（18）

127、其中： ε j为距离权重，

128、 v cell为有限元单元的体积，作为流场计算的动量源项；

129、 f i为第i个网格所施加给颗粒的作用力。

130、为了更好地实现本发明，进一步地，所述步骤s5中，基于dem离散元方法，实现颗粒与其他相邻颗粒的碰撞检测；根据接触模型计算颗粒之间的相互作用力与合力矩，结合流场作用力的合力，对颗粒的状态与位置显示求解。

131、本发明的有益效果如下：

132、（1）传统有限元方法求解采用牛顿法，或者拟牛顿法，求解较为复杂，尤其在与dem耦合层面，难以应用。本发明采用pjfnk方法求解流场的cfd-dem耦合，改进了求解层面以及建模层面。本发明可以提高cfd-dem耦合计算精度，节约计算内存，并具备优良的扩展性与适应性，适用于电池喷热失控喷发颗粒物在复杂的包内环境中运动特性模拟。

133、（2）采用本发明提供的方法，可以避免计算复杂的jacobian矩阵，节约计算内存，并易于并行计算。采用本发明提供的方法可以更精确地计算复杂流场，更准确地模拟流场中的颗粒物行为，更好地模拟电池热失控所喷发的颗粒物在复杂的电池包内流场中的碰撞、运动，从而更精准地分析颗粒物的运动特性与沉积特性。

134、（3）本发明通过有限元方法对流场ns偏微分方程组进行离散，再使用pjfnk方法进行求解，最后采用半解析的dem和cfd耦合方法，具有良好的扩展性，可以对复杂流场进行高精度数值计算，较为准确地模拟流场中的颗粒物行为，可以较好的适用于电池热失控所喷发的颗粒物在复杂的电池包内流场中的碰撞、运动的模拟，分析颗粒物的运动特性与沉积特性。

135、（4）本发明是基于有限元方法的cfd-dem耦合方案，有限元法则可以通过提高基函数插值精度达到任意阶数值模拟精度，且几何适应好，网格划分较为简单。有限元方法在偏微分方程求解中的应用十分广泛，可以基于fem实现统一的多物理场耦合计算框架，考虑物理场的耦合作用。有限元方法在偏微分方程应用领域较广，更容易扩展新的物理场，例如，在流场计算中，引入结构场的计算，实现热-流-颗粒-固耦合，或者引入电磁力，考虑带电颗粒在电场影响下与流体的相互作用，扩展能力大幅度增强。

136、（5）本发明能够克服一般有限元方法中，需要求解各单元上的jacbobian矩阵，耗时且求解难度极大的问题，规避直接求解jacobian矩阵，使得求解更加稳定、鲁棒性更好，也更易于并行计算，是求解层面的创新。