搞懂这些，你就掌握了气固两相流传热模拟的核心

以上几个问题欢迎优秀的你留言回答，当然了，如果嫌弃留言区面积不够，欢迎投稿~

   

   

总之，以上这些都可以用三大传热途径来解释：辐射传热、导热传热和对流传热。上一篇简单介绍了如何模拟气固两相流及相关成果：

👇🏻

关于稠密气固两相流，你get到多少？

接下来介绍稠密气固两相流CFD-DEM模拟中的传热模型。

• 方程1

• 方程2

其中，流体速度、密度和压力分别为u_f，ρ_f 和p_f。重力加速度为g。ε_f为网格空隙率，亦即某个网格内流体所占据的体积分数。S为动量交换源项。牛顿流体粘性应力张量τ_f 的计算方法为：

• 方程3

其中，λ_f 为流体的体积粘度，而μ_f 为流体的剪切粘度，δ_ij 为Kronecker函数。

在稠密气固两相流动中，颗粒的运动遵循拉格朗日框架下的牛顿第二定律，采用积分方法对颗粒的运动方程进行显式推进，从而求解得到每个颗粒的加速度、平动速度、角速度和位置等信息。系统中的颗粒会受到流体曳力、压力梯度力、重力、Saffman升力和虚拟质量力（附加质量力）等，但是由于颗粒/气体密度比值较大（1000 ~ 2500），因此只需要考虑重力、曳力和压力梯度力，其余的力可以忽略不计。对于每一个特定的颗粒，其平动和转动控制方程分别为：

• 方程4
  

• 方程5

其中，m_p 是颗粒质量，v_p 是颗粒平动速度，ω_p是颗粒角速度，I_p 是转动惯量，大小为  。L是颗粒质心到接触点的距离。N是和当前颗粒碰撞的所有颗粒数量。此外，F_gp (=-v_p▽p_f )为压力梯度力，F_drag为曳力，F_c为碰撞接触力。

• 方程7
  

其中，对流换热系数通过无量纲的努塞尔数（Nu_p,i）进行计算，通过转换，表示为如下形式：

• 方程8

其中，κ_f 为流体热导率，d_p,i 为颗粒粒径。对于Nu_p,i 的计算，一般是将其同无量纲的颗粒雷诺数（Re_p,i）和普朗特数（Pr）相关联。本文采用Li和Mason^[1]的关联式进行计算：

• 方程9：
  

其中，普朗特数Pr= μ_fC_p,f /κ_f，C_p,f 为流体比热容。

通常而言，颗粒-颗粒导热在气固流化床中可以假设不计。但是，对于固定床和低速鼓泡床内，颗粒-颗粒碰撞较为剧烈，应当考虑颗粒-颗粒碰撞导热。其导热传热发生在两个相互碰撞颗粒较短时间内的较小碰撞面上。对于众多的颗粒-颗粒导热模型，采用广泛使用的Batchelor和O'Brien模型[2]，其写作：

• 方程10

其中，κ_p,i 和κ_p,j 是颗粒i 和颗粒j 的热导率。

最终的颗粒-颗粒导热传热量为：

除了接触时的直接导热外，颗粒和颗粒还会通过其间的流体进行导热。假设颗粒的周围包裹了一层厚度δ = d的气膜，颗粒-颗粒之间的导热通过两个颗粒共同占有的那部分气膜进行传热，气膜的厚度假设为颗粒粒径的1/5。

 

图1 颗粒-流体-颗粒导热示意图

图1给出了具有不同粒径的颗粒-流体-颗粒导热示意图，存在如下关系式：

• 方程13
  

基于如上关系式，积分上限可以化简为：

• 方程15
  

因此，颗粒-流体-颗粒的导热传热量可以表示为：

• 方程16

其中，积分下限写作：

• 方程17

04

辐射传热

辐射传热发生于颗粒表面和周围的介质（颗粒或流体）之间。具体而言，稠密气固系统中的辐射传热发生于颗粒和壁面/颗粒群，以及颗粒和流体之间。在CFD-DEM计算时，通常假设壁面为绝热条件（亦即，无热能从系统损失或进入）。因此，不考虑颗粒和壁面之间的辐射传热。此外，假设气体为理想透明（亦即，气体对辐射传热无阻碍或贡献）。因此，也不考虑气体和颗粒之间的辐射传热。当系统内的颗粒温度低于600 °C时，颗粒-颗粒间的辐射传热可以忽略不计。然而，当温度升高时，颗粒-环境辐射会逐渐增强。对于流化床内生物质气化燃烧的CFD-DEM数值模拟，显然需要考虑辐射传热带来的影响。

 

