摘要：本文提出了一种快速的离散颗粒模拟方法，即非解析LBM-DEM。流体相和颗粒相分别由格子Boltzmann法(LBM)和离散单元法(DEM)完成计算，并采用浸入运动边界法(IMB)实现流固耦合。通过将IMB方法中的额外碰撞项视为力源项并引入等效曳力来修正IMB中的权重函数以施加正确的耦合。此外，在GPU端实现了LBM-DEM以提升计算效率。基于非解析LBM-DEM对三种典型的气固流化进行模拟，均与实验数据吻合较好。同时，基于GPU的LBM-DEM相较于CPU或CPU-GPU异构的CFD-DEM算法，计算速度提升了1~2个数量级。该方法在分辨率尺度和计算效率上均优于传统的CFD-DEM方法，在流态化和多相流计算中具备广泛应用前景。SIMPOP引言颗粒流体系统广泛存在于自然界及工业生产中，如石油催化裂化、循环流化床燃烧、煤气化、流化焙烧和沙尘暴等。传统的双流体模型TFM将流体相和颗粒相均视作连续介质求解，无法获得颗粒的详尽信息；直接数值模拟DNS，令流体网格比固体网格小一个数量级，直接解析颗粒受力并描述每个颗粒周围的详细流场，但其计算量巨大，难以实现对真实的流动现象进行模拟。离散颗粒模拟DPS介于双流体模型和直接数值模拟之间，其不解析颗粒的受力，而是基于曳力方程关联流体和颗粒，不仅可以获得颗粒运动轨迹等局部信息，也可以获得流体的平均流场信息。其采用的流体网格可比颗粒尺寸大一个数量级，能在计算精度和成本之间达到一个平衡，已成为探索实验室规模气固系统的有力手段。DPS下最常用的方案为CFD-DEM，流体相基于Gidaspow提出的A类和B类模型(体积平均N-S方程的两种不同形式)，采用SIMPLE算法或基于Fluent，MFIX和OpenFOAM等软件平台。颗粒相普遍采用基于软球模型的DEM(Discreteelementmethd)进行更新。在计算层面，CFD-DEM中流体相的求解多为隐式或半隐式算法，并行上存在巨大困难，为了降低CPU的负载，提高计算效率，目前比较流行的方法是将可并行的DEM部分使用GPU计算，而难以并行的流体部分依然采用CPU计算，即CPU-GPU异构算法。然而CPU和GPU并不共用内存，每一步下不同设备间存在复杂的异构通信，尽管较纯CPU计算已取得显著加速，但效率依然较低。格子玻尔兹曼方法LBM（LatticeBoltzmannmethod）是近年来发展的一种新的计算流体力学方法，由于自身的高并行度，非常适合采用GPU计算。目前LBM在颗粒两相流中的应用主要集中在颗粒解析的DNS层次，应用LBM实现颗粒非解析的气固两相流模拟目前鲜有报道。本文提出了一种新的颗粒非解析LBM-DEM算法。其求解方程基于求解LBM-DNS中IMB的形式，但不再采用原IMB方法计算权重函数，而是将IMB方程中用于恢复无滑移条件的额外碰撞项视作力源项，通过定义格点曳力获得权重函数以修正LBM方程。同时，鉴于LBM和DEM均具备较好的并行性，本文在GPU上对该算法进行实现以发挥其最大优势。纯GPU算法避免了CPU-GPU间复杂的异构通信，较目前算法有着巨大的性能提升。LBM-DEM方法01气体运动——格子玻尔兹曼法标准的格子Boltzmann方法(LatticeBoltzmannmethod,LBM)将密度分解为多个方向离散的分布函数，并在格点网格内将各分布函数碰撞迁移，通过对当前格点上各分布函数的组合还原出密度、速度等宏观量。其方程形式为：其中为t时刻位于x格点的第i个分布函数，为格子速度，三维下一般采用19个方向的分布函数，即D3Q19模型，如图1所示。为松弛因子，由粘度确定，为平衡分布函数，由当前格点的密度和速度确定：式中为运动粘度，为密度，u为速度，为格子声速，，c为空间步长与时间步长大小之比，当空间步长与时间步长相同时。为权重系数。流体速度、压力和密度可由下式更新：02颗粒运动——DEM方法单颗粒运动遵循牛顿第二定律：其中为颗粒质量,为颗粒速度，为颗粒受到的碰撞力，为颗粒受到的曳力，为浮力,为重力加速度。在DEM方法下，基于弹簧阻尼假设的软球模型获得颗粒间碰撞力，如图2所示。这里采用简化DEM模型，仅考虑颗粒间法向力计算如下：其中e为弹性恢复系数，描述颗粒碰撞后的能量耗散程度；分别为第i个颗粒的位置、半径和质量。i和j分别为相互碰撞的颗粒i和j，为法向刚性系数，为重叠量，为单位法向矢量,为法向阻尼系数，为颗粒相对速度的法向分量。03流固耦合流体和颗粒间的耦合通过颗粒所受曳力实现。单颗粒的受力可基于Wen-Yu曳力描述，其形式为:其中为颗粒雷诺数,为滑移速度，。r为颗粒半径，为颗粒位置处的流体速度。为颗粒位置处的固相体积分数。浸入运动边界法引入了格子控制体和格子固含率的概念，即计算以流体格点为中心的立方单元内固体体积所占比重，图3为以xy面为例的示意。基于非平衡态反弹和流固分相的思想，在标准单松弛的LBM中引入额外碰撞项描述固体作用，则公式(1)修正为：式中为格点上的固体速度，为格点上的流体速度，即格子固含率，-i代表与i反向。为权重函数，控制流体碰撞算子和额外碰撞项的计算比重，满足。为1时完全由固体控制，公式（15）还原为标准反弹格式，为0时完全流体控制，方程还原为标准的单松弛LB方程。求解DNS问题时，传统IMB的和网格内颗粒对流体作用力形式为：经典的浸入运动边界法在描述颗粒解析的气固两相流系统中取得了重要成功，由于DNS系统下流体网格比颗粒网格小一个数量级，加和颗粒内所有的流体格点的即为颗粒所受合力，则牛顿第三定律天然满足，且基于固含率实现对边界的平滑描述。在颗粒非解析的模拟系统中，颗粒尺寸比流体网格小一个量级，基于曳力公式推导的颗粒曳力和LB方程下推导的流固耦合作用力（18），其计算过程相对独立导致二者大小并不匹配，需要进行统一。本文将公式（15）中的额外碰撞项近似视为力源项，提出一个新的流固耦合方案如下，由图3(b)所示，假设格点内存在一个等效颗粒球，其速度为，定义等效粒径由下式给出：其中为网格体积，可为1。则基于和，通过公式(14)或(15)-(18)计算出等效颗粒曳力，令LB获得的流固作用力与该等效曳力相等，即，则由公式（18）得出形式为：通过上述公式取代公式（17），即可对流场施加流固作用力，同时消除了传统IMB方法计算受力时，权重函数存在粘度项导致粘度依赖的缺陷。GPU实现LBM-DEM整体求解框架如图4所示,通过交互各自的映射信息实现耦合。在流场和颗粒场均初始化完毕后，将数据由CPU送入GPU执行运算。首先DEM部分，在GPU端开辟与颗粒数量相同的线程，每个线程独立计算单个颗粒以同时更新颗粒的速度和位置，随后执行颗粒排序并重构邻居列表，执行颗粒碰撞。进一步，将离散的颗粒速度和固体体积映射到流体网格，同时获得流体网格信息以求解颗粒所受曳力，最后将颗粒速度修正到t+1时刻。LBM下，在GPU端开辟与流体格点数量相同的线程，每个线程独立更新一个流体格点。流体获得格点固体速度和固含率后，将流体速度更新到t+1时刻，如此循环直到结束完成计算。整个过程全部在GPU上完成，CPU仅负责调用和后处理数据，避免了传统CFD-DEM方法中，CFD和DEM分别在CPU和GPU求解，二者耦合带来的复杂异构通信，因此具备极高的效率。LBM-DEM准确性验证01鼓泡流化首先，通过鼓泡流化床与CFD-DEM的结果进行对比。这里基于Lu的CFD-DEM算例，其采用CPU基的MFIX求解器完成流场的计算，基于自主实现的GPU求解器完成DEM的计算并完成二者耦合。本节基于完全GPU实现的LBM-DEM求解器求解该算例，所有物理参数与Lu的设置完全相同。速度入口采用标准平衡分布格式，壁面为标准反弹格式，出口为一阶外插，每两个格子步长更新一次颗粒邻居列表，曳力采用Wen-Yu形式。物理量对应格子单位下的无量纲参数如表1所示。这里对比了本文的耦合方案以及WANG等人提出的不修正IMB耦合方案。不同时刻下的瞬态颗粒分布如图5所示。由图5(a),t=0.05s时顶部颗粒开始向上运动，两侧和底部颗粒由于气速较小颗粒向下沉降；t=0.1s时达到一定高度后，除顶部颗粒继续向上运动外，中部颗粒开始向下沉降，非均匀结构即将形成；t=1.0s后系统基本达到了动态平衡，体系内同时存在着向上和向下运动的颗粒聚团，颗粒分布与Lu的结果定性吻合。基于WANG等人提出的传统IMB耦合方案结果如图5(b)所示，各个时刻下颗粒始终沉降堆积在床层底部，并没有出现流化现象。图6为不同耦合方案下计算的瞬态压降曲线，与Lu的结果对比。原始的修正方案下，由于颗粒没有流化，床层压降远低于文献结果，表明传统的耦合方式在此粘度下已不再适用，而基于本文提出的耦合方案，床层压降在t=1s后保持了动态稳定，平均压降为8.48Pa，相较于CFD-DEM结果(8.37Pa)误差仅1.31%，表明本文提出的LBM-DEM算法可以得到准确的动力学结果，消除了原LBM-DEM的局限。02流域转变进一步，基于微流化床中的流域转变与实验结果进行对比验证。模拟的物理尺寸3mm*3mm*27mm，模拟区域底部为非平衡外推速度边界，顶部为一阶外插出口，其余部分为壁面，采用标准反弹格式。选用=53μm的FCC颗粒，颗粒总数为443000，初始床层空隙率约为0.37，其余参数详见表2。不同入口速度下的床层孔隙率和对应速度下中心截面的时均固含率分布如图7所示。床层孔隙率由床层膨胀高度推导而来，可以较好验证气固算法的准确性。首先将达到动态稳定后的轴向固含率进行时均化处理，取固含率为0.2时的高度作为床层膨胀高度H，则床层孔隙率为：.其中A为底面积，为单颗粒体积。由图7，不同气速下的床层孔隙率均与实验值相吻合。其中和为实验测定的最小流化速度和最小鼓泡速度。当入口速度大于时，床层空隙率开始快速增长，随着气速进一步增大，空隙率的提升量逐渐减小。曲线平均误差为4.13%，最大误差在=5.69mm/s处为9.01%，均满足误差低于10%的要求。由图7，当入口速度在和之间的范围，流动应为均匀的散式流化。由(III)当气速超过后，系统内产生了均匀上升的气泡，表明LBM-DEM模拟的流动结构与实验测定的最小流化速度和鼓泡速度范围基本吻合。随着气速进一步提高，由(IV)流动结构已经由鼓泡流化转向节涌流化；随着气速增大，颗粒运动更加剧烈。由(VI)，流型已出现向湍动流化过渡的趋势。本文提出的LBM-DEM方法复现了颗粒群在微流化床中的流域转变，定性和定量结论均与实验结果一致，表明该算法可作为模拟气固两相流动的可行方案。03快速流化基于该方法对快速流化进行模拟研究。为了更好展现颗粒运动时的非均匀结构，这里基于三维矩形薄板复现快速流化现象。针对流场，底部为速度入口，顶部一阶外插出口，x法向面为标准反弹边界，y法向面为周期边界。对于颗粒，颗粒沿y方向和z方向运动均做周期运动。模拟使用的总颗粒数为2212848，初始固含率约为0.09，采用Wen-Yu曳力形式。模拟的相关参数如表3所示。应用该算法得到了超高解析度的气体-颗粒分布，如图8所示。t=1s后流化床内颗粒结构达到完全演化，整个床层都分布着颗粒聚团，同时床层内存在明显的稀密相结构，密相下颗粒的浓度较高，对应气体速度较低，集中在底部和边壁，而稀相内颗粒浓度较低对应气体的速度较高。颗粒整体分布呈现上稀下浓，中间稀两侧浓的非均匀结构，与Li等人的模拟和实验结果定性一致。基于t=2s后的时均结果得到快速流化下轴向固含率如图9所示。颗粒浓度沿轴向分布为一条S型曲线，与LiandKwauk在10m床高中的实验规律基本吻合，表明了该算法在分析快速流化问题中的适用性。04LBM-DEM速度测试在验证准确性的基础上，进一步检验计算效率。首先对比01节的算例，与Lu报道的计算时间对比。Lu分别测试了CPU程序和CPU-GPU异构计算的CFD-DEM算法，CPU算法基于MFIX实现，使用10核心CPU并行计算，计算时间记为MFIX-CPUDEM;CPU-GPU算法基于MFIX和自主编程的DEMGPU程序实现，同样使用10个CPU核计算流场，使用单张TeslaP100GPU求解颗粒场,计算时间记为MFIX-GPUDEM。本文所有算例使用的计算设备均为单张TeslaV100GPU完成流场和颗粒场的计算，计算时间记为GPULBM-DEM，模拟的物理时间均为5s，结果如图10所示。由Lu的测试，纯CPU计算总时长为76.25h,其中DEM计算时间占总时长的91.2%，运算时长基本取决于颗粒相的计算；通过将DEM转移至GPU计算，总时长缩短至5.17h，相较于纯CPU算法加速比达到14.7，其中DEM部分计算时间占比17.2%，计算瓶颈由颗粒相转移到流体相。本算法流体相和颗粒相均在GPU求解，总计算时长仅为0.74h，相较于CPU-GPU异构实现，带来了6.98的加速比，相较于纯CPU加速比达到了103。其中流体相的计算时间被大幅缩短，得益于LBM的高并行度，流体部分的计算时长由4.25h缩短为0.03h，仅占总时长的4.05%，总时间重新取决于求解颗粒相的用时。进一步，基于本文前三节的算例对该算法进行速度测试，依次命名为caseA(鼓泡流化，颗粒数12,8000)、caseB(流域转变，颗粒数44,3000)和caseC(快速流化，颗粒数2,212,848)，同时将caseC中的流体网格和颗粒规模沿xyz方向均扩大2倍定义为caseD(颗粒数17,702,784)。类比MLUPS定义MPUPS作为DEM部分的速度指标，即每s更新的百万颗粒数，计算公式为，其中Step为更新总步数。统计物理时间1s下LBM-DEM的计算总时长以及LBM部分和DEM部分的MLUPS和MPUPS，结果如图11所示。caseC中包含160万以上的流体网格和220万以上的颗粒，使用该算法仅需2.071h即可完成物理时间1s的更新，展现了极高的计算效率。caseD计算1s的时间为16.25h，为caseC的7.85倍，与其计算规模较caseC扩大8倍基本相当，表明该算法的计算时间是随计算规模线性增长的，增大流体网格和颗粒数量并不会降低程序的计算速度。随着颗粒规模的增大，MPUPS持续增加并逐渐达到峰值速度，caseD下MPUPS为157.74。基于MLUPS和MPUPS可在已知求解规模下大体预估程序的整体求解时间。本文提出的LBM-DEM具备远超传统CFD-DEM的计算效率，较纯CPU算法可带来两个数量级的加速提升，较CPU-GPU异构算法也能实现数倍加速，非常适合求解大规模气固颗粒系统。尽管计算设备存在差异不能单纯比较计算时间，但本文提出的LBM-DEM作为纯GPU实现的气固两相流算法，未来随着GPU硬件更新将具备更加显著的优势和应用前景。结论本文在LBM和DEM的基础上，提出了一种新的未解析LBM-DEM方法。该方法采用修正的IMB进行流固耦合，成功将IMB从DNS系统扩展到DPS系统。与原始的加权函数相比，新的耦合模型显著提高了气固流动模拟的精度。进一步，利用LBM和DEM固有的并行性，成功地实现了用于未解析LBM-DEM的GPU框架。在鼓泡流化、流域转变、快速流化下，利用新的耦合模型均复现了流化现象且定性定量结论与实验基本吻合。在相同规模下，与传统CFD-DEM上的CPU和CPU-GPU异构算法相比，加速比分别为103和6.98，展现了极高的执行效率。作为一种通用的求解框架，基于GPU的非解析LBM-DEM在求解尺度和计算效率方面较传统的CFD-DEM具备显著的优势，未来将具备广泛的应用前景。本文翻译自Fu,S.,Wang,L.,2023.GPU-basedunresolvedLBM-DEMforfastsimulationofgas-solidflows.ChemicalEngineeringJournal465,142898.LMFD（Lattice-basedMulti-FluidsDynamics）是面向多相流体系大规模数值模拟的科研和工程软件，由中国科学院过程工程研究所EMMS团队开发。该软件基于格子玻尔兹曼方法开发，结合自研的耦合算法和计算模块，能够处理包括单相流、气液两相、气固两相以及气液固三相等复杂多相流问题。同时软件设计了全过程GPU计算方案，具备天然的高并行优势，能够利用GPU计算资源实现大规模问题的高效计算。来源：多相流在线
。 
图3 液滴流流场情况，(a)图为数值模拟结果；(b)图为实验结果 
