基于CCDM的流态化制冰的多相流动与传热特性

梁坤峰，阮春蕾，王林，袁竹林

(1 河南科技大学制冷与空调技术研究所，河南 洛阳 471003；2 东南大学能源与环境学院，江苏 南京 210096)

引 言

基于流态化原理实现直接接触式制取流体冰的方法是一种新型的动态制冰方式[1-2]，制冰过程中循环流动载冷液体与液滴的直接接触换热，具有较高的热效率。为了研究流态化制冰过程中的颗粒运 动、相变与碰撞聚并，以及循环载冷液体的流场与温度场分布特征，掌握制冰过程的控制方法，需要探讨流化床内液滴、冰颗粒和循环载冷液体所构成的颗粒-流体多相系统的流动与传热特性。但由于多相系统相间相互作用的复杂性，存在观测困难，单纯采用实验手段进行其内部规律性研究，难以完整阐释。随着计算机的高速发展，数值模拟方法以其较强的再现能力、低耗费、零危险性等优势，在对多相流动与换热体系的研究中发展很快。

近年来，一种改进的基于欧拉-拉格朗日方法的CCDM（combined continuum and discrete model）模型开始应用于颗粒-流体多相系统的研究[3-6]，该模型结合了用于计算颗粒的离散相模型（discrete phase model，DPM）和模拟连续相的计算流体动力学模型，既采用欧拉法对连续相的宏观场特征进行求解，又基于欧拉-拉格朗日法跟踪求解场中的每一个真实颗粒，反映颗粒在流动过程中的微观细节，很好地描述颗粒-流体多相系中颗粒运动、传热与传质的微观特性。自从Tsuji 等[7]首先提出CCDM 模型以来，已在颗粒-流体多相系统模拟研究中得到关注；Swasdisevi等[8]对二维喷动床内的气固传热特性进行了研究，考虑了气固相间的对流传热，并在喷动区和下降管中采用不同的传热系数；Li 等[9]考虑了颗粒-流体间的对流传热和颗粒-颗粒间的热传导，基于CCDM模型研究了水平管中的气固流动与传热问题；Zhou等[10]采用CCDM 模型研究了气固流化床内的煤粉燃烧问题，对气相湍流采用了大涡模拟方法，相间对流传热和颗粒-颗粒的接触导热也得以考虑；而采用CCDM 模型在液-固两相流体系的流动与传热方面的研究较为缺乏，Malone 等[11]采用准三维的CCDM模型研究了液-固流化床内的流动与换热问题，比较了两个不同的相间传热系数模型的应用条件。

本文建立了描述颗粒-流体多相系统的流动、传热与相变的CCDM 模型，基于欧拉法求解连续相的宏观场特征，由欧拉-拉格朗日法跟踪多相场中的离散颗粒相的运动、换热与相变特性；采用两套不同尺度的计算网格准确描述宏观场特征与微观颗粒行为，实现双向耦合作用；获得了流态化床内颗粒-流体多相系统的流动与传热特性，所获结论为流态化制冰过程的运行设计提供了理论指导。

1 CCDM 模型建立

1.1 基本方程

1.1.1 连续相守恒方程

（1）质量守恒方程

（2）动量守恒方程 在流化床内，由于离散相液滴与连续相载冷液体间局部的流体温度梯度将引起局部的流体密度差，考虑到床内多相流动的温度梯度不太大，采用Boussinesq 假设，在动量守恒方程中的重力项中将流体的密度处理为温度的线性函数[11]，控制方程为

（3）能量守恒方程 流场局部平均的流体能量守恒方程为

1.1.2 离散相的运动与热平衡方程 床内离散颗粒相的运动与传热由DEM 模型求解。根据牛顿运动定律，颗粒相的运动方程为

忽略热辐射传递量、颗粒-颗粒碰撞接触的导热量，离散颗粒的热平衡方程为

在较小的时间尺度Δt内，假设传递给颗粒i的热量保持不变，则有

1.2 颗粒-流体相间全耦合机制

