Fluent的气固两相流中离散颗粒的数值模拟

法向/切向弹性系数(N/m)
800
法向/切向恢复系数
0.9/0.35
滑动摩擦系数
0.3
其他参数
网格尺寸(mm）
10(x-direction)×20(y-direction)
时间步长(s)
0.00001
颗粒数
节涌(1440)/鼓泡(1224)
图 2 展示了不同时刻鼓泡床颗粒速度云图。系统内床层主要经历两个阶段，即启动阶段和稳定流化阶段。在
u
f
)
+
ρf ε
( ∂ε ∂t
+uf
⋅ ∇ε ) = 0
(8)
或
∂ρ f ∂t
+ ∇ ⋅ (ρ f u f ) = Sc
(9)
其中：
Sc = − ρ f ε
∂ε (
∂t
+uf
⋅ ∇ε )
动量方程：
ε
∂ ∂t
(ρ
f
u
f
)
+
ε∇
⋅
(ρ
f
u
f
u
f
)
+
ρ
f
u
f
(
∂ε ∂t
+uf
⋅ ∇ε ) = −ε∇p + ∇ ⋅ (ε τ f
颗粒群中颗粒 Reynolds 为：
ε Re =
ρf u
−υ dp
μf
式(12)、相间动量交换系数β及颗粒 Reynolds 计算式中的ε、ν均为网格内的平均值。
流场求解后，在式(12)基础上，根据牛顿第三定律，单颗粒所受曳力为：
Fd
= Vpβk 1−εk
(uk
−υk )
(13)
式中εk、βk 分别代表颗粒 k 的局部孔隙率、局部动量交换系数；uk 为颗粒 k 质心处的虚拟连续相速度。
参考文献
[1] 欧阳洁，孙国刚，李静海. 气固两相流模拟的随机离散模型. 数值计算与计算机应用[J], 2003, 2: 88-100 [2] Tsuji Y, Tannka T, Ishida T. Lagrangian numerical-simulation of plug flow of cohesionless particles in a horizontal pipe[J]. Powder Technology, 1992, 71:239-250 [3] Xu B H, Yu A B. Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics[J]. Chemical Engineering Science, 1997, 52:2785-2809 [4] T. Kawaguchi,T. Tanaka, Y.Tsuji. Numerical simulation of two-dimensional fluidized beds using the discrete element method(comparison between the two- and three-dimensional models).Powder Technology, 1998,96: 129-138 [5] C.N.Wu, Y.Cheng, Y.L.Ding, Y. Jin. CFD-DEM simulation of gas-solid reacting flows in fluid catalytic cracking(FCC) process. Chemical Engineering Science, 2010, 65:542-549 [6] Tsuji Y, Tanaka T, Ishida T. Lagrangian numerical simulation of plug flow of cohesionless particle in a horizontal pipe. Powder Technology, 1992, 71: 239-250 [7] Watano S, Saito S, Suzuki T. Numerical simulation of electrostatic in powder pneumatic conveying process. Powder Technology, 2003, 135/136: 112-117 [8] Li J, Webb C, Pandiella S S, et al. Solid deposition in low-velocity slug flow pneumatic conveying. Chemical Engineering and Processing, 2005,44: 167-173 [9] 李云，顾兆林，王梓入. 沉浸管流化床内气固两相流动数值模拟. 工程热物理年报. 中国工程热物理年会 多相流 编号： 066073 [10] C.L.Wu, J.M.Zhan, Y.S.Li, K.S.Lam. Dense particulate flow model on unstructured mesh . Chemical Engineering Journal, 2006, 61: 5726-5741 [11] C.L.Wu, A.S.Berrouk, K. Nandakumar. Three-dimensional discrete particle model for gas-solid fluidized beds on unstructured mesh . Chemical Engineering Journal, 2009, 152: 514-529 [12] Alberto Di Renzo, Francesco Paolo Di Maio. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes. Chemical Engineering Journal, 2004, 59: 525-541 [13] D.Gidaspow. Multiphase flow and fluidization-continuum and kinetic theory descriptions. Academic Press,San Diego,1994
Vi
= Vi,0
+
dVi dt
Δt
(3)
ωi
= ωi,0
+
dωi dt
Δt
(4)
xyzi = xyzi,0 + Vi Δt
(5)
2.2 连续相控制方程
连续相控制方程采用基于局部平均并考虑流-固耦合作用的双流体方程(Model A)[13]，连续性方程及动量方程
如下：
∂(ερ ∂t
f
)
+
∇
⋅
(ερ
f
u
关键词：FLUENT；离散单元法；计算流体力学
1、引言
对气固流动系统，双流体模型以其大规模模拟的可行性在数值模拟领域占据重要地位，但双流体模型中的颗 粒连续假设从本质上削弱了颗粒流体系统中非均匀特性的描述[1]。与传统两相流数值模拟方法相比，离散单元法 (DEM，Discrete/Distinct Element Method，又称软球模型)可跟踪每个颗粒的运动信息，并考虑颗粒-颗粒、颗 粒-壁面碰撞及颗粒-流体相互作用，在气固非均匀流数值模拟方面具有强大优势。1992 年日本 Tsuji 首次将软球 模型用于流化床模拟，得到与试验一致的最小鼓泡速度、压力波动频率及不同气速下鼓泡床内的颗粒运动信息[2]。 而后众多学者采用该方法对气固流动行为进行了详细研究， XU 和 YU 运用软球模型模拟了流化床中的鼓泡和节涌 现象 [3]；T. Kawaguchi 比较了二维与三维软球模型的鼓泡床中，气速、隔板及摩擦系数对颗粒动态行为的影响[4]； C.N. WU 研究了考虑颗粒-颗粒及颗粒-壁面传热的流化床反应器中颗粒与气体行为[5]；Tsuji[6]、Watano[7]、Li [8] 则将该方法应用于气粒输送领域。至今，离散单元法已被广泛用于气固流动系统，对深入了解气固流动行为具有 重要意义
2、数学模型
2.1 颗粒相控制方程 在气(液)固系统中，颗粒主要受四种力作用—气体对颗粒的曳力、压力梯度力、颗粒重力及颗粒-颗粒间碰
撞力。根据牛顿第二定律，颗粒运动方程可表述为：
mi
dVi dt
= Fd ,i
+ mi g + Σ(Fn,ij
+ Ft,ij ) − V p∇p
(1)
Ii
dωi dt
= ΣTt,ij
3、模拟计算与分析
模拟对象为伪三维流化床系统，即流化床厚度为一个颗粒直径，孔隙率计算沿用 XU 和 YU 的准三维方法[3]。 方程离散采用二阶迎风格式，流化床进、出口分别采用速度入口、压力出口。计算过程中，不考虑前后壁面对颗 粒的作用力;模拟中气体、颗粒主要物性参数见表 1。
发生于小直径床的节涌现象可用软球模型进行模拟，如图 1 所示。在活塞流流化床中，颗粒所受的曳力强烈 依赖于孔隙率，孔隙率较低，颗粒所受曳力越大。随着床层的上升，颗粒群底部空隙率降低，该区域颗粒所受到 的曳力逐渐减小，导致床层变形并且颗粒向下移动。
中国工程热物理学会
传热传质学
学术会议论文
编号：1 2 3 1 1 8
基于 Fluent 的气固两相流中离散颗粒的数值模拟
任立波 1 何海澜 2 韩吉田 1
1、 山东大学能源与动力工程学院，山东济南，250061 2、兰州兰石换热设备有限责任公司，甘肃兰州，730050 Tel：0531-88399060 Email：jthan@sdu.edu.cn
连续相对颗粒作用力包括曳力、压力梯度力、浮力、附加质量力等，本文只考虑曳力及压力梯度力，其中压
力梯度力已在式(1)中给出。本算法基于双流体模型气相方程中的耦合项推广到颗粒轨道模型，再反作用于颗粒相。
流体运动方程中的动量交换源项取为：
S fp = β (u − υ )
(12)
其中相间动量交换系数β为：
本文基于通用计算流体软件 FlUENT，通过用户自定义函数(User Defined Function) 实现了伪三维 DEM-CFD 耦合算法，其中连续相控制方程由通用计算流体软件 FLUENT 求解。构建空隙率标量场，重组基于局部平均并考虑 流-固耦合作用的连续相控制方程，提高了连续相控制方程的求解稳定性。将该算法用于球形颗粒的气固鼓泡床及 节涌工况模拟中，数值模拟计算结果与有关数值模拟及试验结果定性一致。
0.1s
0.2s
0.3s
0.4s
0.5s
0.6s
0.8s
图 1 节涌现象颗粒速度云图
0s
0.1s
0.2s
0.3s
0.4s
0.5s
0.6s
0.7s
关键词：FLUENT；离散单元法；计算流体力学
图 2 鼓泡现象颗粒速度云图
基于 Fluent 的气固两相流中离散颗粒的数值模拟
任立波 1 何海澜 2 韩吉田 1
1、 山东大学能源与动力工程学院，山东济南，250061 2、兰州兰石换热设备有限责任公司，甘肃兰州，730050 (TEL: 13173022361;Email: jthan@sdu.edu.cn)
摘要：基于通用计算流体软件 FlUENT，通过用户自定义函数(User Defined Function) 实现了伪三维 DEM-CFD 耦合算法，其中连续 相控制方程由通用计算流体软件 FLUENT 求解。构建空隙率标量场，重组基于局部平均并考虑流-固耦合作用的连续相控制方程，提 高了连续相控制方程的求解稳定性。最后，将该算法用于球形颗粒的气固两相流数值模拟中，数值模拟计算结果与有关数值模拟及 试验结果定性一致。
对气固流动系统，双流体模型以其大规模模拟的可行性在数值模拟领域占据重要地位，但双流体模型中的颗 粒连续假设从本质上削弱了颗粒流体系统中非均匀特性的描述[1]。与传统两相流数值模拟方法相比，离散单元法 (DEM，Discrete/Distinct Element Method，又称软球模型)可跟踪每个颗粒的运动信息，并考虑颗粒-颗粒、颗 粒-壁面碰撞及颗粒-流体相互作用，在气固非均匀流数值模拟方面具有强大优势。
)+Sf
+ ερ f g
(10)
或
∂(ρ f u ∂t
f
)
+
∇ ⋅(ρ
f
uf
uf
)
=
−∇p
+
∇ ⋅ (τ
f
)
+
Sm
+
ρ
f
g
(11)
其中：
Sm
=
S fp ε
−
ρf g ε
( ∂ε ∂t
+uf
⋅ ∇ε ) + τ
⋅ ∇ε ε
或
Sm
= Scu f
+
1 ε
(S
fp
+τ ⋅ ∇ε )
2.3 连续相-颗粒相耦合模型
泡床模拟中，数值模拟计算结果与有关数值模拟及试验结果定性一致。
至今，离散单元法已被广泛用于气固流动系统，对深入了解气固流动行为具有重要意义。 DEM 算法中颗粒-
颗粒、颗粒-壁面搜索、碰撞算法计算量非常庞大，并且计算时间步长较小，这就限制了其大规模应用。随计算机
硬件水平的不断提高，开发高效率 DEM 并行算法为将会成为离散颗粒数值模拟领域的当务之急。
通用计算流体软件 FlUENT 可作为确定性颗粒轨道模型求解器[9-11]，具有通用性及稳定性好、计算精度高、不 会出现代码理解偏差导致的错误等优点。本文将 DEM 与计算流体力学软件 FlUENT 相结合，基于颗粒尺度碰撞、摩 擦等作用和介观网格尺度的颗粒-流体相互作用，实现伪三维 CFD-DEM 耦合算法，其中 DEM 算法通过 UDF(User Defined Function)实现。
启动阶段，入口附近形成一个椭圆形气泡，气泡向上运动，床层上部演变为活塞流。被抬起的颗粒回落床层，而
后形成稳定的流化态，本文结果与 Tsuji et al[4], Xud Yu[3]试验及模拟结果定性一致。
5 结论
基于通用计算流体软件 FlUENT，通过用户自定义函数(User Defined Function) 实现了伪三维 DEM-CFD 耦 合算法，其中连续相控制方程由通用计算流体软件 FLUENT 求解。构建空隙率标量场，重组基于局部平均并考虑流 -固耦合作用的连续相控制方程，提高了连续相控制方程的求解稳定性。最后，Fra Baidu bibliotek该算法用于球形颗粒的节涌及鼓
0.1s
0.2s
0.3s
0.4s
0.5s
0.6s
0.8s
图 1 节涌现象颗粒速度云图
0s
0.1s
0.2s
0.3
0.4s
0.5s
0.6s
0.7s
图 2 鼓泡现象颗粒速度云图
表 1 模拟主要参数
气体参数
密度(Kg/m3)
1.225
粘度(Pa·s)
1.7894e-05
颗粒参数
形状
球形
密度(Kg/m3)
2700
(2)
式中：mi 表示颗粒质量，Fn,ij、Ft,ij 分别表示颗粒 i 受到颗粒 j 的法向力和切向力；Ti 是由颗粒碰撞切向应力 产生的力矩；Ii、Vi、ωi 分别为颗粒 i 的转动惯量、平动速度及角速度；Vp、▽p 分别颗粒体积及其颗粒处的压力 梯度。颗粒间作用力(Fn,ij、Ft,ij) 采用线性模型[12]；壁面看成静止的颗粒，采用与颗粒间碰撞类似的公式计算。在 上述力和力矩作用下，颗粒的速度、旋转速度及位置变化可由式(3)-(5)表示：
f
)
=
0
(6)
∂(ερ f u ∂t
f
)
+
∇
⋅ (ερ
f
u
f
u
f
)
=
−ε∇p
+
∇ ⋅ (ε
τ
f
)
+
S
fp
+ ερ
f
g
(7)
式中 Sfp 表示流固耦合项。
为增加计算求解稳定性，重组控制方程(6)、(7) (Sc、S m 分别为重组后连续性方程及动量方程的源相)：
连续性方程：
∂ρ f ∂t
+
∇ ⋅(ρ
f
⎧ ⎪β ⎪ ⎨ ⎪ ⎪⎩
(1 − ε = 150 ε
)2
μg
d
2 p
+ 1.75(1 − ε
) ρg dp
u −υ
β
=
3 4
Cd0
ε
(1 − ε dp
)
ρf u
− υ ε −2.65
ε < 0.8 ε ≥ 0.8
[ ] Cd
=
⎪⎧ ⎨
24 Re
1 + 0.15(Re)0.687
⎪⎩
0.44
Re < 1000 Re ≥ 1000