流化床内流动特性的DEM数值模拟研究

摘要流化床燃烧技术作为一种清洁煤燃烧技术，已经被广泛应用于生产。但由于流化床内流动的复杂性，人们对其工作机理并没有清楚全面的认识。为更深入研究流化床气固流动机理，本文采用数值模拟方法对不同射流孔数流化床流动特性进行模拟。气固两相流数值模拟技术可以分为两大类：即欧拉—欧拉颗粒拟流体模型和欧拉—拉格朗日离散颗粒模型。本文采用欧拉—拉格朗日模型从颗粒水平上建立流化床内气固流动数学模型，并对流化床内流动过程进行模拟,对气固流动特性进行了分析。首先采用离散单元法对床内流动特性进行了数值模拟，分别模拟得到单孔、双孔分布板中气泡的形成、上升、破裂的过程，同时给出了不同入射流速对流动特性的影响。模拟结果表明：颗粒随气体的运动而运动，气体运动是颗粒运动的主要动力源，相对单喷口系统，由于相邻气泡之间存在横向和纵向双方向的的聚并，挤压，破裂等现象，因此双喷口系统中床层内颗粒的运动和扩散作用也更加强烈。随后，在单双喷口模拟的基础上，对多喷口射流流化床系统进行模拟，通过分析喷口个数、喷口位置、进气速度、颗粒密度、弹性系数等对气固流动的影响，对比单双喷口系统，认为进气孔数的增加时，弹性系数是颗粒之间的“互锁”现象的主要影响因子；随着进气孔数的增大，颗粒与气体的混合作用也就越强烈，流化床进入充分流化阶段的时间也越短；同时，喷气孔密集程度对气固垂直、水平速度也有较大的影响，速度随密集程度的增大而增大，气固流动程度随密集程度的增大而加深。关键词：循环流化床；流动特性；多孔；数值模拟华北电力大学硕士学位论文AbstractAsanewclean-burningtechnology,circulatingfluidizedbed(CFB)hasbeingwidelyusedintheproduction.Butbecauseofitscomplexnature,peopledon’thaveclearandcomprehensiveunderstandingoffluidizationmechanism.Togetamorein-depthstudyofthemechanismoffluidizedbed,hydrodynamicscharacteristicsoffluidizedbedwithdifferentnozzlecountwasresearchedwithnumericalsimulationmethodinthispaper.Numericalsimulationtechnologyfordensegas-solidflowcouldbedividedintotwocategories:theEuler-Eulerparticlesquasi-fluidmodelandEuler-Lagrangiandiscreteparticlemodel.Inthispaper,thediscreteelementmethodwasused,adescriptionofthefluidizedbedgas-solidflowwasbuild,andsimulatedtheprogressofgas-soildflow,analysisthedynamiccharacteristics.First，discreteelementmethodwasproposedtoimplementnumericalsimulationofflowcharacteristicsinfluidizedbed,throughthesimulation,theprograssofconformation,ascensionandburstofbubblewasobtained,theeffectofflowcharacteristicsindifferententrancevelocitywasalsorevealed.Numericalsimulationresultsshowthatparticlemovewithgas,themotionofgaswasthemainpowersourceofparticle,relativetosingle-nozzlesystem,becauseofthemergence,extrusionandburstbetweencongenialbubbles,themotionandthediffusionwasmorefierceindouble-nozzlessystem.Secondly,onthebasisofthesingleanddoublenozzlessimulation,themulti-nozzlefluidizedbedsystemwassimulated,throughtheanalysisoftheimpactofthenozzlenumber,thenozzleposition,thegasvelocity,particledensity,thecoefficientofelasticityofsolidflow,contrastwithsingleanddoublenozzlesystem,withtheincreaseinthenumberofnozzle,thecoefficientofelasticitywasthemainfactorsaffectingthephenomenonof"interlock"betweenparticles;withtheincreaseofthenumberofinletholes,thetimeofchangingintofullflowstagewasshorter;also,theintensiveofnozzlesimpactedtheverticalandthehorizontalvelocityoftheparticlesandgas,thedegreeofgas-solidflowintensityandthevelocityincreasewiththeincreaseofintensive.Keywords:circulatingfluidizedbed;dynamiccharacteristics;multi-nozzle;numericalsimulation目录目录摘要 IAbstract II目录 III第1章绪论 51.1选题背景及意义 51.2气固流态化系统的组成及其分类 61.3气固流态化的研究方法 61.3.1两相模型 71.3.2欧拉方法的颗粒相拟流体模型 71.3.3欧拉-拉格朗日法的颗粒轨道模型 81.4流化床DEM法研究进展 91.5本文主要研究内容 11第2章流化床气固两相流动DEM模拟数学模型 122.1离散单元法原理 122.1.1时间步长的确定 122.1.2颗粒运动控制方程 132.1.3颗粒碰撞力学模型 142.1.4碰撞对象的搜索算法 152.2流体相数学模型及求解方法 162.2.1流体相数学模型 162.2.2流体相控制方程组的求解 192.3气固两相相互作用力分析 232.3.1气体对颗粒的作用力 232.3.2颗粒对流体的反作用力 252.4本章小结 26第3章循环流化床气固流动特性的数值模拟 273.1自由堆积过程模拟 273.2单孔射流流化床流动过程模拟 303.3双孔射流流化床流动过程模拟 34第4章多喷口流化床气固流动特性模拟 394.1三孔射流模拟对象及参数设定 394.1.1弹性系数对流动情况的影响模拟 404.1.2不对称喷口布置对流动过程影响模拟 454.2四孔射流流化床模拟 504.2.1模拟对象及参数设定 504.2.2四孔射流模拟结果 514.3本章小结 54第5章全文总结及展望 565.1全文工作总结 565.2下一步工作展望 56参考文献 58攻读硕士学位期间发表的学术论文及其它成果 62华北电力大学硕士学位论文第1章绪论1.1选题背景及意义在“绿色发展，建设资源节约型、环境友好型社会”这一“十二五”规划的要求与倡导下，作为世界上最大的能源消费国之一，能源与环境问题成为我国所必须面对的重要问题，在我国的一次能源消费比例中，煤炭资源高达70%左右，预计到2050年，这一比例依然会在50%以上，可以说，煤炭资源是我国最主要的能源，但我国煤炭的重要特点之一是高硫煤所占比例较高，导致燃烧大量产生SO2，2010年全国第一次污染源普查公告显示，电力热力的生产和供应业SO2排放量为1068.70万吨，NOX排放量为733.38万吨[1]。其中SO2会导致呼吸系统疾病，与大气中水分结合生成SO3并形成酸雾或酸雨，对建筑物，植物等具有很强的腐蚀性，而NOX则是酸沉降及光化学烟雾的主要形成因素之一，为解决这一问题，实现能源节约与环境保护，循环流化床以其强燃料适应性和低污染物排放等特点为广大学者广泛研究。作为一种比较成熟的清洁燃烧技术，循环流化床燃烧技术具有氮氧化物排放低、低温燃烧、可实现燃烧过程直接脱硫、燃烧强度大、负荷调节性能好、灰渣便于综合利用[2]等优势，已成为煤炭洁净燃烧的首选炉型。然而循环流化床锅炉也存在着磨损严重、排渣系统出力不足、达不到满负荷出力、炉膛温度偏高导致脱硫效率降低等问题，为研究解决这些问题，有必要对流化床运行过程中的内在机理、具体的流动燃烧过程进行全面的分析研究，而传统的试验方法只能体现流化床宏观的一些特性，不能细致反映瞬时结构和动态过程，改进后的试验方法如高速摄像、图像处理、示踪离子等也由于条件的局限性和结果的不完全性，无法对流动燃烧过程有一个全面透彻的解释，鉴于上述情况，研究者着手从理论层面上建立数学模型，用数值方法详细全面分析循环流化床内的气固流动、燃烧、传热特性。作为循环流化床技术的理论基础，气固流态化技术是研究气固两相之间伴随有热量、质量传递的气固相互作用的一门工程科学，由于影响气固两相流动因素的复杂性，尤其流化床中散布相互作用的高浓度颗粒、湍流气流，同时又涉及流体力学、热力学、燃烧学、传热传质学等相关基础学科，促成气固流动的作用机理还没有被研究者、完全认知，流态化理论也成为热门研究方向之一。1.2气固流态化系统的组成及其分类气固流态化系统的结构通常是由容器及其底部的布风板构成，流体从底部的布风板进入容器，流过并推动容器中固体颗粒，从而使其具有一定的流体性质。由气体流速的不同，固体颗粒的流动状态可以分为以下几类：当气流速度较低，流体作用于颗粒之上的曳力小于固体颗粒本身的重力，颗粒不会被流体曳力推动，该种状态称为固定床状态(fixedbed)；气体流速增大，气体曳力开始大于颗粒重力（此时的流体速度称为最小流化速度），这时颗粒开始随流体蠕动，颗粒进入流态化状态，床层显示出类似流体的性质，该种状态称为散式流态化状态；气体流速继续增大，达到最小鼓泡速度，气体以气泡的形式穿过床层，随着气泡的产生、并聚与破裂，此时整个床层分为两种相态：颗粒集中的乳化相以及气体浓度较高的气泡相，该种状态称为鼓泡流态化状态(bubblingfluidizedbed)，如果容器较小，气泡直径增大至撑满整个容器横截面而形成节涌(slugging)；气体流速进一步增大，大气泡由于表面张力及挤压作用而破碎，气泡尺寸变小，整个容器内流动状态变的比鼓泡床状态更加剧烈，床层边缘变的模糊不清，此时该种状态称为湍流流态化状态(turbulentfluidization)；在湍流流化状态下继续提高风速，床层表面变的更加模糊，气固流动也发生变化，气相由湍流流动时的离散相（气泡）变为连续相，颗粒由连续相（乳态化状态）变为离散相，气固流动进入快速流化床状态(fastfluidizedbed)。继续增大风速，床层颗粒浓度逐渐变稀，气固流动将进入气力输送状态(pneumatictransport)。1.3气固流态化的研究方法随着流态化技术应用的日益广泛，众多国内外研究者对流态化过程做出了大量的试验和理论研究，如近年来有很大发展的快速信号处理技术及非接触式的测量技术如过程层析成像技术、颗粒动态分析仪测量技术、激光多普勒测速仪测量技术等，这些技术为气固流动系统中参数瞬态值的测量提供了有效的手段，但仍不能对稠密两相流中的基本物理量如颗粒浓度、压力及颗粒速度等进行较为精准的测量。为了便于更加深入、细致的认识流态化过程，建立较为系统完整的理论，研究者把目光投向了近来来发展迅速的微机计算技术，研究者以动量方程、质量方程、能量方程以及组分守恒方程为基础，通过数值计算方法建立了一系列的流态化过程研究理论。目前，流态化过程的研究模型可以分为三大类：第一类模型是从流态化过程分析，从流动结构的角度出发建立的两相模型；第二类模型是以欧拉方法为基础的颗粒相拟流体模型；第三类模型是以拉格朗日方法为基础的颗粒轨道模型。1.3.1两相模型两相模型有其独特的优点,如能够根据一些局部信息合理的估计整个流化床的流动特性，比如气固流动形式，各相速度及气固两相混合程度等。两相模型理论是由ToomeyandJohnstone[3]在1952年首先提出。Davisdon和Harrison则在1969年提出了著名的Davidson模型从而解释了上升气泡周围气体以及颗粒的运动规律与压力分布[4]。在以远离气泡的区域压力梯度无扰动及气泡中的压力不随气泡运动而改变这两个边界条件为前提下提出如下理论：气泡运动挤压颗粒团时，颗粒团会像一般流体一样，随气泡的挤压而向两侧流动，同时，气固两相流动相对速度满足Darcy定律。D.K.unni经过研究又在Davidson模型的基础上提出了K-L模型[5]。李静海等人则建立了多尺度能量最小模型，给出了两相流流域转化和不同操作模式的存在条件，并从机理上揭示流域转化特征[6]。由于实验关联式是两相模型的建立必须组成部分，而其精度较低，其应用范围也有一定的限制，因此两相模型只能描述流态化系统的整体运动行为，在分析局部动态特性时，两相模型往往不能达到要求。1.3.2欧拉方法的颗粒相拟流体模型从根本上说，流化床中的气固流动属于稠密的气固两相流范畴，与一般的稀疏气固两相流系统不同，稠密气固两相流中，由于颗粒的质量、体积所占份额较大，颗粒对气体流动的影响，颗粒之间的碰撞作用对整个系统流动特性的影响都比较明显，因此研究一般两相流系统常用的模型对流化床系统不再适用。为了更合理的描述流化床系统的气固流动规律，研究者提出了一些新的模型，其中，颗粒拟流体模型是目前应用较广的一种模型，其原理为：在将流体视为连续相的同时，把颗粒也视为一种类似流态的连续介质，即拟流体。在模型中，根据守恒定理建立气固相的动量控制方程组，空间各点两相共用，两相相互渗透并存在着相互作用。颗粒相拟流体的概念最早是由S.L.Soo[7]等人在20世纪60年代末提出，现已经被很多学者应用到流态化现象的数值模拟中，其中Gidaspow等人预测了鼓泡床内气泡的形成、长大及破裂过程，预测结果与实际实验对比，可信度较高，他们通过计算得到了床内各相的速度分布[8-10]。同时，他们模拟了PYROFLOW型循环流化床锅炉的简化模型，得到了床内各相的速度分布，并得到了固体颗粒相各时刻的浓度场分布[11,12]。Gidaspow等人在计算流化床内壁面与床层之间的传热特性时也成功引入了拟流体模型[13,14]。Kuipers等人采用拟流体模型研究了三维流化床的流动特性[15]，Enwald等人全面的总结了拟流体在气固流态化方面的应用[16]。虽然利用拟流体模型在模拟流态化过程中得到了广泛的应用，也取得了很多的成果，但该模型本身也存在着一些难以解决的问题，比如，将颗粒相处理为连续相后如何确定拟流体的粘性应力和静压力；由于颗粒相的连续性处理，也就导致无法在微观层次上对颗粒进行研究，只能得到多个离散颗粒在局部区域的平均特性。1.3.3欧拉-拉格朗日法的颗粒轨道模型在颗粒轨道模型中，流态相仍然看作连续介质，离散颗粒相就作为离散体处理。早期的颗粒轨道模型只考虑流态对离散颗粒的作用，而颗粒之间以及颗粒对流体相的作用并不考虑，也因此单颗粒动力学模型只能在稀相气固两相流中适用，后经研究者的改进，使之能充分考虑颗粒对颗粒及流体的作用，尤其随着微机计算技术的飞速发展，目前属于颗粒轨道模型范畴的离散单元法(Discreteelementmethod,DEM)越来越多的被应用于研究流态化系统。DEM方法最早在1979年由Cundall等人提出[17]，1992年，Tsuji等人在首次成功应用DEM方法于浓相水平管道气力输送研究中[18]，随后又把DEM方法与流体力学模型相结合，成功模拟了流化床内气固两相流动过程[19]。当前DEM法已经成为流化系统模拟中一种重要的方法。DEM方法的主要思路为：用离散项考虑流化系统中颗粒，通过分析每个颗粒的受力判断计算颗粒的速度、加速度及位移情况，其同时考虑了颗粒间的作用力及颗粒对流体的作用力，所以DEM模型非常适合于模拟流化床这样的稠密气固两相流动系统。由于考虑了颗粒之间的相互作用力，所以DEM方法在处理颗粒碰撞时需要考虑碰撞时颗粒的形变，按假设有无形变，DEM法在处理颗粒相互作用力时分为硬球(hardsphere)软球(softsphere)两种模型。硬球模型中，假设颗粒不会因碰撞变形，视为刚性球体，碰撞在瞬间发生，颗粒之间的相互作用力按照动量守恒确定，与软球模型相比，硬球模型可以采用较大的时间步长，而且可以处理颗粒较多的流态化系统，其不足是不能模拟三个以上颗粒的碰撞，在处理碰撞的传热问题时由于颗粒内部导热热阻及颗粒间静止气膜的存在也会遇到很大的困难。与硬球模型不同，软球模型类似于阻尼-弹簧模型，其假设颗粒碰撞过程中发生形变并持续一定时间，当选定合适的时间步长后，软球模型可以对每一个颗粒的受力、运动状况进行跟踪，不足是由于要跟踪系统中每一个颗粒的受力与运动情况，并且时间步长不能取的太大，所以计算量很大。由于本文研究对象流化床属于稠密型气固流态化体系，所以采用软球模型。1.4流化床DEM法研究进展DEM模型在Cundall提出后，最早是由Y.Tsuji等人用于二维鼓泡流化床的数值模拟中[19]，其只在计算流体对颗粒曳力时考虑流体粘性，采用SIMPLE算法[20]求解N-S方程，成功模拟了床内气泡的形成、成长上升及破裂的过程及床层的压力波动曲线。JieOuyang等人采用DEM硬球模型对流化床进行了较为全面的模拟，他们通过二维鼓泡床中进行模拟，得到了气泡的生成、变大和破裂的过程并在一个小直径的鼓泡床中模拟观察了颗粒的节涌现象[21]，得到了流化床中模拟颗粒团的动态图象，同时也研究了操作条件对循环流化床颗粒团结构和“噎塞”行为(chokingbehavior)的影响[22]，通过模拟上升段内气固流动，得到了颗粒相运动的动态图象以及颗粒运动的无规则轨迹[23]。周浩生等人利用DEM硬球模型模拟了鼓泡流化床内气固两相流动特性[24]用颗粒粒径不同，平均粒径为3.04e-4m，颗粒数为420个，模拟得到了密度相同但直径不同的颗粒在颗粒在床内的流动情况:小颗粒集中在床的上方，大的颗粒则集中在床的下方，同时，也验证了颗粒运动时的B.H.Xu等人利用DEM软球模型对准三维鼓泡流化床内的颗粒运动特性进行模拟[25]，文中在气相处理中忽略了气体湍流影响，但考虑了层流粘性的影响，作者模拟得出了流化床最小流化风速，得到了鼓泡床床层压降与气体流速的关系，通过分析提出了床层压降与气体流速之间存在滞后现象。DegangRong等人利用DEM软球模型对一个有浸没水平管束的鼓泡流化床进行研究[26]，得到了流化速度及管道排列方式对床内气固流动运动规律的影响，并根据管壁附近颗粒运动状态分析了管道壁面的磨损规律。Mikami等人利用DEM软球模型模拟了粘结性颗粒的二维流化床[27]，文中考虑了颗粒表面所带水分产生的液桥力，模拟得到了湿颗粒床层的流动特性，与干燥颗粒条件下的模拟作比较，湿颗粒条件下的压力波动、最小流化风速的值要高。周浩生等利用DEM软球模型对系统为沙粒和煤粒组成的二元流化床体系进行模拟[28]，利用传热温差及相关传热系数计算式计算得到碰撞传热量，提出气固流动特性对颗粒的热特性有较大的影响，同时，研究了传热温差、流化速度、煤粒自身物理性质对煤粒热特性的影响。陆继东等利用DEM软球模型建立了流化床内气固流动、传热及煤粒燃烧的数学模型，将流化床内煤粒的加热、燃烧过程作为一个连续的过程进行了数值模拟[29]，研究了燃烧热、对流热、辐射热及碰撞热这几种不同的热传递方式在煤粒加热和燃烧过程中的贡献，同时研究了颗粒刚度、旋转性及表面摩擦系数对颗粒碰撞率的影响。李晓光等人利用DEM软球模型对二维高温流化床进行研究[30]，所选颗粒直径为3mm大直径颗粒，分别模拟了温度为300K与1000K时床内临界流化风速、压降脉动及床层最大膨胀比，对比发现温度较高时，临界流化风速也相当较大，运行较平稳，床层膨胀比明显降低，并指出高温大颗粒流化床内气泡快速成长同时气泡内部存在颗粒淋落现象。刘向军等人利用DEM软球模型研究了稠密气固两相流中颗粒团运动规律[31]，通过模拟得到了流化床内颗粒团详细运动碰撞的过程及浓度、粒径分布。黎明等人利用DEM软球模型对二维鼓泡床进行研究[32]，文中考虑颗粒与流态间的相互耦合关系，得到了不同入口气速下流化床随时间变化的状态图。田凤国等人对非均匀布风方式流化床的气固流动特性进行DEM数值模拟[33]，研究认为：颗粒在床内存在剧烈的循环流动；气泡在进气口处活动较为活跃，由于气泡的携带牵引颗粒向上运动，同时又由于重力作用颗粒向下运动，两者构成了床内气固流动内循环；风速对颗粒流量变化具有显著的影响；其后又对流化床内物料微观混合过程进行研究[34]，通过引入定量描述混合发张程度的Ashton指数，得到示踪颗粒场变化过程，研究发现：在均匀布风情况下，床内气泡横向运动受到限制，横向混合以扩散方式为主；而对于非均匀布风流化床，颗粒团具有较大的横向浓度梯度、速度梯度。张勇等人利用DEM软球模型对喷动流化床内颗粒混合特性进行模拟[35,36]，得到了床中上下两层颗粒的运动和混合过程，提出了表征颗粒混合质量的特征参数，观察了喷动气速和流化气速对颗粒混合的影响，研究发现：如果进气流速增大，将使气固系统活动更加剧烈，可以更快进入稳定流化状态，而且混合更均匀，同时发喷动气速增加，颗粒在喷射区和环形密相区的平均速度增加，浓度减小，平均碰撞频率减小；流化气速增加，颗粒速度、浓度和碰撞频率与喷动气速度增加具有相同的变化规律。刘安源等人对二维鼓泡流化床内瞬时和局部传热系数进行了模拟[37]，得到了瞬态表面传热系数随气速及局部传热系数随高度的变化规律。他们又通过建立流化床燃烧的数学模型[38]，模拟了流化床锅炉热烟气点火过程。1.5本文主要研究内容目前，流化系统中DEM方法的应用大多集中在流态动力学特性方面，对于热量传递，燃烧等过程研究较少，而由于流态化过程中热质传递及化学反应对流体动力特性有较大影响，因此可以预见，与一般稀疏气固两相流研究中颗粒轨道模型的发展过程一样，对流态化系统的DEM模拟也将会沿着气固流动、热量传递、质量变化的方向发展。本文将在在学习研究国内外现有研究成果的基础上，建立描述流化床内气固两相流动的DEM数学模型，进而对流化床的气固流动进行模拟，并研究分析模拟结果，主要工作内容为：介绍流化床内气固两相流动数学模型，根据数学模型设计出可用于模拟流化床气固流动过程的通用程序；利用所编写程序对流化床内的气固两相流动特性进行数值研究，包括对单双喷气孔两种进气方式的分别模拟及分析比较；在单双喷口研究的基础上，对多喷口流化床的气固流动进行模拟，分析喷口个数、喷口位置、进气速度、颗粒密度、弹性系数等对气固流动的影响。最后对本文工作做出总结，并对以后的研究提出展望。