1.2.1 力的耦合 与气-固流动体系相比，液-固两相流的相间作用力要更复杂，需要考虑的力有：稳态流体曳力、虚拟质量力、Basset 力、剪切提升力、旋转提升力、压力梯度力、浮力。基于文中研究的颗粒-多相系统特征，可忽略Basset 升力和剪切提升力。相间作用力表达式如下

表1 相间作用力系数Table 1 Equation of force between liquid and solid phases

式中系数与相关指数说明见表1。

1.2.2 热量的耦合 对流传热通常是连续流体相中主要的传热模式，传热公式常常是建立在实验数据上的经验或半经验公式。对于颗粒-流体多相流动体系，Haid[16]根据实验数据拟合出了三系数模型，预测的相对误差为34.4%；Bremford 等[17]实验研究了床内黏度相对较大的连续相流体中的相间传热系数，指出固相颗粒的存在对流体黏度和相间传热的影响是不可忽略的；而Agarwal[18]通过引入床层空隙率函数，将描述单个球形颗粒在均匀流场中强制对流的Ranz-Marshall 关联式应用到颗粒-流体多相系统中，实现了流化床体系中颗粒对流传热的Nusselt 数和在相同表观速度流体速度下单个颗粒的Nusselt 数的关联，并以此计算床内液滴的对流传热系数

离散颗粒和周围流体间的对流换热由下式确定

1.2.3 颗粒的碰撞 流化床内液滴与载冷液体之间的流动换热，促使其表层快速冻结，冻结一定程度后的液滴发生碰撞时，难以发生颗粒的聚并，类似弹性刚体间的碰撞，服从胡克定律。由此，文中通过判定液滴冻结程度的深或浅，采用不同的碰撞模型计算，液滴冻结程度浅时，碰撞采用经典的O’Rourke 定律，而冻结程度深时，液滴间的碰撞则采用软球模型。

当连续相为黏度较大的液体介质时，由于液体介质在颗粒表面所产生的流动润滑作用，颗粒之间的碰撞与在气固流动体系存在着较大的差异。Joseph 等[19]以及Malone 等[11]通过对颗粒的恢复系数e进行修正，来考虑黏性液体介质对颗粒碰撞产生的影响

式中，Stc为颗粒碰撞产生反弹的关键Stokes损伤值，计算中取Stc=10[20]；Sti是与黏性流场相关的量

1.3 求解方法

1.3.1 计算条件 数值模拟的对象为实验研究装置中的测试床段，床体的内径为40 mm，考虑到计算时间，根据所考察内容，将床体高分为2 m 与8 m（实验装置中，竖直床体段高2 m，而从液滴雾化形成的床体底部到冰晶与载冷液体混合物出口处的全距离为8m）两种情况。

整个求解区域划分非结构六面体网格，网格总数分别为8000（2 m）和32000（8 m）。床体的入口为速度入口边界条件，出口为压力出口边界条件，而壁面采用无滑移、绝热边界条件处理。颗粒的喷入位置为床体底部0～0.03 m 区间，采用随机方法确定喷入颗粒的初始位置，应用Rosin-Rammler 分布函数确定喷入颗粒的初始粒径分布。计算中所设计到的相关参数及条件见表2。

1.3.2 分割式搜索法 目前，对于离散颗粒受力以及颗粒对流场反作用的求解，较为普遍的是采用一种简单近似、结构网格统计平均的方法，但对于具有复杂外形几何体的非结构网格而言，难免会产生误差，甚至极大地偏离真实情形。文中采用分割式搜索法确定颗粒所在的当地非结构网格，颗粒-流体间相互作用均在该非结构网格中进行计算，实现一对一的双向耦合作用与反馈。以二维网格为例说明如下：将不规则外形几何体，用较大的规则几何体覆盖，如图1(a)所示，不规则椭圆区域为求解区域，规则四边形为覆盖区域。将椭圆区域内的非结构网格采用一维编号技术进行编号[21]，并计算其所在结构网格的编号。在求解区域内，每个结构网格中包含至少1 个不规则非结构网格，然后对其中的非结构网格依据扫描顺序重新进行编号，如图1(b)所示，并建立数组进行存储。在每个求解时间步，根据颗粒物理位置，确定其所在的结构网格。搜索仅在该结构网格及其邻近结构网格中进行，若满足lc=min(l1,l2,…,ln)，即为颗粒所在的当地非结构网格，其中l为非结构网格质心到颗粒的距离。采用该方法，能准确地计算颗粒-流体间的耦合作用与反馈，且搜索效率显著提高。

