关于稠密气固两相流，你get到多少？
什么是稠密气固两相流？
气固两相流的分类通常以气固两相流中固相所占的体积分数（εp）作为分类标准[1]，如下表所示。
气固两相流
固相体积分数
耦合标准
耦合变量
稀相流
εp＜10-6
单向耦合
流体⇒颗粒
中等浓度流
10-6 < εp < 10-3
双向耦合
流体⇔颗粒
密相流
εp＞10-3
四向耦合
流体⇔颗粒⇔颗粒
当固相体积分数小于10-6时，气固两相流中颗粒之间距离较远，颗粒主要受到流体的作用力，而颗粒本身对于流体运动的影响可以忽略不计，因此可以采用单向耦合进行考察，称之为稀相流。
当气固两相流中的固相体积分数进一步增大，颗粒的群体效应增强，颗粒和流体之间的相互作用也进一步明显，此时需要采用双向耦合进行考察，称之为中等浓度流。具体而言，需要同时考虑流体对颗粒的作用力和颗粒对于流体的反馈作用力。
当固相体积分数进一步增大，颗粒之间由于距离较近而发生剧烈的碰撞，且颗粒之间的尾涡相互作用，会显著影响到周围流体的运动，这种以碰撞为主导的气固两相流需要采用四向耦合进行考察，称之为密相流。具体而言，在双向耦合的基础上，需要考虑颗粒-颗粒/壁面的相互碰撞以及它们之间的相互耦合。
稠密气固两相流广泛存在于工业实际中，如循环流化床锅炉、鼓泡流化床气化炉、化学链和炼铁高炉等。相对固定床而言，流化床具有气固接触效率高、化学反应速率强和污染物产量低等诸多优势。流化床主要有两种形式：鼓泡流化床（BubblingFluidized Bed, BFB）和循环流化床（CirculatingFluidized Bed, CFB）。
流化床是一个典型的稠密气固两相反应系统，包含了多尺度相间耦合和多物理过程耦合。就多尺度而言，有连续相湍流涡团的多尺度、离散相颗粒的多尺度、燃烧火焰的多尺度等；就多物理过程而言，存在着湍流流动、颗粒相互作用、传热传质、化学反应、凝并破裂等。需要有针对性地采用多尺度模拟策略。
如何模拟稠密气固两相流？
图1给出了Shaffer等[1]关于不同尺度下流化床内气固运动的示意图。
从空间尺度而言，微观尺度（Micro-scale, mm）上的气体和颗粒相互作用以及颗粒和颗粒的碰撞反映到介观尺度（Meso-scale, cm）上就会呈现出气固流动非均匀性，亦即会发生气泡的破碎、颗粒的聚团等现象，进而决定床内宏观尺度（Macro-scale, m）的气固流动特性，如气固相间接触效率、停留时间分布等。
从时间尺度而言，颗粒在流体中的弛豫时间以及颗粒之间的碰撞时间都是毫秒量级的，而颗粒在设备内的流化周期（如提升管内停留时间、整体循环时间等）又都是秒量级的。因此，稠密气固两相流这种在空间和时间上的复杂多尺度特性对数值模拟提出了较高的要求，需要相应地采用多尺度求解策略。
图1（b）直观给出了基于不同尺度的三种典型模拟方法：颗粒可解析直接数值模拟；欧拉-拉格朗日方法；欧拉-欧拉方法。
三类方法的具体数值模型分别为：
1 | 颗粒可解析直接数值模拟 (Particle-resolvedDirect Numerical Simulation, PR-DNS） |
2 | 计算流体力学-离散单元法（Computational Fluid Dynamics-Discrete ElementMethod, CFD-DEM） |
3 | 颗粒单元体积方法（Multiphase Particle-in-cell, MP-PIC） |
4 | 欧拉-欧拉框架下的双流体模型（Two-fluid Model, TFM） |
三类方法在计算尺度上依次增大，在计算时间花销上依次减小。以下进行具体分析。
1. PR-DNS
PR-DNS模型在欧拉框架下求解气相运动，在拉格朗日框架下求解固相运动。
从颗粒尺度而言，稠密气固两相流中的颗粒显著地受到其周围非稳态气流运动的影响，因此需要对气相运动进行精密的求解。
在实际的计算中，往往采用小于颗粒粒径数倍（< 1/10 dp）的网格来捕捉颗粒周边流场的信息。颗粒表面常假设为无滑移边界条件，通过对作用在其上的应力张量进行积分而获取相间作用力。
PR-DNS在计算流体运动时不需要引入任何模型进行封闭，因此求解精度高。众所周知，数值模拟的计算时间和所用网格数紧密相关。该方法所需的大量网格会造成极大的计算资源需求，使得此方法只能局限于较小尺度的模拟区域，且计算的颗粒数也受到相应限制，一般小于O(103)[3, 4]。
2. CFD-DEM
CFD-DEM模型是基于欧拉-拉格朗日框架下的一种计算方法。其气相运动在欧拉框架下进行求解，而固相颗粒运动在拉格朗日框架下进行求解。和PR-DNS不同，该方法要求网格直径是颗粒粒径的3~5倍[5-9]，网格量的减少造成了计算量的相应减小，其对气相运动的求解精度略低于PR-DNS。
CFD-DEM方法追踪每个实际颗粒的运动轨迹。通过引入空隙率，将每个颗粒的受力通过插值反馈到计算网格上求解气相运动的Navier-Stokes控制方程，从而获取每个网格内气相的速度和压强等重要信息。通过引入离散单元法中的硬球及软球碰撞模型，可以精确获取每个颗粒的位置、速度、旋转和受力等颗粒尺度的信息。
此外，因为每个颗粒单独进行计算，所以可以很容易添加诸如润滑力、磁场力和粘性力等受力模型，同时可以求解多粒径颗粒、非规则几何外形颗粒和反应颗粒等。但是此方法面临的问题是：
这3个问题都会显著降低CFD-DEM的计算效率。目前为止，采用此方法时，计算设备的尺寸受到其内颗粒数量的限制。但随着计算机性能的提升以及大规模并行算法的完善，CFD-DEM方法将会越来越多的被应用到稠密气固两相流的分析和预测之中。
3 MP-PIC
MP-PIC模型采用体积平均的Navier-Stokes方程描述气相的流动，采用计算颗粒思想离散固相。对每个计算颗粒进行追踪从而获得颗粒相在不同时刻的空间分布。此外，将颗粒属性插值到流体网格进行求解，从而使颗粒相同时具有离散相和连续相的性质。
该方法具有以下的特点：
（i）针对CFD-DEM在稠密气固两相流中颗粒碰撞处理方面需要消耗大量计算资源的问题，MP-PIC引入了固相应力模化颗粒碰撞，有效减小了计算量；
（ii）引入计算颗粒（Parcel）概念，每个数值粒子包含了数个具有相同运动形式的真实颗粒（Particle），并使用连续的固相应力模型阻止颗粒过度堆积而超出允许的物理极限；
（iii）将固相群体的特性映射到欧拉网格来计算连续颗粒应力场，而每个计算颗粒都具有独立的速度和位置，通过积分牛顿运动方程进行更新，从而使其同时具有了TFM和CFD-DEM两者的特点。基于这些特点，MP-PIC可以计算大规模颗粒系统，同时还能对固相在颗粒尺度进行描述。
4 TFM
双流体模型是欧拉-欧拉框架下最典型的计算方法。欧拉-欧拉方法采用连续介质假设，将固相也视作一种流体（“拟流体”），在欧拉框架下求解其运动。此方法的关键在于如何描述固相运动的粘度和压力，对此一般通过引入颗粒动理学（Kinetic Theory of Granular Flow, KTGF）理论进行模化。但是KTGF的本构方程十分复杂，目前并没有统一、通用且准确的模型解决该问题。
此外，稠密气固反应系统内的多粒径颗粒分布（Particle Size Distribution，PSD）很难在双流体模型中予以考虑。将固相颗粒视作为流体，客观导致了该方法对于模拟和预测稀相区气固流动存在不足。该方法只需要一套计算网格，每个网格大约是几倍或者十几倍颗粒粒径。相比于PR-DNS，TFM可以大规模减少计算量，近年来已被广泛应用于预测较大尺度工业设备内的气固流动。但是欧拉-欧拉方法存在较为明显的缺陷：对固相运动的求解精度仅限于计算网格尺度，因此无法获得丰富的颗粒尺度信息，从而难以用来研究流态化设备内颗粒的输运机制。
以上列举了三类（四种）最常见的基于不同尺度求解稠密气固两相流的数值模拟方法。这三类方法并不是孤立存在的，它们之间具有紧密的联系。一般而言，求解的尺度越小，则其模型精度越高，所能模拟的计算区域越小。亦即，小尺度对应高精度，而大尺度对应低精度。高精度的模型（PR-DNS）引入的模型更少，可以获取更精确的颗粒尺度信息（如气固相互作用力），从而为低精度模型（CFD-DEM，MP-PIC，TFM）提供封闭模型。
除了以上四种常见数值模拟方法外，针对稠密气固两相流动还存在着离散气泡模型（Discrete Bubble Model，DBM）、直接蒙特卡洛方法（Direct Simulation Monto-Carlo，DSMC）、光滑粒子流体动力学方法（Smoothed Particle Hydrodynamics，SPH）等其它方法。
稠密气固两相流模拟哪家强？
浙江大学能源工程学院计算组是国内较早开展稠密气固两相反应流数值模拟工作的课题组：王则力博士[11]和谭骏华博士[12]主要开展了稠密气固两相流的PR-DRS数值模拟工作，提出了多重直接力算法以及动态内缩算法；桂南博士[13]、张科博士[14]、房明明博士[15]、杨世亮博士[16]、王帅博士、胡陈枢博士先后开展了大量关于稠密气固两相反应流的CFD-DEM数值模拟工作。
课题组相关研究在《Journal of Computational Physics》、《AIChEJournal》、《Chemical Engineering Science》等SCI期刊上共发表50余篇论文，以下简要介绍关于PR-DNS和CFD-DEM方面的工作。
王则力博士
（1）在颗粒全尺度计算的内嵌边界方法基础上提出了内嵌边界多重直接力算法(Multi-direct Forcing Scheme)，并结合不同的计算方式提出了新颖的在笛卡尔均匀坐标系下的基于压力-速度方程的内嵌边界多重直接力算法和基于涡量-速度方程的内嵌边界多重直接力算法。该算法已经成功用于二维空间和三维空间中多相流的全尺度计算研究，并被国内外相关研究广泛采用。
（2）提出了内嵌边界多重直接热源算法(Multi-direct Heat Source Scheme)。该算法结合多重直接力算法已经成功的用于研究复杂边界与流体之间的传热研究以及多相流中运动颗粒与其周围流体之间传热的研究。
图2 全三维颗粒沉降（颗粒数：6000；网格数：1.2亿）
详细内容请参考王则力博士论文《多相流全尺度直接数值模拟方法及应用》。
谭骏华博士
（1）发展了适于大规模稠密气固两相流直接数值模拟的内嵌边界方法。以往针对颗粒计算水力直径的修正都是静态的且是基于系统平均、不区分颗粒状态差异的。本内嵌边界方法基于拉格朗日力点内缩策略和颗粒雷诺数计算方法，可对流场内每一个颗粒的计算水力直径实现实时的修正，这提高了颗粒流体两相流的模拟准确度，同时还具有良好的计算网格无关性，为稠密气固两相流的大规模直接数值模拟建立了基础。
（2）建立了气固流化床内颗粒曳力的直接数值模拟数据库，可为宏观模型发展提供基础数据。分析数据库发现传统曳力模型能从整体上捕捉流化床内颗粒曳力的变化规律，但平均值偏低。基于人工神经网络模型得到了四个物理参量对无量纲曳力的影响贡献，发现颗粒拟温度和空隙率脉动是影响流化床内颗粒曳力的两个因素，但是没有颗粒雷诺数和空隙率的影响大。
（3）开展了目前国际上最大规模的实验室尺度鼓泡流化床PR-DNS模拟，发现颗粒的存在导致密相区涡结构破碎更加剧烈，从而增强湍流耗散率。
详细内容请参考谭骏华博士论文《气固流化床内流动特性和颗粒曳力的直接数值模拟研究》。
桂南博士
主要进行了以下研究：
（1）揭示了惯性颗粒碰撞率的统计特性和惯性颗粒的基本特征，澄清并加深了对颗粒碰撞率的认识。同时，也揭示了外力场的效应和小颗粒的碰撞凝并的机理；
（2）设计了统一的波形鼓形态与运动的控制方程，并研究了波形鼓内颗粒群集的能量特性和角动量特性，分析了波形鼓的边界波数、振幅、转速对波形鼓内颗粒运动的影响。同时，对转鼓内颗粒混合界面的分形特征和混合过程的信息熵变化进行了研究；
（3）研究并揭示了流化床内气-固、固-固相互作用的若干规律，以及脉动条件下对床内相间相互作用的调制效应。同时，研究了颗粒-埋管碰撞和埋管磨损的基本特性以及在脉动和扰动气流对颗粒-埋管碰撞特性的影响。
详细内容请参考桂南博士论文《复杂两相流动中颗粒碰撞的DEM-LES/DNS耦合模拟研究》（2013年全国百篇优秀博士论文）。
张科博士
采用基于欧拉-拉格朗日框架的CFD-DEM方法，先后研究了拟三维鼓泡床、三维圆柱-锥体形喷动床、三维带导流管的圆柱-锥体形喷动床、三维狭缝式矩形截面锥底喷动床以及其双体联合的放大床内的若干重要问题进行了深入系统的研究，主要研究内容涉及模型的对比验证、床内的气固动力学特性、气固速度、流率及空隙率的分布特性、颗粒的混合和扩散、颗粒的循环时间、内构件的影响及其受到的磨损、喷动床的三维效应以喷动床的尺寸放大效应等。图5给出了双体喷动床的CFD-DEM数值模拟研究，颗粒数为250万。
详细内容请参考张科博士论文《复杂稠密气固两相流动的CFD-DEM模拟研究》。
房明明博士
以复杂稠密气固两相流动为研究对象，采用多相流研究中的欧拉-拉格朗日理论为气固流动描述的方法 论、以离散元耦合大涡模拟方法为主要研究手段，在对离散元-大涡模拟方法进行对比验证、对边界条件和模拟参数对模拟结果影响进行分析的基础上，针对概念上极为相似的两种典型的气固循环流态化装置——内循环床（ICFB）和循环流化床（CFB）内若干基础问题进行了较为全面和系统的研究。研究内容主要涉及气固物性对内循环床运行表现（主要以颗粒循环流率和气流短路通量来表征）的影响、相关设计和操作参数对内循环床运行表现的影响、内循环床颗粒混合机理及相关参数（设计、运行及颗粒物性等）对装置内混合过程的影响、腔室内置埋管对内循环床气固流动特征（气固速度及其通量、颗粒在各腔室内停留时间及整体循环时间等）的影响、内循环床中埋管自身的磨损规律、三维循环流化床中气固流动细节特征及提升管截面形状对装置内气固流动特征（固含率、气固速度及其通量的分布等）的影响。
图6 内循环流化床颗粒运动
详细内容请参考房明明博士论文《循环流化床内稠密气固两相流动的DEM-LES模拟研究》。
杨世亮博士
主要进行了以下3方面研究：
（1）研究了系统内密相流化床内气固两相间相互作用的若干规律，考察了两相中换热机理；并结合埋管周围局部气固流动特性，深入揭示了埋管周围传热系数的分布机理以及不同传热机制对管壁换热特性的影响。同时，阐释了管束对密相流动内固相运动的影响以及管壁磨损特性；
（2）研究并揭示了喷动系统内气固间运动及不同区域间相互作用强度，同时分析了不同床体几何结构和操作参数对系统性能的影响。
（3）从颗粒尺度揭示了内循环流化系统内固相循环特性以及其在不同区域内停留及耗散特性。同时，研究了采用双侧返料方式对循环流化床提升管内气固流动不均匀性的改善程度。
详细内容请参考杨世亮博士论文《流化床内稠密气固两相流动机理的CFD-DEM耦合研究》（获浙江大学2014-2015学年优秀博士学位论文提名奖）。
王帅博士
基于欧拉-拉格朗日理论框架，发展了适用于稠密气固两相反应流的大规模并行CFD-DEM计算平台，综合考虑了颗粒-颗粒/壁面碰撞、湍流、传热传质、辐射、颗粒缩核、热解、均相/异相反应、污染物生成及其之间的强烈耦合。采用该计算方法，分别对循环流化床内固相输运机理进行了深入研究。研究内容包括：宏观气固流动特性、传热特性、反应特性、污染物生成特性、固相停留特性、固相耗散特性、固相返混特性、颗粒聚团行为、操作参数和内构件对系统性能的影响等。
详细内容请参考王帅博士论文《流化床内稠密气固两相反应流的欧拉-拉格朗日数值模拟研究》。
胡陈枢博士
方法层面：系统地综述了稠密气固两相流模拟方法与相关模型的准确性、效率、适用范围。同时，基于碰撞、传热修正模型并耦合煤气化反应模型，首次将Coarse-grained CFD-DEM方法拓展到气固反应过程。
机理层面：揭示了加压鼓泡床内气泡驱动下埋管磨损机理，揭示了喷动床内恢复系数对气固流动行为的控制机制，系统性地考察了鼓泡流化床内流动、混合、传热与化反应的耦合机理，揭示了鼓泡床煤燃烧过程中局部过热现象的行程机理。
图9 生物质快速热解过程
详细内容请参见胡陈枢博士关于生物质快速热解的CFD-DEM模拟文章（Hu, C., Luo, K., Wang, S., Sun, L., & Fan, J. 2018,58(3), 1404-1416.），该文被《Industrial& Engineering Chemistry Research》选为封面文章刊登，同时被科技网站Advances in Engineering（https://advanceseng.com）报道，报道指出：In summary, the study by Zhejiang University researchers successfully presented the use of computational fluid dynamics/discrete element methods in the analysis of the biomass pyrolysis processes thus overcoming challenges witnessed in the previous techniques like two-fluid models. Generally, the influence of model parameters and different fluidized bed conditions are uncovered. Altogether, the proposed framework offers a promising solution in enhancing the use of biomass energy as a way ofreducing the dependence on fossils fuels, which is a key environment pollutant.
未完待续
后续将介绍稠密气固两相流的传热模型及反应模型
参考文献
[1] Rizk M A, Elghobashi S E. Atwo-equation turbulence model for dispersed dilute confined two-phase flows[J].International Journal of Multiphase Flow. 1989, 15: 119-133
[2] Shaffer F, Gopalan B, Breault RW, Cocco R, Karri S B R, Hays R, Knowlton T. High speed imaging of particleflow fields in CFB risers[J]. Powder Technology. 2013, 242: 86-99.
[3] Luo K, Tan J, Wang Z, Fan J.Particle-resolved direct numerical simulation of gas-solid dynamics inexperimental fluidized beds[J]. Aiche Journal. 2016, 62: 1917-1932.
[4] Tang Y, Lau Y M, Deen N G,Peters E A J F, Kuipers J A M. Direct numerical simulations and experiments ofa pseudo-2D gas-fluidized bed[J]. Chemical Engineering Science. 2016, 143:166-180.
[5] Clarke D A, Sederman A J,Gladden L F, Holland D J. Investigation of Void Fraction Schemes for Use withCFD-DEM Simulations of Fluidized Beds[J]. Industrial & EngineeringChemistry Research. 2018, 57: 3002-3013.
[6] Freireich B, Kodam M, WassgrenC. An exact method for determining local solid fractions in discrete elementmethod simulations[J]. Aiche Journal. 2010, 56: 3036-3048.
[7] Peng Z, Doroodchi E, Luo C,Moghtaderi B. Influence of void fraction calculation on fidelity of CFD‐DEM simulation of gas-solid bubbling fluidized beds[J].Aiche Journal. 2014, 60: 2000-2018.
[8] Sun R, Xiao H. Diffusion-basedcoarse graining in hybrid continuum–discrete solvers: Theoretical formulation and a priori tests[J].International Journal of Multiphase Flow. 2015, 77: 142-157.
[9] Wang Z, Teng Y, Liu M. Asemi-resolved CFD-DEM approach for particulate flows with kernel basedapproximation and Hilbert curve based searching strategy[J]. Journal ofComputational Physics. 2019.
[10] Amritkar A,Deb S, Tafti D. Efficient parallel CFD-DEM simulations using OpenMP[J]. Journalof Computational Physics. 2014, 256: 501-519.
[11] 王则力.多相流全尺度直接数值模拟方法及应用[D]. 杭州: 浙江大学, 2010.
[12] 谭骏华. 气固流化床内流动特性和颗粒曳力的直接数值模拟研究[D]. 杭州: 浙江大学, 2017.
[13] 桂南.复杂两相流动中颗粒碰撞的DEM-LES/DNS耦合模拟研究[D]. 浙江大学, 2010.
[14] 张科. 复杂稠密气固两相流动的CFD-DEM模拟研究[D]. 杭州: 浙江大学, 2012.
[15] 房明明. 循环流化床内稠密气固两相流动的DEM-LES模拟研究[D]. 杭州: 浙江大学, 2013.
[16] 杨世亮.流化床内稠密气固两相流动机理的CFD-DEM耦合研究[D]. 杭州: 浙江大学, 2014.
