一种流化床中非球形粘性颗粒流动及碰撞特性获取方法与流程

1.本发明涉及一种基于cfd-dem模型的流化床中非球形粘性颗粒流动及碰撞特性获取方法，对颗粒在气固两相流过程中的流动及碰撞过程进行数值模拟和计算，可用于获取流化床内非球形粘性颗粒气固两相流动过程的流动及碰撞特性。


背景技术：

2.粘性颗粒流动过程是一种工业生产中常见的过程，目前在造粒、燃烧、除尘等领域均有较为广泛的应用，在颗粒气固两相流动过程中常伴随颗粒间的碰撞过程。数值模拟方法是研究颗粒流动特性的一种重要方法，然而当前大多数研究都将颗粒形状简化为了球形。虽然将颗粒假设为球形可以缩短模拟时间，提高计算效率，但是在真实的工业过程中，颗粒通常为外观形貌不规则的非球形颗粒，非球形颗粒具有较为明显的各向异性特征，其流动及碰撞过程与球形颗粒的流动过程有显著的差别，如果不考虑颗粒形状的影响，将会使数值模拟计算结果的准确度大幅降低，因此对非球形颗粒流动特性的数值模拟研究具有重要意义。


技术实现要素：

3.针对上述情况，为克服现有技术之不足，本发明之目的就是提供一种流化床中非球形粘性颗粒流动及碰撞特性获取方法。
4.本发明解决的技术方案是：
5.一种流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，包括以下步骤：
6.(1)采用球元拼接法建立不同球形度的非球形颗粒物理模型；
7.(2)建立用于气固两相流动过程的流化床物理模型；
8.(3)建立流化床气固两相流动过程中的气相运动方程组；
9.(4)建立流化床气固两相流动过程中的固相运动方程组；
10.(5)建立气固两相流动过程中颗粒间的粘附性接触模型；
11.(6)确定气固两相流动过程中颗粒碰撞存在的情况，具体包括：
12.a、单个球元与单个球元间的碰撞；
13.b、单个球元与多个球元碰撞；
14.c、多个球元与多个球元碰撞；
15.d、单个球元与壁面碰撞；
16.e、多个球元与壁面碰撞；
17.(7)建立气固两相流动过程的非球形颗粒碰撞核模型；
18.(8)基于上述步骤建立的模型和运动方程组，在相同的基本求解参数下，对流化床内不同形状的非球形颗粒气固两相流动过程进行cfd-dem仿真计算；比较不同球形度颗粒在流化床模型内的运动轨迹及碰撞率，得到颗粒形状对流化床内粘性颗粒流动及碰撞特性
的影响规律；
19.所述基本求解参数包括气相流速、颗粒密度、颗粒数量。
20.与现有技术相比，本发明具有以下优点：
21.(1)本发明采用dem离散单元模型对颗粒相进行理论建模，将气相视作连续相，固相视作离散相，可对气固两相流动过程中任意时间段、任意形状颗粒的流动及碰撞过程进行追踪，便于分析粘性非球形颗粒气固两相流动过程的流动及碰撞特性。
22.(2)本发明综合考虑非球形颗粒的形状以及颗粒相与气相、颗粒相与颗粒相之间的相互作用，使仿真计算结果更贴近实际过程，并更具可靠性。
23.(3)本发明可以通过更改入口参数及边界条件，对比不同工况下非球形颗粒的流动特性，从而对含有粘性颗粒气固两相流动过程的工艺设备进行优化。
24.(4)本发明具有仿真结果直观、计算精度高、可靠性好等显著优点，可以较为便捷地得到实验方法中难以直接获取的颗粒流动及碰撞的规律及参数。
附图说明
25.图1为本发明实施例流化床物理模型的网格划分图。
26.图2为本发明实施例非球形颗粒模型示意图。
27.图3为本发明实施例不同形状颗粒轨迹图。
28.图4位本发明实施例非球形颗粒碰撞率随粒径变化图。
29.图5位本发明实施例碰撞率随颗粒球形度变化图。
30.图6为本发明cfd中的两相流模型耦合计算流程图。
具体实施方式
31.以下结合附图和实施例对本发明的具体实施方式作进一步详细说明。
32.本发明公开了一种流化床中非球形粘性颗粒流动及碰撞特性获取方法，包括以下步骤：
33.(1)采用球元拼接法建立不同球形度的非球形颗粒物理模型；
34.具体为：
35.球元拼接法即采用若干个小球面通过拼接的方式来构建不同形状颗粒的方法；采用球形度来描述非球形颗粒的外观形貌特征，其数学表达式如下：
[0036][0037]
其中，s
ev
为与颗粒体积相同的球体表面积，sc为颗粒的实际表面积，d
ev
为与颗粒体积相同的球体直径，d
es
为与颗粒表面积相同的球体直径；
[0038]
(2)建立用于气固两相流动过程的流化床物理模型；
[0039]
(3)建立流化床气固两相流动过程中的气相运动方程组；
[0040]
气相运动方程组为：
[0041]
ug|时，则发生相对滑移，切向力的表达式改写为：
[0063][0064]
其中，μ为颗粒间的摩擦系数；
[0065]
进一步地，颗粒所受的法向力与颗粒的接触半径和表面能间的关系用下式表达：
[0066][0067]
其中，ri为颗粒的有效接触半径；r
*
为颗粒的有效半径；y
*
为颗粒的有效杨氏模量系数；γi为颗粒的表面能。
[0068]
(6)确定气固两相流动过程中颗粒碰撞存在的情况，具体包括：
[0069]
a、单个球元与单个球元间的碰撞；
[0070]
b、单个球元与多个球元碰撞；
[0071]
c、多个球元与多个球元碰撞；
[0072]
d、单个球元与壁面碰撞；
[0073]
e、多个球元与壁面碰撞；
[0074]
虽然球形颗粒与非球形颗粒的碰撞方式不同，但本质上依然是球面之间的碰撞，即非球形颗粒所受的碰撞力为各球元所受碰撞力的合力；
[0075]
(7)建立气固两相流动过程的非球形颗粒碰撞核模型；即计算非球形颗粒间的碰撞率c
ij
，通过下式计算：
[0076][0077]
其中，c
ij
为颗粒i与颗粒j的碰撞率常数，为颗粒i与颗粒j的碰撞次数，ni为颗粒i的数量，nj为颗粒j的数量，v
all
为床体的体积，t
sim
为仿真计算的时间。
[0078]
(8)基于上述步骤建立的模型和运动方程组，在相同的基本求解参数下，对流化床内不同形状的非球形颗粒气固两相流动过程进行cfd-dem仿真计算；比较不同球形度颗粒在流化床模型内的运动轨迹及碰撞率，得到颗粒形状对流化床内粘性颗粒流动及碰撞特性的影响规律；
[0079]
所述基本求解参数包括气相流速、颗粒密度、颗粒数量。
[0080]
其中，cfd-dem仿真计算方法是一种基于欧拉-拉格朗日计算模型的方法，其基本求解思路是：dem模型在数值模拟过程中将气相看作连续介质，使用欧拉方法对其进行描述，通过求解n-s方程得出气相的运动特性；将固相看作离散相，使用拉格朗日方法进行描述，对每个颗粒的运动轨迹和运动特性进行追踪。球形度为与颗粒体积相同的球体表面积和颗粒的实际表面积之比，因此值为1时，颗粒为球形颗粒；值不为1时，颗粒为非球形颗粒。
[0081]
本发明经实际应用，取得了良好的技术效果，实施例如下，包括如下步骤：
[0082]
s1：如图1所示，在本实施例中，所述物理模型为流化床物理模型，内部有扰流叶片，扰流叶片的作用是加强流场的湍流扰动和颗粒间的碰撞频率。所述流化床模型的长、
宽、高分别为2.5mm、0.5mm和0.15mm；采用正六面体结构化网格对所述流化床模型进行划分，得到63601个网格单元；气相从左侧边界进入聚并装置，边界条件采用速度入口；模型右侧边界为气相和固相出口，边界条件采用压力出口，扰流挡板和壁面均设置为刚性壁面。
[0083]
s2：在本实施例中，在颗粒工厂中分别设置球形度值为0.6、0.67、0.79、0.81、0.85、0.88的非球形颗粒，如图2所示。为方便举例，同时设置球形度为1的球形颗粒；为比较颗粒大小对颗粒碰撞特性的影响，本实施例设置尺寸大小不同但球形度值为0.6的颗粒。
[0084]
s3：在实施例中，基本求解参数为：气相粘度1.8335
×
10-5
pa
·
s，气相流速10m/s，颗粒密度800kg/m3，dem时间步长为5
×
10-8s，cfd时间步长为5
×
10-6s。
[0085]
s4：在实施例中，在cfd中建立气相运动方程组：
[0086][0087][0088]
其中，ε为网格内空隙率；ρg为气相密度；ug为气相速度；p为气相压力；g为当地重力加速度；τ为应力张量。fg为单位网格内气固作用力的矢量和：fg＝fd+fm+fs，fd为曳力，fm为magnus升力，fs为saffman升力，由于升力数量级远小于曳力，可忽略不计；所述fd由下式计算：
[0089]
fd＝β(u
s-ug)
[0090]
其中，β为气固相间曳力系数，由下式计算：
[0091][0092]
其中，μg为气相黏度；为颗粒的当量直径；cd为单个颗粒的曳力系数，由下式计算：
[0093][0094]
颗粒雷诺数re可用下式表示：
[0095][0096]
其中，ug为气相速度，us为固相速度。
[0097]
s5：在实施例中，在dem中建立固相运动方程组，具体为：
[0098][0099]
[0100]
式中：ms为颗粒质量；fc为颗粒间接触力；gs为重力；is为转动惯量；ωs为旋转角速度；ts为力矩。
[0101]
s6：在实施例中，建立颗粒粘附性接触模型，接触力fc分为切向力f
ct
和法向力f
cn
，表达式如下：
[0102]fct
＝-k
tdt-η
t
[(v
r-vr·
n)n+rωs×
n]
[0103]fcn
＝-k
ndn-ηn(vr·
n)n
[0104]
其中，k
t
、kn分别为切向、法向弹性系数；d
t
、dn分别表示切向、法向的相对位移；μ
t
、μn分别为切向、法向的阻尼系数；ωs为颗粒相对角速度；vr为颗粒相对线速度。当f
ct
＞μ|u
s-ug|时，则发生相对滑移，切向力的表达式改写为：
[0105][0106]
其中，μ为颗粒间的摩擦系数。
[0107]
颗粒所受的法向力与颗粒的接触半径和表面能间的关系用下式表达：
[0108][0109]
其中，ri为颗粒的有效接触半径；r
*
为颗粒的有效半径；y
*
为颗粒的有效杨氏模量系数；γi为颗粒的表面能。
[0110]
s7：分析气固两相流过程中颗粒碰撞存在的情况：(1)单个球元与单个球元间的碰撞；(2)单个球元与多个球元碰撞；(3)多个球元与多个球元碰撞；(4)单个球元与壁面碰撞；(5)多个球元与壁面碰撞。虽然球形颗粒与非球形颗粒的碰撞方式不同，但本质上依然是球面之间的碰撞，即非球形颗粒所受的碰撞力为各球元所受碰撞力的合力。
[0111]
s8：在实施例中，在dem中建立非球形颗粒碰撞核模型，非球形颗粒间的碰撞率c
ij
由下式计算：
[0112][0113]
其中，c
ij
为颗粒i与颗粒j的碰撞率常数，为颗粒i与颗粒j的碰撞次数，ni为颗粒i的数量，nj为颗粒j的数量，v
all
为床体的体积，t
sim
为仿真计算的时间。
[0114]
在实施例中，dem和cfd是通过cfd中的两相流模型实现并行计算的，耦合计算流程如图6所示，cfd软件通过计算每一个时间步长内的流场信息，对流场进行求解至收敛后，在dem中进行迭代，然后通过耦合接口将每一个颗粒的位置、运动、体积等信息传递至cfd中，计算气固相间作用，待收敛后固相与气相的相互作用将再次通过耦合接口传递至dem中作为颗粒体积力影响颗粒的运动，而对流体的作用则通过动量源相的方式作用于流体中。逐次循环迭代，实现全过程的瞬态模拟。
[0115]
选取不同球形度的颗粒，其球形度分别为0.6、0.67、0.79、0.81、0.85、0.88，1按照本发明的方法步骤，对上述7种颗粒在流化床中的流动过程进行仿真，得出每组颗粒在流化床中的运动轨迹，如图3所示，从图中可以看出，当颗粒经过聚并叶片后球形颗粒的轨迹相对较平稳，而非球形颗粒的轨迹较紊乱。并且球形度越接近于1的非球形颗粒，运动轨迹也
越平稳，而非球形颗粒的球形度越小，运动轨迹越紊乱，发生碰撞的概率也越大。产生这种现象的主要原因是球形颗粒的形状比较规则，在颗粒各个方向的受力较为均匀，更容易达到纵向的受力平衡从而接近于直线运动。而非球形颗粒由于其形状不规则，各向异性特征明显，在流动过程中具有较强的取向性，在颗粒的各个表面上受到的流场力及颗粒间接触力的大小和方向均不相同，因此在流场中的运动轨迹更容易发生改变。
[0116]
如图4为形状相同大小不同的颗粒碰撞率分布图，粒径范围为30～50μm。从图中可以看出，当颗粒大小相等时，小颗粒间的碰撞率较小，而大颗粒间的碰撞率较大，主要是因为颗粒发生碰撞的概率与颗粒的表面积和体积正相关；对比不同粒径颗粒间的碰撞率可以看出，对于一种颗粒而言，与其发生碰撞的颗粒越大，碰撞率就越大，这是因为大颗粒的表面积与体积较大，而小颗粒的表面积与体积较小，在流场中两者由于所受气相作用力的大小差距较大，容易产生速度差，因此两者发生接触的概率较大，并且大颗粒对小颗粒在范德华粘附力的作用下具有较强的吸附作用，因此两者的碰撞率较大；从图中可以看出，相同粒径颗粒之间碰撞率较小，这是由于相同粒径的颗粒在流场中的受力情况较为接近，曳力与升力的大小和方向较为接近，其运动轨迹也较为相似，不容易产生较大的速度差，因此发生碰撞的概率也较小。
[0117]
如图5所示，为不同球形度颗粒碰撞率分布图，反映了流化床中不同形状的非球形颗粒碰撞率的大小。从图中可以看出，球形颗粒与球形颗粒之间的碰撞率最小，这是由于球形颗粒的形状比较规则，在流场中受力分布相对于非球形颗粒来说较为均匀，受流场湍流扰动作用较小，所以颗粒间的碰撞率也较小；对于一种球形度的非球形颗粒，与其发生碰撞的颗粒球形度越小，碰撞率就越大，这是因为颗粒的球形度越小，其表面越不规则，在流场中的受力情况就越复杂，运动轨迹越容易发生改变，因此相互之间发生碰撞的概率也就越大；当颗粒球形度和与其发生碰撞的颗粒球形度相同时，颗粒间的碰撞率较小，这是因为其外观形貌特征相同，在流场中受力情况比较类似，因此其运动轨迹不容易发生交叉，导致碰撞率较小。
[0118]
有上述情况可以清楚的看出：
[0119]
(1)本发明采用dem离散单元模型对颗粒相进行理论建模，将气相视作连续相，固相视作离散相，可对气固两相流动过程中任意时间段、任意形状颗粒的流动及碰撞过程进行追踪，便于分析粘性非球形颗粒气固两相流动过程的流动及碰撞特性。
[0120]
(2)本发明综合考虑非球形颗粒的形状以及颗粒相与气相、颗粒相与颗粒相之间的相互作用，使仿真计算结果更贴近实际过程，并更具可靠性。
[0121]
(3)本发明可以通过更改入口参数及边界条件，对比不同工况下非球形颗粒的流动特性，从而对含有粘性颗粒气固两相流动过程的工艺设备进行优化。
[0122]
(4)本发明具有仿真结果直观、计算精度高、可靠性好等显著优点，可以较为便捷地得到实验方法中难以直接获取的颗粒流动及碰撞的规律及参数。