表2 模型的计算条件和参数Table 2 Calculation conditions and calculated parameters

图1 非结构网格颗粒搜索示意Fig.1 Non-structural grid used for searching particles

1.3.3 计算流程 采用先分离后耦合的算法，即颗粒相与连续相场分开计算，而耦合计算固相颗粒场 与连续相流场之间的相互作用。计算步骤为：首先，获得一个稳定的初始连续相流场；接着将离散颗粒喷入连续相流场中，根据颗粒所在的位置计算整个流场区域中每个控制体积的空隙率（即连续相所占的体积分数）；然后计算离散颗粒所受的作用力和颗粒-流体之间的对流传热量，颗粒所受力包括流体作用力、颗粒间碰撞接触力、重力等；随后依据牛顿第二定律，在多个固相颗粒计算时间步内，更新颗粒的速度、位置、温度等属性，并在固相颗粒场计算结束时，再次计算流场中每个控制体积的空隙率，获得固相颗粒对流体的反作用力和总的换热量；最后以源项的形式将其加入连续相场的求解控制方程中，迭代求解连续相流场，进入下一循环。

2 实验装置

流态化制冰实验系统原理图如图2所示，整个实验系统由流化床、制冷机和测试装置3 部分组成。流化床的主体由有机玻璃管制成，流化床底部中心轴线处设置1 个喷射制冰液体的喷头，喷头以下 300 mm 处为另一种非相溶的循环流动的载冷液体入口。制冷机（5 hp 压缩机，1 hp=745.699 W）通过蒸发器给吸收制冰液体相变热的载冷液体提供 足够的冷量。测试装置包括温度[温度传感器采 用Pt100 型铠装热电阻，精度等级为A，允许偏 差±(0.15+0.002|t|)]、压力与流量（压力测量选用绝压变送器，静态精度为0.5%FS；载冷液体的流量采用LC-40 型齿轮流速计，精度等级为0.5，制冰液体流量由LZJ-10 型玻璃转子流量计测量）的测试与采集。

制冰过程在流化床内部进行，床的底部形成两个并列两路液体入口：循环载冷液体流动的油路与制冰液体喷射雾化的水路。循环载冷液体流路中，经制冷装置蒸发器降温后的载冷液体在离心泵的作用下，自蒸发器流经过滤器，最后流入流化床内。当载冷液体的温度、流速稳定后，打开增压泵将水强制喷入载冷液体中，水在非相溶的载冷液体中液-液雾化形成细小液滴，并与其一起向流化床上部同向流动。此时，可通过测试装置进行实验测量，如通过设置在床高0、0.5、1.0、1.5、2.0 m 处的温度传感器进行温度测量。

3 结果与分析

3.1 颗粒碰撞模型的判别准则

当液滴在雾化形成的初始阶段，由于液滴冻结程度浅，液滴大部分都处于液态，此时，液滴-液滴的碰撞极有可能引起液滴聚并现象的发生，可采用O’Rourke 经典的液滴碰撞模型进行求解；而当液滴冻结到一定程度后，液滴外壳坚硬，即使颗粒之间发生碰撞，但聚并的概率极小，此时可遵循刚性球体间碰撞机理进行计算。如何确定液滴冻结的程度，并以此作为两种碰撞模型使用的判别准则？文中提出冰晶冻结程度关键值ξ的定义：ξ=rc/R，式中，rc为等效半径，即在某一时刻颗粒表面的对流换热量等效于其由外向里冻结所需冷量时的颗粒半径；R为颗粒的半径。

