流态化是利用流动流体的作用将固体颗粒群悬浮起来的一种物理现象[1]．流态化技术已经应用于工业生产[2]，最常见的流态化生产设备是流化床．目前常见的两类气固流化床为循环流化床(circulating fluidized bed，CFB)和鼓泡流化床(bubbling fluidized bed，BFB)．其中，鼓泡流化床具有较高的传热速率和良好的气固接触，在矿物筛选加工、固体干燥、药品生产等领域发挥着重要作用[3]．研究BFB中的流动特性对大型工业设备的设计优化和结构改进具有重要意义．对于气固流化床，一般认为颗粒受到曳力、升力、压力梯度力、重力和碰撞力[4]．Lu等[5]认为，曳力在上述力中占据主导地位．考虑到速度波动是影响BFB中曳力的重要因素[6-7]，一些学者尝试将速度波动添加到曳力模型中，最经典的代表速度波动强度的指标为颗粒拟温度．参考雷诺数的定义，对颗粒拟温度进行无量纲处理，使用颗粒拟温度雷诺数ReT作为曳力模型的添加项，ReT=ρgdpθ/μg，式中：ρg为气体密度；dp为颗粒直径；θ为颗粒拟温度；μg为气体动力黏度．相对应地，颗粒雷诺数为Rep=ρgdpug-us/μg，us为颗粒速度；ug为气体速度．Wylie等[8]发现颗粒的平均曳力随着ReT的增加而显著增大．Tang等[9]在TPKKV模型基础上，将PR-DNS计算的曳力减去TPKKV模型的曳力得到一个差值，认为这一差值与ReT和固含率成正比，与空隙率的二次方成反比，将ReT添加到TPKKV模型中，即为TPK模型．Huang等[10]发现，颗粒速度波动对曳力的影响非常显著，并根据PR-DNS数据提出了HWZL模型．Luo等[6]根据PR-DNS数据，提出了一种基于人工神经网络的曳力模型(LWJ模型)，发现添加颗粒拟温度和空隙率波动ς后曳力模型的计算精度大幅提高，其中颗粒拟温度的作用效果远大于空隙率波动．目前，有关ReT改进曳力模型的研究都是基于PR-DNS数据拟合新的曳力模型[11-13]．这种方式计算资源耗费大、适用性有限，不同流化床都须进行曳力模型修正，因此迫切须探索能满足需求的曳力模型改进方法．采用CFD-DEM计算结果改进曳力模型能兼顾计算精度和计算量，是一种具有前景的改进方案．然而，现有的CFD-DEM模拟中曳力模型的研究尚不充分，缺乏验证．为此，本研究将参考TPK和LWJ曳力模型的改进，根据先前得到的CFD-DEM数据[11]，在BVK曳力模型的基础上添加ReT改进曳力模型．1 正交试验正交试验法是一种常用的分析多因素多水平的方法，广泛应用于土木建筑[14]、自动控制[15]、数值模拟优化[16]等领域，具有变量分布均匀、数据可比性强的特点，仅需少量试验就能获得科学准确的结果，可代替完整的试验[17]．TPK模型可以把添加的颗粒拟温度项改写成ReT(2.98εg-2-2.98εg-1)，εg为孔隙率．LWJ模型可以把添加的颗粒拟温度项改写成ReT(2.296εg-2-3.387εg-1)．将颗粒拟温度项写成ReT与空隙率多项式相乘是可行的方法，为防止多项式线性回归中出现过拟合问题和模拟工作量的限制，空隙率的范围不能太大，在BVK模型的基础上进行改进，具体模型如下FBVKM=FBVK+ReT(Aεg-3+Bεg-2+Cεg-1+D+Eεg),式中：FBVKM为改进BVK模型的曳力；FBVK为BVK模型的曳力；A~E为系数．将空隙率的-3，-2，-1，0，1次项分别视为5个因素，分别给定5个水平，即可制作L16(54)的正交试验因素水平表．模型和模拟参数设置与文献[11]保持一致．须说明的是，对BVK模型的改进增加了颗粒所受的曳力．为了避免颗粒逸出模拟域，将y方向高度适当增加至0.18 m．正交试验的难点在于找到一个合适的范围，使模拟得到的结果是正常的．经过大量试算，最终确定正交试验因素水平如表1所示.10.13245/j.hust.240668.T001表1因素水平表水平ABCDE139.5-230.5490.5-454.0152.5240.0-230.0491.0-453.5153.0340.5-229.5491.5-453.0153.5441.0-229.0492.0-452.5154.0表2和表3为正交试验的方案和结果，评价实验方案的指标是时均相间能量传输速率(P)．与文献[11]中的计算类似，选取10.0~30.0 s模拟得到的相间能量传输速率的平均值．表3中：K¯ij为第j列第i个因素4个水平的平均值，K¯ij=∑i=14(aij/4）；Rj为第j列K¯ij的极差值，Rj=max{K¯ij}-min{K¯ij},Rj越大，该因素对试验结果影响越大．10.13245/j.hust.240668.T002表2正交试验方案工况ABCDEP/mW1111114.6822122224.0443133333.3544144442.8945212343.9186221433.8227234123.9078243213.8419313423.59210324313.85111331243.79812342133.80813414233.75314423143.53515432413.54716441323.78610.13245/j.hust.240668.T003表3正交试验结果参数ABCDEK¯1j3.7443.9874.0223.9833.980K¯2j3.8723.8133.8293.8593.832K¯3j3.7623.6513.5813.7273.684K¯4j3.6553.5823.6013.4643.537Rj0.2170.4050.4410.5190.44310-3由表2和3可得：最优参数组合为A1B4C4D4E4，即A=39.5，B=-229.0，C=492.0，D=-452.5，E=154.0．此时时均相间能量传输速率与PR-DNS结果(P=2.16 mW)相比，相对误差为33.98%，远小于原始BVK模型(P=3.87 mW)的相对误差(79.17%)，对应的改进曳力模型命名为BVKM-1，              FBVKM-1=FBVK+ReT(39.5εg-3-229.0εg-2+492.0εg-1-452.5+154.0εg)．极差大小排列为R4R5R3R2R1，所以5个因素中D影响最大，E次之，A最小．根据极差大小的排列，可以适当放大D的取值范围，缩小A的取值范围，在A=39.5，B=-229.0，C=492.0，D=-452.5，E=154.0附近寻找更优的参数组合．经过大量测试可得：当A=39.7，B=-230.1，C=490.9，D=-453.6，E=152.7时，时均相间能量传输速率P=1.474 mW，相对误差的绝对值更小(-31.76%)，对应的改进曳力模型命名为BVKM-2，                FBVKM-2=FBVK+ReT(39.7εg-3-230.1εg-2+490.9εg-1-453.6+152.7εg)．2 曳力模型改进前后结果对比将BVKM-1，BVKM-2，BVK模型和实验、PR-DNS数据及文献[11]进行对比．2.1　截面时均颗粒轴向速度y=15，25和35 mm高度截面上时均颗粒轴向速度(uy)分布如图1所示．可以看到：BVKM-1和BVKM-2的uy呈现床层中心高、两侧近壁面低的分布，uy最高值出现在x=22 mm附近．此外，与BVK模型相同，除了y=15 mm截面，两个改进模型的uy在床层中心位置始终低于实验和PR-DNS数据．10.13245/j.hust.240668.F001图1不同高度截面时均颗粒轴向速度uy分布尽管添加了颗粒拟温度项，但是BVKM-1，BVKM-2与BVK模型、PR-DNS一样，无法正确预测截面时均颗粒轴向速度在x=4，40 mm附近先减小后增加的趋势．然而，BVKM-1和BVKM-2显著提高了床层中心和近壁面的uy，尤其是BVKM-2，在3个高度截面上近壁面结果的改善都非常明显．在y=15 mm截面x=0.88 mm处，BVK，BVKM-1和BVKM-2与实验值的相对误差分别为63.95%，60.67%和26.39%，y=25 mm截面x=0.80 mm处分别为208.36%，119.53%和65.56%，y=35 mm截面x=1.01 mm处分别为93.73%，73.67%和14.03%，说明向曳力模型中添加颗粒拟温度项虽然计算量有一定增大，但可以有效改善uy的结果．BVKM-1与BVKM-2在床层中心的计算结果几乎没有区别，但BVKM-2在近壁面与实验值的相对误差明显小于BVKM-1，在截面时均颗粒轴向速度的预测上表现更好．2.2　压降压降为气流入口与出口面的压力差．当10~30 s时，BVK，BVKM-1和BVKM-模型的压降时均值分别为167.71± 43.26，174.45±32.42，171.17± 46.96 Pa，没有明显区别；BVKM-1和BVKM-2的时均值略高于BVK模型．3种曳力模型的压降频域信号如图2所示，其中：f为频率；I为压降能谱功率．压降能谱功率先上升后下降，最后在均值附近波动．3种曳力模型的主频都能在低频区找到，BVKM-1和BVKM-2的主频分别为3.348和3.249 Hz，都低于BVK模型的主频．10.13245/j.hust.240668.F002图23种曳力模型的压降频域信号总的来说，两种改进模型的压降在时域和频域上的结果与BVK模型区别不大，说明向曳力模型中添加颗粒拟温度项不会对压降产生显著影响．2.3　空隙率图3给出了z=0.05 mm截面上BVK，BVKM-1和BVKM-2的瞬时空隙率分布，其中每幅分图从左往右依次为BVK，BVKM-1和BVKM-2模型．BVK模型的云图与两种改进模型同宽但更短，这是因为改进曳力模型模拟时适当扩大了计算域，y方向长度从0.133 2 m增加到了0.18 m．10.13245/j.hust.240668.F003图3BVK，BVKM-1和BVKM-2的瞬时空隙率分布与压降特性不同，3种曳力模型瞬时结果区别较大，没有两种模型在所有时刻结果都非常相似，说明曳力模型中的颗粒拟温度项对空隙率产生明显影响．但在部分时刻，空隙率分布有一定的相似性．如t=1.0 s时，3种模型都预测了床层中心气泡的形成，其中BVK和BVKM-1的气泡位置和大小很接近．t=2.0 s时3种模型同时模拟出了气泡，但气泡大小区别很大．t=10.0 s时，BVK和BVKM-2模拟的气泡大小类似，但BVK的气泡贴着左墙，而BVKM-2的气泡贴着右墙．BVKM-1和BVKM-2的颗粒拟温度项系数上区别不大，但t=0.5 s时，两者预测的空隙率分布几乎完全相反，BVKM-1预测床层中心有颗粒团簇，而BVKM-2则预测有气泡生成．这说明了曳力模型的改进对参数非常敏感，系数上略有区别，模拟结果可能存在较大的差异．表4给出了3种曳力模型在不同高度截面上的时均空隙率．随着高度的增加，截面时均空隙率都在增大．在y=15 mm截面上，BVKM-1和BVKM-2预测的截面时均空隙率大于BVK．BVKM-1在所有截面上预测的时均空隙率都大于BVK，而BVKM-2在y=25和35 mm截面上时均空隙率均小于BVK．10.13245/j.hust.240668.T004表43种曳力模型不同高度截面的时均空隙率y/mmBVKBVKM-1BVKM-2150.5900.6550.624250.6650.6790.652350.7230.7240.708总的来说，3种模型预测的瞬时空隙率分布在部分时刻相似，整体模拟结果区别较大，但不同高度截面上的数值比较接近，整体的变化趋势都是随着高度的增加而增大．瞬时空隙率分布对曳力模型的选择非常敏感，而时均空隙率则对曳力模型不敏感．向曳力模型中添加颗粒拟温度项未改变上述结论．2.4　角速度图4给出了3种曳力模型不同高度截面y方向角速度ωy分布．以y=15 mm截面x=0 mm附近为例，随着x坐标的增大，BVK模型的ωy增大，BVKM-1的ωy减小，BVKM-2的ωy先减小后增大．在y=25和35 mm截面上也能观察到同样的情况．10.13245/j.hust.240668.F004图43种曳力模型不同高度截面y方向时均角速度表5给出了3种曳力模型截面时均角速度．BVK和BVKM-1都是先减小后增加，BVK和BVKM-2在y=35 mm截面上时均角速度都很大．但添加颗粒拟温度项后对预测的时均角速度杂乱无章，且与原始模型区别较大．这主要与改进后的曳力模型与原始模型计算结果存在差异，且随高度增加差异会被放大，以及当前对于时均角速度的研究不全面有关．10.13245/j.hust.240668.T005表53种曳力模型截面时均角速度y/mmBVKBVKM-1BVKM-215-0.055 840.272 190.662 0525-0.271 04-0.691 061.210 31353.372 310.142 312.748 02rad∙s-12.5　颗粒拟温度在向BVK曳力模型中添加了颗粒拟温度项后，分析两种改进模型的颗粒拟温度结果就变得非常重要．3种曳力模型在y=15，25和35 mm高度截面上的时均y方向颗粒拟温度θyy如图5所示．所有曳力模型的θyy最大值出现在x=6，38 mm附近，最小值出现在x=22 mm附近且较为接近，θyy总体呈M型分布．然而，添加颗粒拟温度项使得两种改进模型与BVK模型的θyy有较大区别．在y=15 mm截面x=0处，BVKM-1的θyy=0.030 87 m2/s2，而BVK和BVKM-2的θyy分别为0.040 83和0.040 66 m2/s2．10.13245/j.hust.240668.F005图53种曳力模型不同高度截面y方向颗粒拟温度此外，BVKM-2在x=6和38 mm附近的峰最高．在y=25 mm截面上，BVKM-1和BVKM-2的θyy都小于BVK，BVKM-1下降幅度大于BVKM-2，尤其是近壁面处．当x=0 mm时，BVK的θyy=0.057 45 m2/s2，而BVKM-1的θyy=0.034 59 m2/s2，为BVK的60.21%，BVKM-2也略低于BVK，θyy=0.051 51 m2/s2，为BVK的89.66%．在y=35 mm截面上，x=0和44 mm处，BVK模型的值最小，BVKM-2最大．具体而言，当x=0 mm时BVKM-2的θyy是BVK的1.423 2倍，而当x=44 mm时BVKM-2的θyy是BVK的1.209 4倍．另外，与y=15 mm截面的情况相似，BVKM-2在x=6和38 mm附近的峰是最高的．纵观3个截面，近壁面处θyy的变化较大，然而在床层中心，始终为BVKM-1BVKM-2BVK．尽管如此，3种模型在床层中心的θyy区别很小，可以认为不同模型床层中心的颗粒速度普遍较大，速度波动较小，导致床层中心的θyy较小，对曳力模型的选择不敏感．3种曳力模型时均颗粒拟温度(θ1)和全局颗粒拟温度(θ2)如表6所示．可以看到：对比每种模型在不同高度上的数据，3种曳力模型的θ1都随着截面高度的增加先增大后减小．值得注意的是，BVK在y=15和35 mm截面上的θ1几乎一样，而在y=25 mm截面上θ1=0.077 77 m2/s2，远大于y=15和35 mm截面．与BVK模型不同，BVKM-1和BVKM-2在3个截面上的θ1变化不大．BVKM-1在3个截面上的θ1都是最小的．BVKM-2在y=15和35 mm截面上的θ1都是最大的，BVK模型则在y=25 mm截面上最大．10.13245/j.hust.240668.T006表63种曳力模型时均颗粒和全局颗粒拟温度参数BVKBVKM-1BVKM-2θ1y=15 mm0.062 820.051 130.065 91y=25 mm0.077 770.056 950.074 29y=35 mm0.063 120.055 850.070 12θ20.002 090.002 360.002 32m2•s-23种模型的θ2分布则与θ1的分布大相径庭．BVKM-1在3个截面上的θ1都是最小的，但θ2最大，为0.002 36 m2/s2，比BVK(0.002 309 m2/s2)高出12.92%．两种改进模型的θ2都大于BVK模型的结果，这说明向曳力模型中添加颗粒拟温度项可能会在部分区域引起计算结果中的颗粒拟温度减小，但全局颗粒拟温度增加．此外，3个模型的θ2数值都为1×10-3数量级，而θ1为1×10-2数量级．实际上，讨论的y=15，25和35 mm截面都处于密相区，属于速度波动较大的区域，随着高度的进一步增加，颗粒数量将逐渐减少，速度波动也随之减少，这与θ1比θ2大1个数量级的情况相符合．3 结论a．相比BVK模型，两种改进曳力模型预测的床层中心和近壁面处截面时均颗粒轴向速度显著增加，大幅减小了与实验和PR-DNS数据的误差．在y=35 mm截面x=1.01 mm处，BVKM-2模型与实验值的相对误差仅为14.03%，改进效果非常显著．另外，改进曳力模型压降的时域和频域结果与BVK模型区别很小，空隙率时均分布也是如此．但添加颗粒拟温度项后瞬时空隙率、时均y方向颗粒拟温度及时均全局颗粒拟温度的预测与BVK模型区别较大．b．BVKM-2模型的表现比BVKM-1更好，尤其是在截面时均颗粒轴向速度上．整体上BVKM-2模型的结果比BVKM-1模型更接近原始BVK模型．向曳力模型中添加颗粒拟温度项时，瞬时结果及颗粒拟温度和角速度的时均结果会对多项式的参数非常敏感．