技术特征：
1.一种流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，包括以下步骤：(1)采用球元拼接法建立不同球形度的非球形颗粒物理模型；(2)建立用于气固两相流动过程的流化床物理模型；(3)建立流化床气固两相流动过程中的气相运动方程组；(4)建立流化床气固两相流动过程中的固相运动方程组；(5)建立气固两相流动过程中颗粒间的粘附性接触模型；(6)确定气固两相流动过程中颗粒碰撞存在的情况，具体包括：a、单个球元与单个球元间的碰撞；b、单个球元与多个球元碰撞；c、多个球元与多个球元碰撞；d、单个球元与壁面碰撞；e、多个球元与壁面碰撞；(7)建立气固两相流动过程的非球形颗粒碰撞核模型；(8)基于上述步骤建立的模型和运动方程组，在相同的基本求解参数下，对流化床内不同形状的非球形颗粒气固两相流动过程进行cfd-dem仿真计算；比较不同球形度颗粒在流化床模型内的运动轨迹及碰撞率，得到颗粒形状对流化床内粘性颗粒流动及碰撞特性的影响规律；所述基本求解参数包括气相流速、颗粒密度、颗粒数量。2.根据权利要求1所述的流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，所述步骤(1)采用球元拼接法建立不同球形度的非球形颗粒物理模型具体为：采用球形度来描述非球形颗粒的外观形貌特征，其数学表达式如下：其中，s
ev
为与颗粒体积相同的球体表面积，s
c
为颗粒的实际表面积，d
ev
为与颗粒体积相同的球体直径，d
es
为与颗粒表面积相同的球体直径。3.根据权利要求1所述的流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，所述步骤(3)气相运动方程组为：所述步骤(3)气相运动方程组为：其中，ε为网格内空隙率；ρ
g
为气相密度；u
g
为气相速度；p为气相压力；g为当地重力加速度；τ为应力张量；f
g
为单位网格内气固作用力的矢量和：f
g
＝f
d
+f
m
+f
s
，f
d
为曳力，f
m
为magnus升力，f
s
为saffman升力；所述f
d
由下式计算：f
d
＝β(u
s-u
g
)其中，β为气固相间曳力系数，具体表达式如下：
其中，μ
g
为气相黏度；为颗粒的当量直径；c
d
为单个颗粒的曳力系数，通过以下公式求得：颗粒雷诺数re用下式表示：其中，u
g
为气相速度，u
s
为固相速度。4.根据权利要求1所述的流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，所述步骤(4)固相运动方程组可通过下式描述：所述步骤(4)固相运动方程组可通过下式描述：式中：m
s
为颗粒质量；f
c
为颗粒间接触力；g
s
为重力；i
s
为转动惯量；ω
s
为旋转角速度；t
s
为力矩。5.根据权利要求1所述的流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，所述步骤(5)建立气固两相流动过程中颗粒间的粘附性接触模型具体为：在气固两相流动过程中，粘性颗粒间存在接触力，进一步地，建立颗粒粘附性接触模型，接触力f
c
分为切向力f
ct
和法向力f
cn
，可由下式来表达：f
ct
＝-k
t
d
t-η
t
[(v
r-v
r
·
n)n+rω
s
×
n]f
cn
＝-k
n
d
n-η
n
(v
r
·
n)n其中，k
t
、k
n
分别为切向、法向弹性系数；d
t
、d
n
分别表示切向、法向的相对位移；μ
t
、μ
n
分别为切向、法向的阻尼系数；ω
s
为颗粒相对角速度；v
r
为颗粒相对线速度。当f
ct
＞μu
s-u
g
时，则发生相对滑移，切向力的表达式改写为：其中，μ为颗粒间的摩擦系数；进一步地，颗粒所受的法向力与颗粒的接触半径和表面能间的关系用下式表达：其中，r
i
为颗粒的有效接触半径；r
*
为颗粒的有效半径；y
*
为颗粒的有效杨氏模量系数；
γ
i
为颗粒的表面能。6.根据权利要求1所述的流化床中非球形粘性颗粒流动及碰撞特性获取方法，其特征在于，所述步骤(7)建立气固两相流动过程的非球形颗粒碰撞核模型，即计算非球形颗粒间的碰撞率c
ij
，通过下式计算：其中，c
ij
为颗粒i与颗粒j的碰撞率常数，为颗粒i与颗粒j的碰撞次数，n
i
为颗粒i的数量，n
j
为颗粒j的数量，v
all
为床体的体积，t
sim
为仿真计算的时间。

技术总结
本发明涉及一种流化床中非球形粘性颗粒流动及碰撞特性获取方法，本发明采用DEM离散单元模型对颗粒相进行理论建模，将气相视作连续相，固相视作离散相，可对气固两相流动过程中任意时间段、任意形状颗粒的流动及碰撞过程进行追踪，便于分析粘性非球形颗粒气固两相流动过程的流动及碰撞特性；综合考虑非球形颗粒的形状以及颗粒相与气相、颗粒相与颗粒相之间的相互作用，使仿真计算结果更贴近实际过程，并更具可靠性；通过更改入口参数及边界条件，对比不同工况下非球形颗粒的流动特性，从而对含有粘性颗粒气固两相流动过程的工艺设备进行优化；具有仿真结果直观、计算精度高、可靠性好等显著优点。好等显著优点。好等显著优点。


技术研发人员：王子涵 官卜瑞 耿彪 李冰天 原峥
受保护的技术使用者：河南京能滑州热电有限责任公司
技术研发日：2023.04.24
技术公布日：2023/9/20