图3为4 种冻结程度关键值下颗粒粒径的轴向分布情况。为方便比较不同关键值对颗粒粒径的影响，不同工况下冰晶平均粒径起始值按该床高处的实验值，统一起始点。由图可知，在相同床高处，冻结程度关键值越大，颗粒的粒径则越小，且随着冻结程度关键值的增大，液滴粒径沿床高而增大的幅度明显减小；表明冻结程度关键值较大时，液滴表面的换热能力强，液滴相变进行得更快，此时液滴间碰撞发生聚并的概率就小，由此颗粒的粒径不会有较大幅度的增加。为了确定计算中冻结程度关键值的大小，计算了几个冻结程度关键值下颗粒平均粒径沿床高的变化曲线，通过与相同工况下的实验结果对比，发现采用冻结程度关键值ξ=10 时的床内颗粒平均粒径沿床高的变化趋势与试验结果比 较吻合，因此，在之后的计算分析中统一采用该颗粒冻结程度关键值来计算流化床内离散颗粒相与连续油介质之间的流动换热。

图3 冻结程度关键值对颗粒粒径沿轴向分布的影响Fig.3 Effect of particles freezing degree on axial size distribution (Tcont=-7 ,℃Gcont=0.4 L·s-1,Gdisp=30 ml·min-1)

图4 载冷液体轴向温度分布Fig.4 Axial temperature distribution of coolant

3.2 循环床内多相流动与传热特性

3.2.1 载冷液体温度分布 图4为流化床从起始到稳定的不同时刻载冷液体沿轴向的温度分布曲线。在雾化液滴进入流化床后的很短时间内（t=0.5 s），离散液滴还没来得及分布于整个流化床内，与连续相间只在床的底部换热，床内载冷液体的温度呈先增大后减小，继而保持初温不变；运行5 s 时间后，受离散液滴在床内扩散运动的影响，载冷液体的温度分布曲线峰值朝床高方向移动；运行50 s 时间后，床内液-液-固三相流动与换热体系基本达到稳定，载冷液体沿床高的温度增长趋势有所平缓。

此外，图4还给出了相同工况条件下的实验测试结果，对比发现模拟与实验结果的总体变化趋势一致，最大误差在13%左右，基本验证了CCDM 模型对于流态化制冰过程的模拟能力。分析局部出现误差的原因有两方面。一是实验数据的处理方面，实验的原始数据去掉了流化床跑冷的影响，而分析时采用的跑冷系数本身是根据实验数据拟合得到的[22]，在不同的阶段和工况下易造成误差的放大和减小；二是传热模型的应用方面，文中采用了Agarwal[18]提出的传热关联式，虽然该模型引入了床层空隙率函数将流化床体系中颗粒对流传热的Nu和在相同表观速度流体速度下Ranz-Marshall 关联式中单个颗粒的Nu关联起来，但应用于低温低黏度载冷液体与颗粒的表面传热背景，会对计算结果会产生一定影响。

图5 载冷液体径向温度分布Fig.5 Radial distribution of temperature of coolant

图6 离散颗粒相轴向温度分布Fig.6 Axial distribution of temperature of discrete particle

图5为连续相载冷液体在不同床高处沿径向的温度分布曲线。由图可知，在床的中心区域，连续 相沿床高方向的温度分布基本呈逐渐升高趋势；而在近壁区域，载冷液体的温度分布则比较紊乱，尤其是在床高1 m 处。原因在于床内离散颗粒相的紊乱运动，尤其是在床的上半段，由于颗粒向上运动速度的减慢，颗粒数密度增大，增加了颗粒间碰撞的概率，而此时大部分颗粒处于未冻结完全或冻结程度很浅的状态，颗粒之间的碰撞极易引起颗粒间聚并现象的发生，颗粒间的频繁碰撞与聚并将使其在床内的运动变得紊乱，颗粒的运动速度以及相间的能量耦合作用，都呈较大幅度的波动，这些因素将导致连续相温度场分布的不均匀。

