专利名称：一种测量颗粒流体两相流反应器内流动参数分布的方法
技术领域：
本发明涉及颗粒流体两相流反应器，特别涉及一种测量颗粒流体两相流反应器内流动参数分布的方法。
背景技术：
颗粒流体两相流广泛存在于过程工业系统中。随物性和操作参数的变化，两相流动呈现出复杂的结构和参数分布变化，这对于工业反应器内的压降分布、传质、传热效率及反应转化率有直接的，有时甚至是数量级的影响。因此，如何准确描述颗粒流体两相流反应器内的流动结构和参数分布对过程工业的设计、运行、控制与放大是一项重要工作。两相流的相间曳力预测方法则是实现这一工作目标的核心内容之一。
目前两相流中常用的相间曳力系数源自固定床压降外推关系或者散式流态化条件下的实验数据关联形成的有效曳力系数，其本质上是对单颗粒曳力系数引入了假想均匀分布的颗粒群的平均扰动影响。这种计算方式被绝大多数研究者和实践者采纳，并被国际上通行的计算流体力学(CFD)软件CFX，Fluent，Phoenics，STAR-CD等推荐(参见参考文献1Gidaspow，D.，Jung，J.，Singh，R.K.，Powder Technology，Vol.148，P.123(2004))。这种方法适用于流动变化相对均匀、缓慢的体系。对于复杂的快速流动，如石油工业中常见的流化催化裂化(FCC)反应器，电力行业中的循环流化床锅炉(CFBB)等核心单元内的流动，由于颗粒成团结构大小不定，微观和介观尺度界限模糊，无法找到满足连续介质统计要求而且内部均匀无结构的微元空间，因此结构对计算的影响不可避免。在此情况下，不考虑结构影响的计算值与实际值相比，差异很大，结果处于不同的数量级，其直接后果是设计或者仿真计算时的流场颗粒浓度及反应器内压降分布失真。必须在两相流计算过程中引入结构的复杂影响，才能准确预测两相流反应器内的流动分布，并真正掌握反应器内部规律，从根本上保障设计、运行、控制与放大。
两相流复杂结构的变化行为与曳力相互耦合，难以实验测量和模型预测，主要体现在●非均匀团聚结构变化迅速，无确定几何外形，其几何与力学参数难以表征，很难直接测量非均匀动态结构下的曳力变化；●曳力系数所代表的相间交换强度与多个尺度下的物性、操作参数相关它不仅与反应器平均的两相表观速度、浓度有关，与局部速度和结构的变化也密切相关。这种跨尺度的关联很难用先缩小研究区域，再从微观局部变化延伸至宏观的方法得到。
前人工作注重于区别团聚物颗粒和单个颗粒的受力差异，介观团聚结构采用简化的，而且往往是静态的，形状不变的大颗粒代替。已有的少量探索性工作采用实验关联或经验的方法获得团聚物的特征尺寸，将颗粒团聚物的动态结构表述为一种静态大颗粒和原始颗粒的混合物，并以此为基础计算相间曳力(参见参考文献2Heynderickx，G.J.，Das，A.K.，De Wilde，J.，Marin，G.B.，Industrial & EngineeringChemistry Research，Vol.43，P.4635，(2004))。这种静态、经验性的估计方法，本质上无法从原有实验条件外推应用至工业反应器中。
已有工作只注意了亚结构微观尺度对介观结构的影响(参见参考文献3Agrawal，K.，Loezos，P.N.，Syamlal，M.，Sundaresan，S.Journal of Fluid Mechanics.Vol.445，P.151.(2001))，试图从样本微元空间的计算结果直接关联得到复杂的介观结构对曳力的影响，而没有考虑跨尺度关联对介观结构的影响。实际过程的样本空间很难选取，在工业反应器的应用中是不可行的。
发明内容
本发明的目的是为了克服现有技术没有考虑跨尺度关联对介观结构的影响，流动参数分布的计算结果与实际情况相差较远的缺点，从而提供一种可精确测量颗粒两相流反应器内的流动参数分布的方法。
为了实现上述目的，本发明提供了一种测量颗粒流体两相流反应器内流动参数分布的方法，包括1)、根据实际的物性和反应器的操作条件情况，初始化流场及边界条件；2)、利用双流体模型计算代码或商业化计算流体力学软件，计算每个空间微元内的质量守恒方程、动量守恒方程，得到速度、浓度分布(ug，up，εg)，所述的ug代表微元内的气体真实速度，所述的up代表微元内的固体真实速度，所述的εg代表微元内平均空隙率；3)、判断当前的速度和浓度分布是否满足全场的收敛标准，若满足全场的收敛标准，输出当前的速度和浓度分布，并结束操作，否则，执行下一步；4)、根据已知的反应器整体操作条件(Ug，Up)，计算并得到满足全局动力学的非稳态、非线性方程组的所有变量根的组合；所述的Ug代表气相表观速度；所述的Up代表颗粒相表观速度；5)、从所有变量根的组合中，寻找满足极值条件的最优根，并保存其中的介观结构参数εc，εf，dcl，所述的εc表示团聚物密相空隙率，所述的εf表示稀相空隙率，所述的dcl表示团聚物直径；6)、由步骤2)得到的速度、浓度分布(ug，up，εg)和步骤5)得到的介观结构参数εc，εf，dcl，在反应器的每个微元空间内求解满足微观局部动力学非稳态、非线性方程组的所有变量根的组合；所述的变量根包括Uf，Upf，Uc，Upc，f，ai，ac，af，其中Uf代表稀相流体速度，Upf代表稀相颗粒速度，Uc代表团聚物密相流体速度，Upc代表团聚物密相颗粒速度，f代表密相体积分数，ai代表介观界面加速度，ac代表密相加速度，af代表稀相加速度；7)、由步骤6)得到的变量根的组合，从中寻找满足极值条件的最优根；8)、由步骤5)和步骤7)所得到的微元内所有变量值(Uf，Upf，Uc，Upc，εc，εf，f，dcl，ai，ac，af)的组合，计算曳力和曳力系数；9)、将步骤8)得到的曳力系数代入步骤2)中所采用的双流体模型计算代码或商业化计算流体力学软件中，替换原有的曳力系数，然后转到步骤2)，做迭代计算。
上述技术方案中，所述的非稳态、非线性方程组包括单位体积微元空间内，浓相团聚物颗粒曳力Fdc平衡方程Fd,c=34CDc1-ϵcdpρg[Uc-ϵcUpc/(1-ϵc)]2=(1-ϵg)(ρp-ρg)(ac+g)]]>单位体积微元空间内，稀相颗粒受曳力Fdf平衡方程Fdf=34CDf1-ϵfdpρg[Uf-ϵfUpf/(1-ϵf)]2=(1-ϵf)(ρp-ρg)(af+g)]]>单位体积微元空间内介观相界面的曳力Fdi平衡方程Fd,i=34CDifdclρg[Uf-ϵfUpc/(1-ϵc)]2(1-f)2=f(ϵg-ϵc)(ρp-ρg)(ai+g)]]>
流体和颗粒相的质量守恒方程Ug＝fUc+(1-f)UfUp＝fUpc+(1-f)Upfεg＝fεc+(1-f)εf团聚物直径dcl公式dcldp=(Nst)max-(Nst)mfNst-(Nst)mf]]>其中，ρ，密度；Ug，气相表观速度；Up，颗粒相表观速度；Uf，稀相流体速度；Upf，稀相颗粒速度；Uc，团聚物密相流体速度；Upc，团聚物密相颗粒速度；εg，微元内平均空隙率；εc，团聚物密相空隙率；εf，稀相空隙率；f，密相体积分数；dcl，团聚物直径；dp，单颗粒直径；ai，介观界面加速度；ac，密相加速度；af，稀相加速度；g，重力加速度；下标g，气相；下标p，固相；下标c，浓相；下标f，稀相。CDc，密相平均单颗粒曳力系数；CDf，稀相平均单颗粒曳力系数；Cdi，介观相间平均单个团聚物曳力系数；上述参数的定义式如下CD,f=(24+3.6Ref0.68)/(Ref·ϵf4.7),]]>Ref＝dpρgUsf/μgCD,c=(24+3.6Rec0.687)/(Rec·ϵc4.7),]]>Rec＝dpρgUsc/μgCD,i=(24+3.6Rei0.687)/[Rei·(1-f)4.7],]]>Rei＝dclρgUsi/μg各表观滑移速度定义为Usc＝Uc-εcUpc/(1-εc)Usf＝Uf-εfUpf/(1-εf)Usi＝Uf-εfUpc/(1-εc)(1-f)上述技术方案中，所述的步骤5)和步骤7)中，所述的极值条件是指单位质量的颗粒悬浮输送能量趋于最小，所述的单位质量的颗粒悬浮输送能量用Nst表示，其定义为Nst=ρp-ρgρp[Ug-ϵf-ϵg1-ϵgf(1-f)Uf].]]>
上述技术方案中，所述的步骤8)中，所述的曳力用Fd表示，其计算公式为Fd＝fFd，c+(1-f)Fd，f+Fd，i；所述的曳力系数用β表示，可由曳力得到，所述的曳力系数的计算公式为β=Fdug-up,]]>其中ug代表微元内的气体真实速度，up代表微元内的固体真实速度。
上述技术方案中，在所述的步骤2)中，所述的商业化计算流体力学软件采用CFX或Fluent或Phoenics或STAR-CD。
上述技术方案中，所述的颗粒流体两相流反应器是流态化煅烧炉或煤粉炉或循环流化床反应器或气固并流向上提升管反应器。
本发明的优点在于1、结构参数真正在微元层次与计算流体力学相耦合，得到了反映实际情况的相间动量交换系数或者曳力系数，纠正了常规方法中曳力系数偏高的情形；2、获得了任意微元内的结构分布，在此基础上可以正确计算颗粒流体两相流反应器内的反应、传热、传质过程，准确实现工业条件下两相流反应器的反应热态计算，从而指导相关过程控制、优化设计乃至工艺放大。
图1为本发明的测量颗粒流体两相流反应器内流动参数分布方法的流程图；图2为包含了结构因素的真实曳力系数β与商业化软件中均匀条件下的曳力系数β0的比值的示意图；图3为实施例2中的锅炉实体示意图；图4为在实施例2中利用本发明方法计算出的炉膛内压力分布图；图5为本发明方法计算的石化新型炼油多段反应器内分布板对颗粒浓度分布的影响；▲有分布板○无分布板图6为本发明方法计算的实验室中试规模冷态循环流态化变径反应器内流场分布图。
●预测结果■实验值▲商业软件计算值
具体实施例方式
下面结合附图和具体实施方式
对本发明的测量颗粒流体两相流反应器内流动参数分布的方法作详细说明。
颗粒流体两相流反应器是流态化煅烧炉、煤粉炉、循环流化床反应器和气固并流向上提升管反应器。本发明的方法是一种结合极值多尺度结构模型与双流体模型的计算方法，在对本发明所述方法进行说明前，首先对极值多尺度结构模型和双流体模型进行解释。
本发明方法中的极值多尺度结构模型是指模型通过在多个特征时间尺度和特征空间尺度上的物理量来计算结构参数和流动参数，物理量之间通过守恒方程和极值条件进行关联。其中，极值多尺度结构模型在所有尺度上所涉及的动力学守恒方程又采用双流体模型。
在颗粒流体两相流的任意微元控制体内都存在非均匀多尺度结构，这种结构可以等效为团聚物浓相和离散颗粒稀相两种双流体介质，及其界面结构并存的形式。各自的动量守恒方程经过转化，可以写成力平衡的形式。不失一般性，以气体颗粒两相流为例单位体积微元空间内，浓相团聚物颗粒曳力Fdc平衡方程Fd,c=34CDc1-ϵcdpρg[Uc-ϵcUpc/(1-ϵc)]2=(1-ϵg)(ρp-ρg)(ac+g)---(1)]]>单位体积微元空间内，稀相颗粒受曳力Fdf平衡方程Fd,f=34CDf1-ϵfdpρg[Uf-ϵfUpf/(1-ϵf)]2=(1-ϵf)(ρp-ρg)(af+g)---(2)]]>由于浓相团聚物界面处颗粒受力与其内部不同，可以得到单位体积微元空间内介观相界面的曳力Fdi平衡方程Fd,i=34CDifdclρg[Uf-ϵfUpc/(1-ϵc)]2(1-f)2=f(ϵg-ϵc)(ρp-ρg)(ai+g)---(3)]]>流体和颗粒相的质量守恒方程Ug＝fUc+(1-f)Uf(4)Up＝fUpc+(1-f)Upf(5)εg＝fεc+(1-f)εf(6)
其中，ρ，密度；Ug，气相表观速度；Up，颗粒相表观速度；Uf稀相流体速度；Upf，稀相颗粒速度；Uc，团聚物密相流体速度；Upc，团聚物密相颗粒速度；εg，微元内平均空隙率；εc，团聚物密相空隙率；εf，稀相空隙率；f，密相体积分数；dcl，团聚物直径；dp，单颗粒直径；ai，介观界面加速度；ac，密相加速度；af，稀相加速度；g，重力加速度；下标g，气相；下标p，固相；下标c，浓相；下标f，稀相。CDc，密相平均单颗粒曳力系数；CDf，稀相平均单颗粒曳力系数；Cdi，介观相间平均单个团聚物曳力系数。定义式如下CDf=(24+3.6Ref0.687)/(Ref·ϵf4.7),]]>Ref＝dpρgUsf/μgCD,c=(24+3.6Rec0.687)/(Rec·ϵc4.7),]]>Rec＝dpρgUsc/μgCD,i=(24+3.6Rei0.687)/[Rei·(1-f)4.7],]]>Rei＝dclρgUsi/μg各表观滑移速度定义为Usc＝Uc-εcUpc/(1-εc)Usf＝Uf-εfUpf/(1-εf)Usi＝Uf-εfUpc/(1-εc)(1-f)而团聚物直径dcl如下定义dcldp=(Nst)max-(Nst)mfNst-(Nst)mf---(7)]]>其中，下标mf代表εg＝εmf时的广义最小流化状态，εmf为最小流化状态时的空隙率，下标max代表εg＝εmax时的稀疏极限状态，εmax为存在颗粒聚团的最大空隙率，单位质量的颗粒悬浮输送能量Nst定义如下Nst=ρp-ρgρp[Ug-ϵf-ϵg1-ϵgf(1-f)Uf]---(8)]]>上述的公式(1)-(7)为动力学守恒方程，在这些方程中包含有11个待定参数Uf，Upf，Uc，Upc，εc，εf，f，dcl，ai，ac，af。单纯的动力学守恒方程不足以确定复杂两相流动中的结构因素，还需要极值多尺度结构模型中的极值条件。
通过对颗粒和流体两相运动机制的协调分析，发现流体运动的趋势是选择阻力最小的路径，而颗粒的运动是在边界条件控制的范围内趋于势能最小或者聚团。两种运动机制相互协调产生了极值多尺度结构模型的极值条件，即单位质量的颗粒悬浮输送能量趋于最小(Nst→min)。
这样，极值多尺度结构模型中11个流动及结构参数可以通过7个守恒方程和1个极值条件组成的极值问题得以数学求解。即在给定的物性和流动参数(Ug，Up)等条件下，求解满足动力学方程(1)-(7)的所有可能根中同时满足Nst→min的根，数学表达为MinNst(X)s.t.Fi(X)＝0，i＝1，2，……7. (9)如图l所示，本发明的测量颗粒流体两相流反应器内流动参数分布的方法包括以下步骤步骤10、根据实际的物性和反应器的操作条件情况，初始化流场及边界条件。所述的初始化流场及边界条件就是要根据实际的物性和操作条件给定颗粒流体两相流反应器的结构尺寸、颗粒粒径、颗粒密度、流体物性、反应温度、表观流体速度和颗粒流动通量等。
步骤20、利用双流体模型计算代码或商业化计算流体力学软件，计算每个空间微元内的质量守恒方程、动量守恒方程，得到速度、浓度分布。本步骤所述的计算流体力学软件可采用CFX，Fluent，Phoenics，STAR-CD等。速度和浓度分布可用(ug，up，εg)表示。ug和up代表微元内的气体和固体真实速度。
步骤30、判断当前的速度和浓度分布是否满足全场的收敛标准，若满足全场的收敛标准，输出当前的速度和浓度分布，并结束操作；否则，执行下一步。
步骤40、根据已知的反应器整体操作条件(Ug，Up)，计算并得到满足全局动力学的非稳态、非线性方程组的所有变量根的组合。在本步骤中，非稳态、非线性方程组是指动力学方程(1)-(7)，在这7个方程中，有11个变量，由7个方程得到变量根(Uf，Upf，Uc，Upc，εc，εf，f，dcl，ai，ac，af)的组合。
步骤50、从所有变量根的组合中，寻找满足极值条件的最优根，并保存其中的介观结构参数εc，εf，dcl。在本发明中，极值条件是指单位质量的颗粒悬浮输送能量趋于最小，即Nst→min。
步骤60、由步骤20得到的速度、浓度分布(ug，up，εg)和步骤50得到的介观结构参数εc，εf，dcl，在反应器的每个微元空间内求解满足微观局部动力学非稳态、非线性方程组(1)-(7)的所有变量根(Uf，Upf，Uc，Upc，f，ai，ac，af)的组合。
步骤70、由步骤60得到的变量根的组合，从中寻找满足极值条件的最优根。所述的极值条件是指单位质量的颗粒悬浮输送能量趋于最小，即Nst→min
步骤80、由步骤50和步骤70所得到的微元内所有变量值(Uf，Upf，Uc，Upc，εc，εf，f，dcl，ai，ac，af)的组合，计算曳力和曳力系数。曳力的计算公式如公式(10)所示。
Fd＝fFd，c+(1-f)Fd，f+Fd，i(10)根据所得到的曳力，还可进一步得到曳力系数β，曳力与曳力系数间的关系如公式(11)所示。
β=Fdug-up---(11)]]>无结构影响的曳力系数用β0表示β0=34ϵg(1-ϵg)dpρgCD0|ug-up|ϵg-2.7]]>步骤90、将步骤80得到的曳力系数代入步骤20中所采用的双流体模型计算代码或商业化计算流体力学软件中，替换原有的曳力系数，然后转到步骤20，做迭代计算。
在本发明方法计算得到的微元内浓稀相结构和流动参数的基础上，可以计算微元内浓、稀相各自的传热速率h、传质速度k和反应速率r，然后再进行加和。由于传热、传质和反应大多是非线性的，具有结构参数的计算结果符合实际，而均匀模型与实际不符。比如对反应而言，r实际＝r(浓相参数)+r(稀相参数)+r(界面参数)≠r(平均参数)将本发明方法耦合至质量、能量守恒方程相中，可以进一步得到颗粒流体两相流反应器中的质量浓度、热量分布的信息，从而更加准确得实现两相流反应器的反应热态计算。
由图2可知，包含了结构因素的真实曳力系数与商业化软件中均匀条件下的曳力系数不同，非均匀结构在大多数情况下都远低于均匀流动时的曳力。
本发明的测量颗粒流体两相流反应器内流动参数分布的方法在工业上有很好的应用。
实施例1将本发明方法应用于年产1000吨煤系高岭土快速流态化煅烧炉内的流动预测。原料高岭土平均粒径0.002mm，颗粒密度800kg/m3，煅烧炉膛内径0.45m，炉膛高10.5m。炉膛温度控制在900℃左右。发生的主要反应是C+O2→CO2，Al2O3·2SiO2·H2O→Al2O3·2SiO2。计算时在商业化CFD软件提供的通用平台上，将本发明方法按照计算框图计算的曳力系数编制成接口程序以代替软件内部的方法。首先对反应装置进行1∶1的几何构体，选择计算方式为双欧拉坐标下的双流体模型；其次设定温度并按照实际物性体系设定颗粒、气体性质；炉膛内物料初始设置为自然堆积状态，堆积空隙率为0.6，气体从炉膛底部入口进入，表观气速1.5m/s，颗粒物料由顶部侧面出口排出并循环至炉膛中部回料口。将计算接口文件调入和编译运行。计算按照与时间相关的动态模拟方式进行，计算实际流动时间20秒。本方法得到的空隙率分布、固体通量分布、固体速度分布具有动态非均匀结构，符合实际操作状态；而商业化软件计算得到的流场均匀，无明显的非均匀结构，这与实际情况不符。对于实际床层平均浓度为5％的情况，本发明方法计算的炉膛总压降(3557Pa)与商业化软件的计算结果(5608Pa)相比，更接近实际测试值(3058Pa)。
实施例2将本发明方法应用于每小时产蒸汽75吨的武汉石化厂循环流化床锅炉内的流动预测以辅助检测故障。如图3所示，锅炉炉膛主截面5.4×3.99m，炉膛高约22.5m，炉膛内循环灰颗粒平均粒径0.21mm，颗粒密度2000kg/m3，炉膛温度控制在900℃左右，设计一次风量49701Nm3/h，二次风量33134Nm3/h。炉膛内主要发生燃烧反应C+O2→CO2。在商业化CFD软件提供的通用平台上，将本发明方法计算的曳力系数写成用户自定义形式的外部接口程序以代替软件内部的方法。计算时首先对锅炉炉膛进行1∶1的几何构体，选择计算方式为双欧拉坐标下的双流体模型；其次设定温度并按照实际物性体系设定颗粒、气体性质；炉膛内物料初始设置为自然堆积状态，堆积空隙率为0.6，炉膛底部表观气速按照设计的一次风量计算，二次风通过在对应高度的气相守恒方程中加入质量源项来实现，颗粒物料由顶部侧面出口排出并循环至炉膛侧面回料口。调入曳力系数的计算接口文件并编译运行。计算实际流动时间60秒。图4显示了本发明方法计算出的炉膛内压力分布图，图中各点是在炉膛相应高度的截面上40个均匀分布的径向位置的压力值，全部点重合为一点说明该高度上不存在径向压力分布，40个径向位置数据重叠为一段线说明该高度上径向存在压力分布，流动混合剧烈、不稳定。本发明预测出的炉膛灰循环倍率为11，工厂操作现场长期测试总结的经验值在10左右，而商业化软件计算结果为46。本发明方法优于目前的商业化软件。
实施例3将本发明方法应用于中石化新型炼油多段反应器内流场分布的预测以辅助优化设计。反应器主体段直径3.4m，进口气速15m/s，进口固体量1500t/h，所用颗粒平均粒径0.075mm，颗粒密度1300kg/m3，炉膛操作温度约400℃。将本发明方法计算的曳力系数写成外部接口程序以代替商业化软件内部的方法。计算时首先对反应器进行1∶1的几何构体，选择计算方式为双欧拉坐标下的双流体模型；其次设定温度，并按照实际物性体系设定颗粒、气体性质；底部表观气速和固体速度按照设计值给定，颗粒物料由顶部出口排出。计算按照与时间相关的动态模拟方式进行，计算实际流动时间120秒。图5显示了本发明方法计算出的分布板对颗粒浓度分布的影响(▲-有分布板；○-无分布板)。结果显示反映其内的分布板影响反应器内物料分布形式，但是对扩径段反应器总体平均浓度(保持在2.4％左右)影响不大。这与工厂实际情况相符。
实施例4将本发明方法应用于中试规模循环流态化变径反应器内流场分布以放大设计工业规模反应器。反应器上段直径0.09m，下段直径0.045m，上段高约6.2m，颗粒平均粒径0.082mm，颗粒密度1450kg/m3，反应器底部进口气速10m/s，固体通量95kg/m2s，侧面进口固体通量12.6kg/m2s。将本发明方法计算的曳力系数写成外部接口程序以代替商业化软件内部的方法，并编译运行。图6显示本发明方法的预测结果(●)与实验值(■)吻合得很好，而现有的商业软件计算值(▲)预测结果偏稀，与实验值偏离较大。
权利要求
1.一种测量颗粒流体两相流反应器内流动参数分布的方法，包括1)、根据实际的物性和反应器的操作条件情况，初始化流场及边界条件；2)、利用双流体模型计算代码或商业化计算流体力学软件，计算每个空间微元内的质量守恒方程、动量守恒方程，得到速度、浓度分布(ug，up，εg)，所述的ug代表微元内的气体真实速度，所述的ug代表微元内的固体真实速度，所述的εg代表微元内平均空隙率；3)、判断当前的速度和浓度分布是否满足全场的收敛标准，若满足全场的收敛标准，输出当前的速度和浓度分布，并结束操作，否则，执行下一步；4)、根据已知的反应器整体操作条件(Ug，Up)，计算并得到满足全局动力学的非稳态、非线性方程组的所有变量根的组合；所述的Ug代表气相表观速度；所述的Up代表颗粒相表观速度；5)、从所有变量根的组合中，寻找满足极值条件的最优根，并保存其中的介观结构参数εc，εf，dcl，所述的εc表示团聚物密相空隙率，所述的εf表示稀相空隙率，所述的dcl表示团聚物直径；6)、由步骤2)得到的速度、浓度分布(ug，up，εg)和步骤5)得到的介观结构参数εc，εf，dcl，在反应器的每个微元空间内求解满足微观局部动力学非稳态、非线性方程组的所有变量根的组合；所述的变量根包括Uf，Upf，Uc，Upc，f，ai，ac，af，其中Uf代表稀相流体速度，Upf代表稀相颗粒速度，Uc代表团聚物密相流体速度，Upc代表团聚物密相颗粒速度，f代表密相体积分数，ai代表介观界面加速度，ac代表密相加速度，af代表稀相加速度；7)、由步骤6)得到的变量根的组合，从中寻找满足极值条件的最优根；8)、由步骤5)和步骤7)所得到的微元内所有变量值(Uf，Upf，Uc，Upc，εc，εf，f，dcl，ai，ac，af)的组合，计算曳力和曳力系数；9)、将步骤8)得到的曳力系数代入步骤2)中所采用的双流体模型计算代码或商业化计算流体力学软件中，替换原有的曳力系数，然后转到步骤2)，做迭代计算。
2.根据权利要求1所述的测量颗粒流体两相流反应器内流动参数分布的方法，其特征在于，所述的非稳态、非线性方程组包括单位体积微元空间内，浓相团聚物颗粒曳力Fdc平衡方程Fd,c=34CDc1-ϵcdpρg[Uc-ϵcUpc/(1-ϵc)]2=(1-ϵg)(ρp-ρg)(ac+g)]]>单位体积微元空间内，稀相颗粒受曳力Fdf平衡方程Fd,f=34CDf1-ϵfdpρg[Uf-ϵfUpf/(1-ϵf)]2=(1-ϵf)(ρp-ρg)(af+g)]]>单位体积微元空间内介观相界面的曳力Fdi平衡方程Fd,i=34CDifdclρg[Uf-ϵfUpc/(1-ϵc)]2(1-f)2=f(ϵg-ϵc)(ρp-ρg)(ai+g)]]>流体和颗粒相的质量守恒方程Ug＝fUc+(1-f)UfUp＝fUpc+(1-f)Upfεg＝fεc+(1-f)εf团聚物直径dcl公式dcldp=(Nst)max-(Nst)mfNst-(Nst)mf]]>其中，p，密度；Ug，气相表观速度；Up，颗粒相表观速度；Uf，稀相流体速度；Upf，稀相颗粒速度；Uc，团聚物密相流体速度；Upc，团聚物密相颗粒速度；εg，微元内平均空隙率；εc，团聚物密相空隙率；εf，稀相空隙率；f，密相体积分数；dcl，团聚物直径；dp，单颗粒直径；ai，介观界面加速度；ac，密相加速度；af稀相加速度；g，重力加速度；下标g，气相；下标p，固相；下标c，浓相；下标f，稀相。CDc，密相平均单颗粒曳力系数；CDf，稀相平均单颗粒曳力系数；Cdi，介观相间平均单个团聚物曳力系数；上述参数的定义式如下CD,f=(24+3.6Ref0.678)/(Ref·ϵf4.7),]]>Ref＝dpρgUsf/μgCD,c=(24+3.6Rec0.678)/(Rec·ϵc4.7),]]>Rec＝dpρgUsc/μgCD,i=(24+3.6Rei0.678)/[Rei·(1-f)4.7],]]>Rei＝dclρgUsi/μg各表观滑移速度定义为Usc＝Uc-εcUpc/(1-εc)Usf＝Uf-εfUpf/(1-εf)Usi＝Uf-εfUpc/(1-εc)(1-f)。
3.根据权利要求1所述的测量颗粒流体两相流反应器内流动参数分布的方法，其特征在于，所述的步骤5)和步骤7)中，所述的极值条件是指单位质量的颗粒悬浮输送能量趋于最小，所述的单位质量的颗粒悬浮输送能量用Nst表示，其定义为Nst=ρp-ρgρp[Ug-ϵf-ϵg1-ϵgf(1-f)Uf].]]>
4.根据权利要求1所述的测量颗粒流体两相流反应器内流动参数分布的方法，其特征在于，所述的步骤8)中，所述的曳力用Fd表示，其计算公式为Fd＝fFd，c+(1-f)Fd，f+Fd，i；所述的曳力系数用β表示，可由曳力得到，所述的曳力系数的计算公式为β=Fdug-up,]]>其中ug代表微元内的气体真实速度，up代表微元内的固体真实速度。
5.根据权利要求1所述的测量颗粒流体两相流反应器内流动参数分布的方法，其特征在于，在所述的步骤2)中，所述的商业化计算流体力学软件采用CFX或Fluent或Phoenics或STAR-CD。
6.根据权利要求1所述的测量颗粒流体两相流反应器内流动参数分布的方法，其特征在于，所述的颗粒流体两相流反应器是流态化煅烧炉或煤粉炉或循环流化床反应器或气固并流向上提升管反应器。
全文摘要
本发明公开了一种测量颗粒流体两相流反应器内流动参数分布的方法，包括根据反应器的结构尺寸和物性、操作条件，由双流体模型计算每个空间微元内的守恒方程得到速度、浓度分布；由反应器整体操作条件计算满足全局动力学的非稳态方程组的所有变量根的组合，并寻找满足极值条件的最优根以得到介观结构参数；在反应器每个微元空间内求解满足微观局部动力学方程组的所有变量根的组合，并寻找满足极值条件的最优根；计算颗粒流体相间曳力，并替代双流体模型的曳力关系做迭代计算，得到反应器任意微元内的浓度、速度流场的空间分布和时间演化。本发明方法可得到反映实际情况的曳力系数，可准确计算反应器内的流动、反应、传热、传质过程。
文档编号G01M10/00GK1979111SQ20051012622
公开日2007年6月13日 申请日期2005年11月30日 优先权日2005年11月30日
发明者王维, 葛蔚, 杨宁, 李静海 申请人:中国科学院过程工程研究所