图2 辐射传热示意图

图2给出了辐射模型示意图。假设颗粒只和其周围的一个子域（红色虚线）内的颗粒发生辐射传热，则可以得到下列公式：

• 方程18

其中，ε_p,i为颗粒i的发射率，σ为斯特潘-玻尔兹曼常数，大小为5.67032×10^-8 W/(m²·K⁴)。T_p,i和T_env分别为颗粒温度和环境温度。许多的研究者在其CFD-DEM工作中采用上述公式计算辐射传热。最具争议性的问题集中于环境温度T_env的计算方式。一种方式是取系统内所有颗粒温度的平均值作为环境温度。这种处理方式表明颗粒i和系统内其它所有颗粒均发生了辐射传热。但是这种处理方式具有一个缺陷：同一时刻内处于系统不同位置的每个颗粒都有几乎相同的环境温度。假设在一个冷颗粒群和一个热颗粒群中有两个温度相同的颗粒，因为其环境温度相同，那么其辐射传热相同，这和现实情况是不符的。因此，如图2所示，将颗粒i周围1.5倍粒径范围内的局部颗粒温度平均值和局部气体温度（T_f,Ω）的综合值作为环境温度，引入流体空隙率进行关联：

• 方程19

其中，N_p,Ω是质心处于此局部范围内的总颗粒数，T_p,j为此局部范围内颗粒j的温度。

对于流体而言，采用温度（T_f）表示的能量守恒方程如下：

• 方程20
  

其中，ρ_f 为流体密度，ΔH_rf 为化学反应热，γ_Rf 为气相辐射传热系数，T_Rf 为气相辐射温度。方程左侧第一项表示了在控制体内热能的变化，第二项表示通过控制体表面传递的热能量。右侧第一项表明了对流传热量，通过傅里叶定律进行求解；第二项表示气固对流传热量；第三项表示了化学反应产生的热量；第四项表示气相的辐射传热量。

4

模型验证

Patil等^[3]做了一系列关于鼓泡床内气固传热的实验。因此，主要采用其实验数据来验证传热模型在复杂密相系统内的适用性。颗粒的运动学参数和热物性参数（例如，恢复系数、粘度和热导率等）参照实验设定。  

 

图3 准二维鼓泡床几何模型示意图

表1 准二维鼓泡床验证算例工况设置

物理参数

值

几何尺寸, x × y × z

8.0cm × 1.5 cm × 25 cm

网格数, Nx ×Ny × Nz

32× 6 × 100

物料质量, m

75g, 125 g

粒径, d_p

0.8mm, 1.0 mm, 1.2 mm

流体速度, U_bg

1.20m/s, 1.54 m/s, 1.71m/s

流体时间步长, Δt_g

1.0×10^-4s

固相时间步长, Δt_s

1.0×10^-6s

图3给出了Patil等[3]实验中所采用的准二维鼓泡床几何模型示意图。其长、厚和高分别为8 cm、1.5 cm和25 cm。对于拥有75 g物料颗粒的工况，颗粒数为57296个。温度为363.15 K的颗粒初始堆积在床层下部，随后引入温度为293.15 K的气体流化颗粒。对于速度边界条件而言，入口设置为均匀速度，壁面为无滑移；对于压力边界条件而言，出口设置为环境大气压。具体的研究工况设置列于表1中。根据Patil等[3]的实验，对于速度的影响，研究1.20 m/s、1.54 m/s和1.71 m/s三组工况；对于粒径的影响，研究0.8 mm、1.0 mm和1.2 mm三组工况。计算区域按照32 × 6 ×100划分网格，网格尺寸约为2.5 mm，大约是颗粒粒径的三倍。这样的网格和粒径比满足传统CFD-DEM计算的要求。流体相和固相时间步长分别设置为1.0 ×10-4 s和1.0 ×10-6 s。每个算例计算10 s，最后5s做统计平均。

 

图4 鼓泡床内瞬态颗粒运动和温度分布云图，Ubg = 1.20 m/s，m = 125 g：（a, c）实验结果；（b, d）本文模拟结果