3.2.2 颗粒相轴向温度分布 图6为离散颗粒相轴向温度分布曲线，反映了不同床层高度内颗粒平均温度的变化特征。由图可知，床内离散颗粒相温度分布的波动较为紊乱，但总的变化趋势仍与水的结晶过程曲线较为一致。颗粒相进入流化床后的很短时间内其平均温度便迅速下降，由初温5℃降至 0.5℃，这是由于颗粒入床时，颗粒速度与温度均处于最大值，相间较大的温度差以及较大的流速比，加剧了相间的传热，使颗粒温度大幅度降低。到达床高2 m 时，床内部分颗粒已达到相变温度。在床高2～7 m 区间，床内大部分颗粒均处于相变传热阶段，在床高7 m 以后，床内大部分颗粒基本上已完全冻结，此时，床层离散颗粒相的温度降至冰点温度以下，温度下降的幅度也比较平缓。

3.2.3 颗粒相轴向速度分布 图7为不同粒径时颗粒轴向速度分布曲线，选取了6 种粒径档的颗粒作为示踪颗粒，考察了其在床内流动速度的变化情况。此时制冰液体的喷射流率较小，为3.3 ml·min-1，由于喷射流率较小，床内颗粒数密度低，颗粒间发生碰撞聚并的概率较低。由图可知，在粒径档小于3.12 mm 时，颗粒在床内的运动速度呈先急剧减小，后逐渐增大，最终达到其在床内流动的终端速度；随着颗粒粒径的增大，这种变化趋势逐渐减弱，当颗粒粒径档大于3.62 mm 以后，颗粒的轴向速度则呈持续增大的变化趋势；但同时需注意不同粒径档颗粒的终端速度差异极小。出现这种情况的原因主要在于颗粒的流化过程中，喷射初始动量、自身重力、载冷液体的携带作用之间的相互作用所致；颗粒在床内雾化形成后，其运动速度在很短的时间内急剧减小，随后由于其速度的降低，其受到周围载流液体的流动携带作用，运动速度又逐渐增大，最后达到颗粒的终端速度。

图8为离散颗粒相轴向平均速度沿床高的变化曲线，此时喷射流率较大，为35 ml·min-1。由图可知，离散相喷射流率较大时，由于床层的固相颗粒数密度分布增大，床内颗粒间碰撞、聚并与结团等现象的发生概率增加，致使颗粒相在床内的运动变得紊乱，其轴向平均速度沿床高方向的分布相应 地也变得紊乱；尤其是在床高2 m 以下，由于颗粒运动速度的减慢，离散颗粒相轴向平均速度变化幅度较大。但从总体变化趋势看，离散颗粒相的轴向平均速度呈现先增大，随后在同一水平线上下波动的趋势。

图7 不同粒径时颗粒轴向速度Fig.7 Axial velocity of particles of different radius

图8 颗粒相轴向平均速度Fig.8 Axial average velocity distribution of particles phase

图9 颗粒相数密度轴向分布Fig.9 Axial number density distribution of particle phase

3.2.4 颗粒相浓度分布 图9为颗粒相数密度的轴向分布，图中曲线代表离散数据点的变化趋势。由图可知，床内固相颗粒的数密度沿床高而降低，在床高0.5～2.5 m 区间内，颗粒数密度变化幅度较大；达到床高6 m 以后，颗粒相数密度逐渐趋于平缓，基本保持不变。原因在于该高度区间内液滴未冻结完全，且运动速度减小，大量颗粒汇集在该区域，增加了颗粒间碰撞并发生聚并的概率，从而使该区域离散颗粒相数密度波动幅度较大。在床高2.5 m 以上，根据图8颗粒相轴向平均速度的变化趋势看，已基本趋于稳定，且已有部分颗粒的冻结程度较深，颗粒相的运动相对比较稳定，数密度也趋于稳定。