第2章流化床气固两相流动DEM模拟数学模型由于DEM方法可以从单个颗粒的层次上去描述气固流动过程，对颗粒的形状、大小都可以按照要求设定，同时能够精确的知道每时每刻气固流动的具体情况并可以通过模型设置使颗粒产生物理或化学变化，对研究对象具有很强的适应性，虽然对计算机计算能力的要求比较高，但随着近年来计算机技术的发展，有越来越到的研究者利用计算机采用DEM方法研究气固流动系统的运动规律。本文也将采用离散单元法与计算流体力学相结合的方法来模拟研究流化床内的气固两相流体流动、燃烧及传热特性。2.1离散单元法原理DEM方法按照处理碰撞的不同方式可以分为硬球模型及软球模型两种，这两种方法在第一章已经做了详细的论述和比较，这里不再展开，基于气固流化系统中颗粒浓度较大，硬球模型无法处理，同时由于DEM软球模型在处理颗粒碰撞和传热过程上的方便，本文模拟将采用DEM软球模型。流化床DEM软球模型中，考虑碰撞两球之间的形变，在已知弹性系数的情况下，可以得到其所受碰撞力，再利用牛顿第二定律可以得出球的加速度，如果时间步长选取的足够小，算出颗粒每时每刻的位移，再由上一时刻的位置便可更新得到下一时刻的位置。在以上所述方法中，时间步长的选取必须要合适才能保证稳定的计算结果，由于在碰撞过程中，一个时间步长内颗粒所受作用力、当前速度和加速度都假定为是常数，而实际颗粒碰撞过程中，颗粒的运动参数随时在变化，同时为了节省计算时间，假定一个时间步长内颗粒只受其紧邻的颗粒作用力的影响，更远的颗粒则认为对其没有影响，所以时间步长的必须取的很小才能正确有效模拟出碰撞过程。2.1.1时间步长的确定对于时间步长的选取，通常认为步长取的越小，计算精度就越高，同时计算也越稳定，但由于受计算机性能及计算时间所限制，研究者通过实验来确定在保证精度与稳定的情况下，所能采用的最大时间步长：Cundall等人对物块通过弹簧连接壁面这一单自由度系统进行研究，分析得出了其最大临界时间步长[17]：(2-1)其中，m为的物块的质量，k为弹簧的弹性系数，认为如果所取时间步长大于此临界时间步长则计算是不稳定的。Tsuji等人通过试验发现上述条件并不能概括所有的计算情形，同时，他们还提出了一种改进的方法[18]。他们对作用力和位移之间的变化关系进行分析，得出了物块-弹簧-壁面系统的非线性运动方程：(2-2)他们通过对初始速度为1m/s的颗粒-弹簧-壁面单个颗粒碰撞系统进行试验，计算得到的临界时间步长为3×10-5s，对于多个颗粒组成的系统，通过计算得到最大临界时间步长为2×10-5s。本文为保证计算稳定性，选取时间步长为1×10-5s。2.1.2颗粒运动控制方程由动力学理论可知，颗粒的运动是由其整体的平移及自身的转动构成。(1)颗粒的平移，有：(2-3)其中和分别为颗粒的质量和速度，为该颗粒所受到的总作用力。(2)颗粒的转动，有：(2-4)其中和分别为颗粒i的转动惯量和角速度，对于球形颗粒，Ri为颗粒i的半径，为由于颗粒间碰撞而对颗粒i产生的扭转力矩的总和。在实际流化床中，颗粒受力极其复杂，包括自身重力、所受流体的作用力及颗粒之间的作用力，其中流体作用力又可分为流体曳力、压力梯度力、附加质量力、Basset力、Magnus力、Saffman力、浮力及热泳力等；颗粒之间相互作用力包括碰撞作用力、摩擦力、范德华力、静电力及液体桥力等。考虑模型复杂性及可行性，本文只考虑颗粒本身重力、流体曳力及颗粒间的碰撞和摩擦力，其它作用力相对较小，选择忽略不计。2.1.3颗粒碰撞力学模型颗粒间相互作用力包括弹性力及粘性阻尼力，法向力fn,ij分解为法向弹性力fcn,ij和法向粘性阻尼力fdn,ij，切向力ft,ij分解为切向弹性力fct,ij和切向粘性阻尼力fdt,ij。图2-1颗粒碰撞受力模型图2-1为颗粒碰撞软球模型，即假设颗粒在碰撞中发生了形变，产生了一个碰撞力，而碰撞力的大小跟形变量有关系。图中弹簧模拟了形变效应，阻尼器模拟了衰变效应，偶联器作用是在碰撞时连接相互接触的两个颗粒而在碰撞后又产生排斥力而使颗粒开始分离。对于两个颗粒组成的碰撞系统，其总作用力为：(2-5)在实际的流化床系统中，由于其高颗粒浓度的特性，每个颗粒都会与大量颗粒发生碰撞，此时，颗粒i所受的作用力为所有与之碰撞的颗粒对其作用力之和：(2-6)式中，ni为同一时间与颗粒i相碰撞的颗粒的个数。当颗粒的碰撞对象为壁面时，颗粒所受作用力的计算方法与两个颗粒碰撞相同，只需令第二个颗粒的半径、质量都为无限大。下面具体介绍颗粒间的法向力、切向力计算，颗粒的弹性碰撞采用线性作用力模型，当颗粒j与颗粒i碰撞时，颗粒间的法向弹性力与两颗粒在法向相对位移成正比，即：(2-7)式中，k为颗粒弹性系数，下标n表示法向方向，下标c代表弹性方向，为颗粒间相对位移，(2-8)式中，为颗粒i对颗粒j的相对速度矢量，，为颗粒i沿法线方向的单位矢量，方向有球心指向碰撞接触点；为时间步长。碰撞阻尼力由下式计算：(2-9)式中，为阻尼系数，与颗粒物性相关，下标d代表碰撞阻尼力，为颗粒碰撞法线相对速度，。颗粒碰撞总的法向力为颗粒法向阻尼作用力和法向弹性碰撞作用力之和：(2-10)颗粒所有切向力可分为有无滑移两种情况，当颗粒间不发生相对滑移时，与法向力的计算类似，颗粒切向作用力可以分为切向弹性力和切向阻尼力。作用在颗粒i上的切向弹性力为：(2-11)作用在颗粒i上的切向阻尼力为：(2-12)颗粒总切向力为：(2-13)当颗粒碰撞出现滑移时，由力学知识可知，此时切向的作用力超过颗粒的最大静摩擦力，此时切向作用力方向与无滑移时颗粒所受总切向力方向相同，其绝对值大小为：(2-14)式中为固体颗粒表面滑动摩擦系数。2.1.4碰撞对象的搜索算法由于在计算时如果对每个颗粒都对其他所有颗粒进行碰撞计算，一一判断是否发生碰撞，其计算量会非常大，因此如何建立合适的碰撞对象搜索算法，对减小计算量，节省计算时间非常重要。借鉴欧阳洁[39]及Hoomans[40]等人在硬球DEM模拟中使用的最近邻居表方法（thenearestneighborlist），本文在模拟时首先将计算区域分为若干小区域，考虑颗粒之间碰撞时，将颗粒视为等直径的球体单元，当要计算颗粒所受作用力时，先确定颗粒的位置，然后只需搜索该颗粒所在网格及临近的网格内其他颗粒即可。2.2流体相数学模型及求解方法对流体相需要采用Navier-Stokers方程求解，而该方程的适用条件需要计算区域是连续的，但本文研究对象为具有强烈运动、混合、相互作用的流化床气固两相流，很难对这一对象进行非常精确的模拟[41]。为解决这一问题，Anderson和Jackson提出了局部平均的方法[42,43]，该方法把连续计算区域划分为多个网格，每个网格包含多个固体颗粒，在网格间速度变化远小于整个系统速度宏观变化这一前提下，提出用每个离散网格内的局部平均变量代替连续相流体的点变量。2.2.1流体相数学模型采用局部平均法时，流体相连续性方程及动量方程如下：连续性方程：(2-15)式中，为颗粒质量源项，如果不考虑燃烧则＝0。动量方程：(2-16)(2-17)式中：(2-18)(2-19)式中，为湍流粘度，，为相互作用力在x,y方向上的分量。我们通常把方程中各个项分别定义为非稳态项，对流项，扩散项，源项，具体见下表：表2-1动量方程中各项定义名称相非稳态项，对流项，，，扩散项,,,源项气相湍流模拟采用方程模型描述，控制方程如下：(2-20)(2-21)式中：(2-22)(2-23)(2-24)带入得：(2-25)(2-26)同动量方程，方程各项定义如下表：表2-2k方程中各项定义名称项非稳态项对流项扩散项产生项耗散项表2-3方程中各项定义名称项非稳态项对流项扩散项产生项耗散项其中，，l是湍流脉动的特征长度。2.2.2流体相控制方程组的求解流体相控制方程可改写成对流—扩散方程的通用形式：(2-27)上式中，表示不同的变量，如，等，对于连续性方程，取1，为变量的交换系数，为源项。利用解析法求解流体相控制方程组，会因为其不具有的线性特性而且存在着与颗粒相运动方程之间强烈的耦合作用，求解会遇到很大障碍，由此我们一般通过数值方法求解。偏微分方程组数值求解一般包括以下几个过程：控制方程的离散；计算区域的离散；线性代数方程组的求解；控制方程组的求解算法。控制方程的离散现在常用的微分方程离散方法主要有：(1)有限元法是把计算区域分成一组离散的单元，然后通过积分控制方程得出离散方程；(2)边界元法是通过格林函数及适当的权函数用边界上的积分方程代替求解域上的偏微分方程；(3)有限分析法是利用求解线性化了的微分方程，把相近积分网格点上的值关联起来，来实现方程的离散；(4)有限差分法是将网格节点的集合代替求解区域，对这些节点用泰勒公式或控制容积法去离散方程；(5)有限容积法是对流动与传热守恒型方程在控制容积上进行积分，从而获得控制方程的离散形式。本文采用有限差分法来离散流体相控制方程组，式(2-27)可改写为：(2-28)(2-29)(2-30)式中J表示包括对流项和扩散项作用的总通量。如图2-2所示控制容积对对流扩散通用方程进行离散。图2-2控制容积示意图方程离散结果如下：(2-31)其中：(2-32)(2-33)(2-34)(2-35)(2-36)(2-37)(2-38)(2-39)其中，源项，F为网格界面的流量；，代表网格界面处对流项与扩散作用的对比。计算区域的离散常用的网格划分方法有内节点和外节点两种[20]。本文采取内节点法，即先确定控制容积再确定节点的方法来划分计算网格，如图2-3。图2-3计算网格划分示意图为了解决由于一阶导数项的存在而使在一般的网格上检测不出不正确压力场的存在的问题，人们提出了交错网格的方法[20]，即把速度与压力分别存储在不同的网格系统中，如图2-4所示，压力控制容积为主控制容积，的控制容积是以主控制容积的东西界面为中心，的控制容积是以主控制容积的南北界面为中心。图2-4交错网格示意图离散方程组的求解流体相动量方程离散后，计算区域内所有内节点均有式(2-31)的表示形式，此形式的线性代数方程组可以采用三对角阵算法(Tridiagonalmatrixmethod，TDMA)结合逐行迭代法来求解[44]。流体相动量方程的求解采用在计算流体力学中广泛应用的SIMPLE(Semi-ImplicitMethodforPressureLinkedEquation)方法[44]，即求解压力耦合方程的半隐方法。由于在求解速度场时，会发生相邻两次迭代之间的速度变化过快而导致迭代发散的情况，本文在进行SIMPLE算法的过程中，采用亚松弛技术(under-relaxation)，用通用变量表示如下：(2-39)式中为上次迭代所求解，为松弛因子，将上式改写后可得：(2-40)2.3气固两相相互作用力分析由于对固相采用的是DEM方法求解每一个离散颗粒，而对流体相是采用CFD方法求解控制体，所以当进行气固两相流的模拟时，气固两相之间的耦合尤为重要，本文利用牛顿第三定律实现气固两相间作用力的耦合，固体受气体的曳力，而气体也受到固体的反作用力，体现为作用在气体控制体上的体积作用力。2.3.1气体对颗粒的作用力本文模型中，气体对颗粒的作用力在单颗粒水平上确定，气体对颗粒曳力与两相间的相对速度及物性有关，其具体表达式为:(2-41)由曳力系数选取的不同，目前有以下几种不同的曳力模型：Syamlal-O’Brien模型:(2-42)(2-43)XuandYu模型[25]：(2-44)(2-45)(2-46)Gidaspow模型[45]：当空隙率较低，≤0.8时，利用Ergun等式得：(2-47)当空隙率较大，＞0.8时，采用Wen&Yu的修正方程得：(2-48）本文采用Gidaspow的曳力模型，则颗粒在流体中所受的曳力为：(2-49)(2-50)式中为颗粒直径，为流体密度，为流体速度，为颗粒速度，为空隙率，为流体剪切粘度，为颗粒体积，为有效曳力系数。(2-51)其中n为空隙率修正因子，一般取值为-4.65，为曳力系数，计算公式如下：当时，(2-52)当时，(2-53)为颗粒雷诺数，计算式为：(2-54)空隙率对气体及颗粒之间作用力的计算有很大影响，本文计算选用区域为二维平面，但为计算及描述方便，假定颗粒为球形，且从流化床实际运行角度出发，三维更接近与实际情况。本文采用Xu&Yu提出的准三维概念[25]，按三维的方法确定空隙率：(2-55)式中：为颗粒体积；为划分区域的控制体积；为所划分区域中的颗粒的数目，数值模拟计算中，空间的厚度设置为等于颗粒的直径D。2.3.2颗粒对流体的反作用力按牛顿第三定律，颗粒受到流体曳力，相对的流体也会受到颗粒对其的反作用力，这个力与曳力大小相等，方向相反。颗粒对流体的体积作用力为：(2-56)其中为气体对颗粒的作用力。2.4本章小结本章DEM数学模型及其求解方法做了具体介绍，同时也介绍了本文在进行数值模拟时对方程离散、求解的做法主要内容包括：1、对离散固体颗粒相采用属于Lagrangian范畴的DEM法来模拟，对连续流体相采用属于Euler方法范畴的Navier-Stokes方程模拟；2、采用内节点法进行求解区域的计算，对速度变量采用交错网格法；采用逐行迭代的TDMA方法求解动量离散化后得到的线性代数方程组；3、采用Di.Felice提出的经验关联式计算流体对颗粒的曳力，床层空隙率采用准三维方法确定；固体颗粒对流体相的反作用力采用牛顿第三定律计算。

第3章循环流化床气固流动特性的数值模拟本章将采用DEM方法对流化床内的气泡特性进行模拟研究。本章首先选对床内颗粒自由堆积过程进行模拟。进而在不同入口风速下模拟得到了单孔射流流化床内颗粒运动情况。同时在单孔射流基础上本章模拟还得到了双孔射流流化床内气固流动的运动情况，最后对单双孔流化床内颗粒流动情况做出对比分析。3.1自由堆积过程模拟本文采用重力沉降法确定颗粒的初始堆积位置，具体为：在零时刻利用随即函数分配颗粒初始的位置，即假定颗粒初始无规律散布整个空间，此时设定气流流速为零，颗粒只在重力的作用下进行自由落体运动，考虑颗粒相互间碰撞力，颗粒运动至床底部自由堆积，记录此时的颗粒位置，作为颗粒初始位置。下图为模拟计算网格示意图：图3-1堆积状态模拟计算网格示意图对图中所示计算区域，边界条件进行下述假定：(1)给定入射气流流速，壁面上采用无滑移边界条件，出口处采用自由滑移边界条件；(2)假定出口法向压力为零；(3)颗粒和壁面发生碰撞时，假定壁面为直径和质量无穷大的一个特殊颗粒，并取与一般颗粒相同的刚度系数，然后再利用前面的颗粒—颗粒碰撞计算模型进行计算。对于相关初始参数设定如下表：表3-1自由落体过程参数表参数名称参数取值颗粒密度/(kg/m3)2700弹性系数/(N/m)800摩擦系数0.3粘性系数恢复系数0.9颗粒直径(m)0.004颗粒数目2400床高(m)0.9床宽(m)0.15时间步长/(s)0.00001重力沉降过程模拟图：0s0.2s0.3s0.5s1.0s图3-2颗粒自由落体过程图图3-3颗粒最终沉积状态细节图3.2单孔射流流化床流动过程模拟在颗粒自由落体模拟的基础上加入气相，使气流以一定流速从流化床底部中间位置射入床内，模拟各参数见下表：表3-2单孔射流模拟参数表参数名称参数取值颗粒密度/(kg/m3)2700弹性系数/(N/m)800摩擦系数0.3粘性系数恢复系数0.9床宽(m)0.15床高(m)0.9进气孔宽度(m)0.1气体粘度/(kg/(ms))1.810-5气体密度/(kg/m3)1.205时间步长/(s)0.00001进气孔速度设为35m/s，模拟流动过程时间为2秒，单孔射流流动过程模拟结果及颗粒，气体垂直、水平速度分布见下图：0s0.3s0.6s0.9s图3-4(a)单孔射流流化过程模拟1.2s1.5s1.8s2.0s图3-4(b)单孔射流流化过程模拟图3-5(a)单口颗粒垂直速度分布图3-5(b)单口颗粒水平速度分布图3-5(c)单口气体垂直速度分布图3-5(d)单口气体水平速度分布由图3-4中可以看出气体从布风板进入床层后，进气口上方颗粒受到气体流动作用力而运动。在局部，颗粒群内部形成气泡，在整体，颗粒群上方形成凸形曲面，随着时间的推移，进入床层的气流使气泡不断形成，成长，移动，最终在床层表面由于局部外部压力减小而破裂，同时颗粒随气泡的上升而运动，从底部被携带至床层上方再由两侧落下，同时两侧颗粒由于卷吸作用到达中部，如此循环往复，最终达到稳定流化状态。而图3-5则描述了流化1s时颗粒、气体分别在水平、垂直方向上的速度分布，在水平方向上，由于卷吸作用，两者都是以喷口为中心向两侧流动；在垂直方向上，由于进气孔处气体的速度较大同时带动了颗粒有较大的速度分布。为验证本文模拟正确性，下取Xu&Yu等人对单孔射流模拟结果[25]进行对比:图3-6(a)Xu&Yu等人对单孔射流模拟结果图3-6(b)Xu&Yu等人对单孔射流模拟结果与文献[25]模拟结果(图3-6)相比，本文的模拟结果(图3-4)具有很大的相似度，可以证明本文模拟结果具有一定的可信度。由图3-4可见，在气流刚进入床内时，床层出现松动，颗粒之间产生较小气泡，气泡往上流动，同时带动周围颗粒一起向上运动，随时间的继续气泡不断增大，同时伴有气泡的破裂，气泡间的挤压、吞并等现象发生，进而带动更多颗粒一起运动，到稳定流化阶段，气体散布整个床层，颗粒团完全流化，颗粒进入上升-下降-上升循环往复的运动阶段。3.3双孔射流流化床流动过程模拟现有研究普遍都是进行单喷口流化床的研究，与实际现场中的多喷口工况有一定偏差，本文引入双喷口模型，在相同速度下与单喷口模型模拟结果进行对比，并对双喷口给出不同入口速度，分析研究相邻气泡之间相互作用，进而探讨流化床内气固流动机理特性。为更好分析模拟结果，减少因进气孔位置不同而引起的实验偏差，双孔射流模拟时将流化床底部四等分，进气孔设置在第1、3等分点处。对称布置与流化床底部参数设置如下表：表3-3双孔射流模拟参数表参数名称参数取值颗粒密度/(kg/m3)2700弹性系数/(N/m)800摩擦系数0.3粘性系数恢复系数0.9床宽(m)0.15床高(m)0.9进气孔宽度(m)0.1进气孔个数2气体粘度/(kg/(ms))1.810-5气体密度/(kg/m3)1.205时间步长/(s)0.00001本文模拟分别取双口相同射流速度及不同射流速度两种工况模拟，其中相同射流速度时取射流速度与单孔射流模拟相同，为35m/s。模拟结果见图3-7：0s0.3s0.6s0.9s图3-7(a)双孔射流相同流速流化过程模拟1.2s1.5s1.8s2.0s图3-7(b)双孔射流相同流速流化过程模拟图3-8(a)双口颗粒垂直速度分布图3-8(b)双口颗粒水平速度分布图3-8(c)双口气体水平速度分布图3-8(d)双口气体垂直速度分布如图3-7所示，本文模拟了双喷口流化床内气固两相流动的运动规律，由于两喷口所取气体入口速度相同，可以清楚看出气体从进气口喷出，在进气口形成两个气泡，进而两气泡成长，上升，在0.9s时两气泡相互接触使得局部压力变小而破裂、融合，进而推动整个颗粒群的运动，颗粒在两股气流共同作用下运动，最终达到稳定流化状态。而在图3-8中可以看出水平方向上，与单喷口流化床相比，颗粒与气体的速度存在着强烈的非均匀性，在垂直方向上，速度仍然是随高度的上升而降低。双孔射流速度不同时模拟结果见图3-9，其中左进气孔入口射流速度为30m/s，右进气孔入口射流速度为25m/s，其他参数与双孔相同射流速度保持不变。0s0.3s0.6s0.9s图3-9(a)双孔射流不同流速流化过程模拟1.2s1.5s1.8s2.0s图3-9(b)双孔射流不同流速流化过程模拟如图中所示，除与图3-7中具有共同的流动特征外，在0.3s可以很清楚的看出左侧气泡体积大于右侧气泡，随时间推移，颗粒群在左侧的高度也明显高于右侧，在1.2s左右，颗粒整体达到最大高度，而后下落一定距离后达到最终稳定状态。可以看出，进气孔处流速对床层高度有直接的影响关系，流速越大床层高度越高，虽然颗粒有强烈的混合、扩散，但这种高度差基本保持不变。3.4本章小结在模拟颗粒自由落体堆积过程的基础上，本章首先对单孔射流流化床内的气固流动特性进行了模拟研究，对比单孔射流，又在相同、不同射流速度下对双孔射流流化床内的流动特性进行了模拟，综合上述单喷口、双喷口同速及双喷口不同速度的流化模拟结果，我们可以发现：1、颗粒的大范围位移都发生在气泡附近，颗粒随气体的运动而运动，气体运动是颗粒运动的主要动力源；2、双喷口与单喷口相比，颗粒群整体所能达到的最高高度与最终稳定状态都明显较高，但由于颗粒被吹穿，气体穿过颗粒群，气体作用力并没有完全作用于颗粒，所以速度为35m/s的双喷口高度并非为相同速度下单喷口高度的两倍；3、对比单、双喷口流化系统中的气泡形态，可以发现，双喷口流化系统中，气泡的形状会因其他气泡的接近而改变，同时相邻气泡之间存在横向和纵向双方向的的聚并，挤压，破裂等现象，由此，床层内颗粒的运动和扩散作用也更加强烈。

第4章多喷口流化床气固流动特性模拟在实际生产中，循环流化床锅炉往往分布多个进气口，为分析流化床中多喷口间，本章将在单、双喷口流化床研究的基础上，继续分析多喷口流化床的流动特性，详细分析喷口个数、喷口位置、进气速度、颗粒密度、弹性系数等对气固流动的影响，对比单双喷口系统，研究多喷口系统对颗粒混合，气固运动剧烈程度的影响。4.1三孔射流模拟对象及参数设定三孔射流模拟对象见图4-1：图4-1三孔射流计算网格示意图如图中所示进气孔分布在第5、第9以及第13网格处，模拟所建数学模型与单双孔相似，模拟参数设定见表4-1：表4-1三孔射流模拟参数表参数名称参数取值颗粒密度/(kg/m3)2700摩擦系数0.3粘性系数恢复系数0.9床宽(m)0.15床高(m)0.9进气孔宽度(m)0.1进气孔个数3气体粘度/(kg/(ms))1.810-5气体密度/(kg/m3)1.205时间步长/(s)0.000014.1.1弹性系数对流动情况的影响模拟进行三孔射流模拟过程中，开始选用单双孔模拟时所用弹性系数，发现并不能得到理想的颗粒流动情况，射流速度设为30m/s，模拟情况如图4-2所示：0.2s0.4s0.6s0.8s1.0s图4-2(a)弹性系数800N/m时三喷孔流化过程模拟1.2s1.4s1.6s1.8s2.0s图4-2(b)弹性系数800N/m时三喷孔流化过程模拟由图4-2中，可以看到在流化初始阶段，床层下方颗粒整体被托起，同时进入充分流化阶段后，床层下部有大量颗粒致密排列，并没有发生流化现象，这一现象被称为“互锁”现象，为分析互锁现象的影响因素，在其他条件不变的情况下，改变弹性系数，分别达到了1000N/m，500N/m，300N/m下，初始及流化充分阶段气固流动图：图4-3(a)0.5s时弹性系数分别为1000N/m，500N/m，300N/m时颗粒流动细节图图4-3(b)1.5s时弹性系数分别为1000N/m，500N/m，300N/m时颗粒流动图如图中所示在流化初期(0.5秒时)，当弹性系数取的较大时，床层并没有被射流气体击穿，而是整体被托起，模拟所取弹性系数越大，床层下部空出区域越大；在完全流化后(1.5秒时)，床层下部也没有出现应有的三孔气流击穿床层的情形，而床层上方颗粒则有较为充分的混合现象发生，这主要是由于当选用的刚度越大时，床层底部颗粒存在的“互锁”现象越明显，弹性系数越大床层下部的不动区厚度就越大，气流通过颗粒间空隙向上运动，上部颗粒由于压力较小，所以被流化，导致形成下部致密，上部分散流化情况。当弹性系数设为180N/m，速度设为28m/s时，三孔射流流动过程模拟结果及颗粒，气体垂直、水平速度分布如图4-4、图4-5所示：0.2s0.4s0.6s0.8s1.0s图4-4(a)三喷孔流化过程模拟如图4-4所示，当弹性系数取180N/m时，入射气体可突破底层颗粒“互锁”作用，气体从三个进气孔进入床内，生成三个气泡，气泡生长并向上运动，同时床层膨胀，在1.2s左右，三个气泡发生吞并，融合现象，1.5s后整个床层进入完全流化状态，整个床层不再增高，同时，颗粒充分混合，三个进气孔无明显的单个大气泡生成，进气孔处与其他底部壁面无明显界限。1.2s1.4s1.6s1.8s2.0s图4-4(b)三喷孔流化过程模拟图4-5(a)三喷孔颗粒垂直速度分布图4-5(b)三喷孔颗粒水平速度分布图4-5(c)三喷孔气体垂直速度分布图4-5(d)三喷孔气体水平速度分布图4-5为流化1s时三喷孔颗粒、气体位置为床高0.05m、0.10m及0.20m处的的垂直、水平速度分布，从图中可以看出在靠近进气孔处(0.05m)，颗粒的垂直速度受到气流的影响很大，随着床层高度的增加同一高度处颗粒速度差逐渐缩小，在床高0.2m处，颗粒速度基本一致，气体垂直速度也是随高度的增加，由于推动颗粒运动，动能减小，速度逐渐减小同时趋于同一值。在图4-5(b)、图4-5(d)中我们可以看到，在底部受气体影响较大，喷口处颗粒水平速度较小，随床层高度的增高，颗粒流化程度也越深，颗粒在水平方向上运动活跃。4.1.2不对称喷口布置对流动过程影响模拟为分析喷口布置位置不同对气固流动的影响，本文对三孔射流进气孔位置进行设定，其他参数不变，首先设置三个进气孔分别为底部布风板第3，7，15处（记为左侧布置）。速度设为28m/s，其气固流动图及颗粒、气体垂直、水平速度见图4-6、图4-7：0.2s0.4s0.6s0.8s1.0s图4-6(a)进气孔左侧布置时流动过程1.2s1.4s1.6s1.8s2.0s图4-6(b)进气孔左侧布置时流动过程图4-7(a)进气孔左侧布置颗粒垂直速度分布图4-7(a)进气孔左侧布置颗粒水平速度分布图4-7(c)进气孔左侧布置气体垂直速度分布图4-7(d)进气孔左侧布置气体水平速度分布当三个个进气孔分别为底部布风板第3，11，15处（记为右侧布置），速度同样设为28m/s，其气固流动图及颗粒、气体垂直、水平速度见图4-8、图4-9：0.2s0.4s0.6s0.8s1.0s图4-8(a)进气孔右侧布置时流动过程1.2s1.4s1.6s1.8s2.0s图4-8(b)进气孔右侧布置时流动过程图4-9(a)进气孔右侧布置颗粒垂直速度分布图4-9(a)进气孔右侧布置颗粒水平速度分布图4-9(c)进气孔右侧布置气体垂直速度分布图4-9(d)进气孔右侧布置气体水平速度分布综合图4-5到图4-8，可以看出喷口位置对气固流动有较大的影响，对比图4-5，图4-7所示流动过程，由于中间喷气孔位置的不同，所得模拟结果正好相反，图4-5中流化开始，床层由左侧突起，左侧颗粒垂直速度也较大，两个喷气孔有效的推动了左侧颗粒的流化过程，同时进入充分流化阶段后，左侧形成气泡区也明显大于右侧，图4-7中右侧先被气流吹起，进入充分流化阶段后，右侧气泡也明显大于左侧，同时，两图中气体垂直速度与颗粒速度一致，在水平方向上，喷气孔的布置对速度分布也有明显的影响，喷气孔密集处，颗粒与气体的水平速度都较大，综合垂直速度，可以认为，喷气孔密集程度对气固垂直、水平速度都有较大的影响，速度随密集程度的增大而增大，气固流动程度随密集程度的增大而加深。4.2四孔射流流化床模拟4.2.1模拟对象及参数设定四孔射流模拟计算区域见图4-10：图4-10四孔射流计算网格示意图如图中所示进气孔分布在第3、第7、第11以及第15网格处，模拟所建数学模型与单双孔相似，模拟参数设定见表4-2：表4-2四孔射流模拟参数表参数名称参数取值颗粒密度/(kg/m3)2700摩擦系数0.3粘性系数恢复系数0.9床宽(m)0.15床高(m)0.9续表4-2进气孔宽度(m)0.1进气孔个数4气体粘度/(kg/(ms))1.810-5气体密度/(kg/m3)1.205时间步长/(s)0.000014.2.2四孔射流模拟结果与三孔模拟相同，四孔模拟先考虑弹性系数的选取，当选取弹性系数与三孔模拟一样为800N/m时，模拟结果见图4-11：0.2s0.4s0.6s0.8s1.0s图4-11(a)弹性系数800N/m时四喷孔流化过程模拟1.2s1.4s1.6s1.8s2.0s图4-11(b)弹性系数800N/m时四喷孔流化过程模拟如图4-11所示，在弹性系数相同的情况下，四孔模拟在流化初期颗粒被较高托起，在流化充分阶段，下部也有小厚的颗粒层没有被流化，这主要是由于四孔均匀布置在底部布风板，产生曳力均匀作用于底部颗粒，颗粒团更容易被托起而不是击穿，同时在充分流化阶段，气体也会渗透颗粒团而使上部的颗粒先流化起来。经过调整模拟试验，当弹性系数取100N/m时所得气固流动图较为合理，四孔射流气固流动图及颗粒、气体垂直、水平速度见图4-12、图4-13：0.2s0.4s0.6s0.8s1.0s图4-12(a)四喷孔流化过程模拟1.2s1.4s1.6s1.8s2.0s图4-12(b)四喷孔流化过程模拟图4-13(a)四喷孔颗粒垂直速度分布图4-13(b)四喷孔颗粒水平速度分布图4-13(c)四喷孔气体垂直速度分布图4-13(d)四喷孔气体水平速度分布在图4-12中，我们可以清楚看到，与三孔模拟相似，在流化初期，四个进气孔之间相互影响不大，分别生成四个气泡，不同的是，气泡在相互影响吞并、融合时的尺寸明显小于三孔模拟时气泡破裂的尺寸，同时四孔射流进入稳定流化阶段的时间也比三孔射流稍短，在图4-13中可以看到，颗粒垂直、水平速度分布与前面几种模拟所得规律一致，同时颗粒在水平方向上的速度也较大，颗粒运动较为剧烈。4.3本章小结综合本章多喷口流化系统模拟，可以看出随着流化进气孔数的增加，弹性系数的选取对模拟结果有较大影响，孔数越多，床内颗粒所受气体曳力也就较平均，在弹性系数选取较大时，颗粒之间的“互锁”现象就越明显，同时随着进气孔数的增大，颗粒与气体的混合作用也就越强烈，流化床进入充分流化阶段的时间也越短，在对比不同喷口位置的模拟结果后，认为喷气孔密集程度对气固垂直、水平速度都有较大的影响，速度随密集程度的增大而增大，气固流动程度随密集程度的增大而加深。

第5章全文总结及展望鉴于流化床，气固两相流动越来越广泛的应用，随着科技的进步，研究手段的多元化，人们肯定会对流化床的流动、燃烧及传热特性有更深入更细致的研究，而能够在颗粒层面上反映流化床运动、燃烧、传热机理的DEM方法也一定会受到广大研究者的重视。5.1全文工作总结本文在学习研究国内外现有研究成果的基础上，建立了描述流化床内气固两相流动DEM数学模型并分别对单孔、双孔、三孔以及四孔流化床的气固流动做出了模拟，通过对比分析模拟结果得到了下列研究成果：(1)通过模拟得到了单孔射流流化床内的气固流动过程，分析流化床内气固流动的起因与具体流动过程。(2)对双孔射流流化床进行模拟得到了床内床内的气固流动过程，并对比分析单孔射流时结果，发现双喷口流化系统中，气泡的形状会因其他气泡的接近而改变，同时相邻气泡之间存在横向和纵向双方向的的聚并，挤压，破裂等现象，由此，床层内颗粒的运动和扩散作用也更加强烈。(3)对多喷口射流流化床系统进行模拟，通过分析喷口个数、喷口位置、进气速度、颗粒密度、弹性系数等对气固流动的影响，对比单双喷口系统，认为进气孔数的增加时，弹性系数是颗粒之间的“互锁”现象的主要影响因子；随着进气孔数的增大，颗粒与气体的混合作用也就越强烈，流化床进入充分流化阶段的时间也越短；同时，喷气孔密集程度对气固垂直、水平速度也有较大的影响，速度随密集程度的增大而增大，气固流动程度随密集程度的增大而加深。5.2下一步工作展望(1)将文中所做二维流化床研究扩展到三维情况下，建立三维流化床数值模拟模型，使研究结论更接近实际；(2)进一步完善流化床模型的精确性，如考虑旋转、颗粒大小不同，考虑颗粒燃烧直径的变化及破碎等因素，考虑颗粒形状不规则等工况，提高模型的适用性；(3)希望能够探索出一些能够加快计算速度的研究算法，如建立流化床流动、燃烧与传热的并行计算系统等。

参考文献[1]环境保护部,国家统计局,农业部.第一次全国污染源普查公报[R].2010.02-10[2]刘亚妮.循环流化床锅炉数学模型及数值模拟研究[D].武汉大学硕士学位论文,2005[3]R.D.Toomey,H.F.Johnstone.Gaseousfluidizationofsolidparticles[J].Chem.Engrg.Prog,1952,48(5):20-226[4]J.F.Davidson,D.Harrison.Fluidizedparticles[M].NewYork:CambridgeUniversity[5]D.Kunni,O.Levenspiel.FluidizationEngineering[M].NewYorkl：Wiley,1969[6]李静海,郭慕孙,L.Reh.循环流化床能量最小多尺度作用模型[J].中国科学(B辑),1992,35(11):1127-1136[7]S.L.Soo.Fluiddynamicsofmultiphasesystems[M].MA