图4给出了鼓泡床内瞬态颗粒运动和温度分布的实验与模拟对比图。从气泡大小以及床层膨胀而言，图4（b）的模拟结果和图4（a）的实验结果对比很好。床层底部颗粒很快的被流化气体冷却而床层顶部的颗粒依旧保持较高温度。此外，可以观察到鼓泡床下部中心区域存在一个狭窄的低温区域，此区域随着高度的增加逐渐消失（图4（c, d））。模型很好地捕捉到了这个颗粒温度变化的特征，验证了模型的合理性。

如图5所示，鼓泡床内的一个大气泡被拉伸撕裂成为两个紧靠的气泡，低温区域依旧存在于床体下部的中心位置。由于壁面的阻碍和限制，颗粒从边角处运动至中心和低温气体接触从而传递热量。低温气体穿过高温颗粒层时通过对流传热机制被显著加热。因此，气流在两个气泡中有较高温度，这和实验观测是一致的（图5（c, d））。

在定性对比模拟和实验所得的瞬态和时均云图之后，定量对比了不同工况下实验和模拟中鼓泡床内平均颗粒温度的演变（图6）。如图6（a）所示，增大流体速度会增大每个颗粒的对流传热速率，因此床体内的颗粒温度降低较快。此外，较大的流体速度使得颗粒运动更为剧烈，颗粒碰撞频率增加，从而颗粒-颗粒导热以及颗粒-流体-颗粒导热也相应增加，带走系统内更多的热量，颗粒较快降温。随着物料颗粒质量从75 g增加到125 g，相同速度下的颗粒降温明显减缓。这是因为对于拥有更多颗粒的系统，其固定的气体流量导致颗粒需要更长的时间来降温。

(a) 

(b)

5

模型应用

以实验室尺度的循环流化床（高度为1200 mm）为例，提升管下部的截面积为53 × 53 mm²，上部的截面积为58 × 58 mm²。旋风分离器的高度为230 mm，料腿的高度为955 mm。网格尺寸约为4.5 mm，是颗粒粒径的三倍，满足传统CFD-DEM方法关于网格和粒径比值的要求。初始时，在系统中产生200,000个颗粒，其中料腿和提升管中各放置一半数量，颗粒温度为298 K。提升管的中上部分为加热区域，设置为固定壁温373 K。将温度为298 K的冷气流从提升管底部引入循环流化床内。具体的气固物性参数及运行工况列于表2中。因算例中温度较低，因此未考虑辐射传热影响。

表2 流体和颗粒物性参数设置

流体属性

值

颗粒属性

值

密度[kg/m³]

1.225

密度[kg/m³]

1500

粘度[kg/(m·s)]

1.8×10^-5

粒径[mm]

1.6

温度[K]

298

温度[K]

298

比热容[J/(kg·K)]

1007

比热容[J/(kg·K)]

800

热导率[W/(m·K)]

0.0256

热导率[W/(m·K)]

1.0

流化风速[m/s]

6.0

杨氏模量[Pa]

5.0×10⁷

返料风速[m/s]

1.0

泊松比[-]

0.33

恢复系数[-]

0.90

图7给出了流化风速为6.0 m/s时循环流化床内颗粒混合以及温度和固含率演化云图。

   

   

   

图7 循环流化床内气固传热过程

低温的气体从提升管底部和料腿底部通入床内。由于强烈的湍流和气固动量交互，提升管内的气固运动混乱无章，在壁面附近形成条带状的颗粒团。此外，床体下部的颗粒温度较低。一般而言，循环流化床内存在着导热传热和对流传热。前者是通过颗粒-颗粒之间的相互碰撞接触传递热量，而后者主要是通过颗粒和流体之间的接触传递热量。提升管内颗粒-壁面的导热传热较小而对流传热（流体-壁面，颗粒-流体-壁面）则占据主导地位。

参考文献

[1] Li J, Mason D J. A computational investigation of transient heattransfer in pneumatic transport of granular particles [J]. Powder Technology.2000, 112: 273-282.

[2] Batchelor G K, O'Brien R W. Thermal or electrical conduction through agranular material [J]. Proceedings of the Royal Society of London. A.Mathematical and Physical Sciences. 1977, 355: 313-333.

[3] Patil A V, Peters E A J F, Sutkar V S, Deen N G, Kuipers J A M. A studyof heat transfer in fluidized beds using an integrated DIA/PIV/IR technique[J]. Chemical Engineering Journal. 2015, 259: 90-106.