图10为颗粒相体积分数轴向分布，图中曲线代表离散数据点的变化趋势。与图9相对应，在床高0.5～2.5 m 区间内，由于颗粒的汇集，床内颗粒相体积分数达到最大，这一点与图8的颗粒相平均速度最小的表现一致；而在床高2.0 m 以上，由于离散颗粒流动的稳定，床层固相颗粒的体积份额也相应地达到稳定。

图10 颗粒相体积分数轴向分布Fig.10 Axial volume fraction distribution of particle phase

图11 不同粒径档颗粒的停留时间Fig.11 Residence time for particles of different radius

3.2.5 颗粒停留时间 颗粒在流化床内的停留时间是循环流化床制冰运行的重要参数。不同粒径档的颗粒在床内的运动特性存在差异，决定了其停留时间的差异；而在同一工况的相同操作参数下，必然会影响颗粒的冻结程度，并进而影响颗粒的碰撞聚并。因对不同粒径档的颗粒在床内停留时间的分析，有利于对操作参数及系统的优化设计。图11为不同粒径档的颗粒在流化床内停留时间的统计平均值，计算时间为100 s，其中，图中曲线代表离散数据点的变化趋势。由图可知，颗粒在床内的停留时间随 其粒径的变化虽较为紊乱，但从总的变化趋势上看，仍呈现出随粒径的增大，颗粒停留时间先增大后减小的变化趋势，并以粒径档3 mm 为分界点。粒径档小于3 mm 的颗粒，由于其体积较小，冻结过程所需冷量也少，在流化床下半段的流动过程中与载冷液体换热即可冻结到较深的程度，颗粒在随后的流动过程中，即使碰撞，发生聚并的概率很小，从而在一定程度上减小了小粒径颗粒流动紊乱性；在此粒径档范围内，颗粒的停留时间主要取决于颗粒的大小，大颗粒较小颗粒而言，不易被携带，故其停留时间也就较长。然而当颗粒的粒径档大于3 mm时，其冻结时间增加很多，使得液滴流过其初始动量较强的流化床下半段，其冻结程度仍较浅；越过该高度后，由于液滴流速的降低，床内局部液滴浓度增大，增大了颗粒间的碰撞概率，从而对于冻结程度相对较浅的大颗粒而言，发生聚并、粘连、分散等现象的概率较大，使其流动变得较为紊乱与复杂。

3.2.6 相间平均传热系数 图12为相间平均传热系数的轴向分布，通过统计床高0.1 m 高度内所有颗粒的数量和表面传热系数，并对这段高度内的所有颗粒的表面传热系数进行加和平均获得，体现了流化床内相间传热系数随时间由非平衡而逐渐达到平衡的过程。在流化床运行4.2 s 时，床内离散颗粒相流动未稳定，床内相间平均传热系数沿床高呈先增大后减小的变化趋势，在床高1 m 处最小；随后平均传热系数在床高1 m 以上的床层间，沿床高而增大的幅度逐渐减小，最终在t=7.2 s 时达到稳定，床内平均传热系数沿床高基本不变。分析认为，计算的初始时刻，由于床内固相颗粒的数密度较小，初始喷入颗粒的运动受到颗粒间碰撞聚并现象的影响极小，不存在速度突增或突减的情形；在床的上 半段，由于床内连续相载冷流体与离散颗粒相间的换热量也很少，其温度处于较低的水平，相间温差随着床高而增大，由此表现为初始时刻平均传热系数在床的上半段沿床高增大。然而随着时间的变化，床内不同高度处的离散相颗粒浓度分布增大且达到稳定，颗粒相与连续相流动稳定，使床内连续相载冷液体沿床高的温度分布达到稳定，而相间的平均传热系数也相应地达到稳定。

图12 床层相间平均传热系数随时间的变化Fig.12 Change of average heat transfer coefficient with time

4 结 论

（1）基于所建立的CCDM 模型求解了流化床制冰过程所涉及的颗粒-多相系统的流动和传热特性，该模型充分考虑了颗粒碰撞对流化床内液固相间的动量与能量交换的影响，并基于一种高效的非结构网格搜索法，快速地搜索每一计算时刻离散颗粒在流场中所在的非结构网格，实现了连续相与离散相间一对一的相间耦合作用。

（2）由于颗粒间碰撞、聚并、结团等现象的存在，水的喷射流率与颗粒的数密度和颗粒发生碰撞形成大颗粒的概率的变化趋势一致，而颗粒的粒径对于颗粒的运动速度和停留时间具有重要影响，从而导致相间传热系数沿床高发生变化；尤其是在流化床的下半段，颗粒的聚并引起了颗粒运动速度呈2～3 倍的增加，而相间传热系数增大了近百倍，不仅如此，此阶段的床内固相颗粒的数密度和体积分数都比较大。由此可见，在流化床制冰研究中难以忽略颗粒间的碰撞与聚并现象对于流化床的多相流动与换热特性的影响。

（3）在流化床的底部入口处，处于液态的颗粒被喷入床中时，颗粒与冷媒之间存在较大的温度差与流速比，强化了相间的传热，颗粒相平均温度在进入流化床后的很短时间内就迅速下降；在随后的相变换热阶段，颗粒相温度降低的幅度比较平缓。当流化床处于稳定制冰工作时，较大的相间平均传热系数出现在床高0.5 m 以下，超过这个高度，沿床高几乎不再变化。通过对比载冷液体轴向温度分布的模拟结果与实验结果，发现仅有13%左右的相对误差，对所建立的CCDM 模型在模拟流化床制冰过程的能力方面进行了验证，同时表明文中提出的冻结程度关键值，对于判定颗粒间碰撞模型和颗粒相表面传热量的等效计算简单有效。

符 号 说 明

Ac——循环床层截面积，m2

CA,i——虚拟质量力系数

CD,i——流体曳力系数

CS,i——旋转提升力系数

cf——连续相流体的比定压热容，J·kg-1·K-1

cp,l——水的比定压热容，J·kg-1·K-1

cp,s——冰的比定压热容，J·kg-1·K-1

dd——雾化液滴粒径，mm

F——颗粒的受力，N

Ffp——流体-颗粒间动量交换项，N

fA,i——虚拟质量力，N

fB,i——相间密度差引起的浮力，N

fD,i——流体曳力，N

fP,i——流体压力梯度力，N

fS,i——旋转提升力，N

Gdisp——离散相流率，ml·min-1

Hb——循环床高度，m

h——相间对流传热系数，J·m-2·K-1

Ip——颗粒转动惯量，J·kg-1·m-2

Mp——颗粒所受合力矩，N·m

mp——颗粒质量，kg

Nu——Nusselt 数

Pr——Prandtl 数

p——连续相压力，Pa

Qfp——流体与颗粒能量交换项，J

q——相间传热率，W

R——颗粒半径，mm

Rep,i——颗粒流动Reynolds 数

St——Stokes 数

T——温度，K

Tc,in——冷油介质入床温度，K

Tf——连续相流体温度，K

T*f——参考温度，K

Tp——离散颗粒温度，K

Δt——迭代时间步长，s

uf——连续相流体速度矢量

up——离散颗粒速度矢量

vc,in——载流介质入床流速

vp——颗粒运动速度，m·s-1

vp,in——颗粒相入口速度，m·s-1

β——流体热膨胀系数，m·K-1

μ——动力黏度，kg·m-1·s-1

γ——亚松弛系数

ε——网格空隙率

λf——连续相流体热导率，W·m-1·K-1

ξ——冻结程度关键值

ρ——密度，kg·m-3

ρ*f——参考温度下的流体密度，kg·m-3

τ——时间步长，s

ωp——颗粒转动的角速度，rad·s-1

Ω——量纲1 角速度

[1]Zheng Keqing (郑克晴),Zhang Xuejun (张学军),Tian Xinjian (田新建),et al.Heat transfer performance of a single air bubble in direct contact ice slurry generator [J].CIESC Journal(化工学报),2010,61 (S2):58-61

[2]Peng Zhengbiao,Yuan Zhulin,Liang Kunfeng,et al.Ice slurry formation in a cocurrent liquid-liquid flow [J].Chinese Journal of Chemical Engineering,2008,16 (4):552-557

[3]Li Hongchang,Li Yaoming,Gao Fang,et al.CFD-DEM simulation of material motion in air-and-screen cleaning device [J].Computers and Electronics in Agriculture,2012,88:111-119

[4]Kevin F Malone,Bao H Xu,Michael Fairweather.The combined-continuum-and-discrete-model (CCDM) for simulation of liquid-particle flows [J].Computer Aided Chemical Engineering,2006,21:227-232

[5]Chu K W,Wang B,Xu D L,et al.CFD-DEM simulation of the gas-solid flow in a cyclone separator [J].Chemical Engineering Science,2011,66 (5):834-847

[6]Liu Xiangjun (刘向军),Shi Lei (石磊),Xu Xuchang (徐旭常).Activities of dense particle-gas two-phase flow modeling in Eulerian-Lagrangian approach [J].Chinese Journal of Computational Mechanics(计算力学学报),2007,24 (2):166-172

[7]Tsuji Y,Tanaka T.Cluster patterns in circulating fluidized beds predicted by numerical simulation (discrete particle model versus two-fluid model) [J].Powder Technology,1998,95 (3):254-264

[8]Swasdisevi T,Tanthapanichakoon W,Charinpanitkul T,et al.Prediction of gas-particle dynamics and heat transfer in a two-dimensional spouted bed [J].Advanced Powder Technology,2005,16 (3):275-293

[9]Li J T,Mason D J,Mujumdar A S.A numerical study of heat transfer mechanisms in gas-solids flows through pipes using a coupled CFD and DEM model [J].Drying Technology,2003,21:1839-1866

[10]Zhou H,Flamant G,Gauthier D.DEM-LES of coal combustion in a bubbling fluidized bed (Ⅰ):Gas-particle turbulent flow structure [J].Chemical Engineering Science,2004,59 (20):4193-4203

[11]Malone K F,Xu B H.Particle-scale simulation of heat transfer in liquid-fluidized beds [J].Powder Technology,2008,184 (2):189-204

[12]Felice R D.The voidage function for fluid-particle interaction systems [J].International Journal of Multiphase Flow,1994,20 (1):153-159

[13]Oesterle B,Dinh T B.Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers [J].Experiments in Fluids,1998,25:16-22

[14]Tsuji Y,Morikawa Y,Mizuno O.Experimental-measurement of the Magnus force on a rotating sphere at low Reynolds numbers [J].Journal of Fluids Engineering Transactions of the ASME,1985,107:484-488

[15]Cai X L,Wallis G B.A more general cell model for added-mass in two-phase flow [J].Chemical Engineering Science,1994,49 (10):1631-1638

[16]Haid M.Correlations for the prediction of heat transfer to liquid-solid fluidized beds [J].Chemical Engineering and Processing,1997,36 (2):143-147

[17]Bremford D J,Müller-Steinhagen H,Duffy G G.Heat transfer to kraft black liquor in a liquid-solid fluidized bed [J].Heat Transfer Engineering,2000,21:57-69

[18]Agarwal P K.Transport phenomena in multi-particle systems (Ⅱ):Particle-fluid heat and mass transfer [J].Chemical Engineering Science,1988,43 (9):2501-2510

[19]Joseph G G,Zenit R,Hunt M L,et al.Particle-wall collisions in a viscous fluid [J].Journal of Fluid Mechanics,2001,433:329-346

[20]Gondret P,Lance M,Petit L.Bouncing motion of spherical particles in fluids [J].Physics of Fluids,2002,14 (2):643-652

[21]Cai Jie (蔡杰),Yuan Zhulin (袁竹林).Numerical simulation on solids fraction distribution in straw burned circulating fluidized beds [J].Proceedings of the CSEE(中国电机工程学报),2007,27 (20):71-75

[22]Peng Zhengbiao (彭正标).Experiment and modeling of ice slurry production in liquid-liquid calculating fluidization [D].Nanjing:Southeast University,2009