本发明涉及颗粒流体两相流反应技术领域，尤其涉及一种颗粒流体两相流反应器内流动参数分布的测量方法。
背景技术：
颗粒流体两相流广泛存在于能源、化工、冶金和环境工业中。颗粒流体两相流中的非线性非平衡相互作用，如颗粒和流体间的曳力相互作用、颗粒与颗粒之间的碰撞和摩擦等导致两相流动中呈现出气泡和聚团等复杂的非均匀结构。这些非均匀结构对于工业反应器内的压降分布、颗粒通量、传质和传热效率及反应转化率等有明显的影响。传统未考虑非均匀结构的气固相互作用曳力系数通常比真实的曳力系数高1至3个数量级。因此，准确预测颗粒流体两相流反应器内的流动结构和参数分布对颗粒流体两相流的设计、运行、控制与放大等至关重要。两相流的相间曳力预测方法是实现这一工作目标的核心内容之一。
传统的颗粒流体两相流中的相间曳力系数通常采用基于固定床压降外推关系的ergun曳力关系式或者基于散式流态化的实验经验关联的wen&yu曳力系数。这种计算方式被绝大多数研究者和实践者采纳，并被国际上通行的计算流体力学(cfd)软件ansyscfx，ansysfluent，phoenics，star-cd等推荐。这些方法由于未考虑非均匀结构，仅适用于散式流态化体系。对于气体颗粒的流态化，尤其是细颗粒的流态化行为，如石油工业中常见的流化催化裂化(fcc)反应器，电力行业中的循环流化床锅炉(cfbb)等核心单元内的流动，体系中非均匀结构明显。在此情况下，不考虑结构影响的计算值与实际值相比，差异很大，结果处于不同的数量级，其直接后果是设计或者仿真计算时的流场颗粒浓度及反应器内压降分布失真。
cn101285753a公开了一种测量颗粒流体两相流反应器内流体相组分浓度分布的方法，通过将多尺度流动结构引入到颗粒流体传质测量过程中，充分考虑了局部传质的不均衡性，但该方法计算考虑了加速度的影响，且该方法计算复杂，通用性差。
cn101042310a公开了一种测量颗粒流体两相流系统非均匀结构参数分布的方法，所述方法通过关联不同尺度的相互作用，一定程度上提高了计算速度，但该方法中仍然考虑了加速度的影响，不仅计算量大、通用性差且与实际操作情况差距仍然较大。
cn1979111a公开了一种基于能量最小多尺度方法获取非均匀结构曳力，并用于两相流反应器中计算流体力学分析的方法，但该方法在守恒方程中引入了加速度等变量构成非稳态非线性方程组，并通过两次使用极值条件来获取最优解。这一方法虽然能获取考虑结构的非均匀曳力系数，然后计算方法复杂，且只能获取某一操作条件下的非均匀曳力系数，不具有通用性。同时，这一方法仅使用于循环床流域，不适用于鼓泡床等流域的计算，在湍动床流域的计算精度低，不具有流域的普适性。
因此，亟需开发一种新的颗粒流体两相流反应器内流动参数分布的测量方法，提高现有测量方法的计算速度、准确性和通用性。
技术实现要素：
鉴于现有技术中存在的问题，本发明提供一种颗粒流体两相流反应器内流动参数分布的测量方法，所述方法在两相流计算过程中引入结构的复杂影响，将极值多尺度结构模型与多相流模型相结合，并基于力平衡构建了稳态下的非线性方程组，不再引入加速度等变量，不仅计算过程更简便，而且具有操作条件的通用性和流域的普适性，同时具有计算精度高的优点，能够准确预测两相流反应器内的流动分布，并掌握反应器内部规律，从根本上保障设计、运行、控制与放大。
为达此目的，本发明采用以下技术方案：
第一方面，本发明提供一种颗粒流体两相流反应器内流动参数分布的测量方法，所述方法包括如下步骤：
(1)根据流体的物性参数和操作条件，初始化流场和边界条件；
(2)计算各空间微元内的质量守恒方程和动量守恒方程，读取所述空间微元内的速度分布和浓度分布；
(3)根据反应器内的气相表观速度ug和颗粒相表观速度up，计算得到满足稳态下质量非线性方程组和动量非线性方程组的变量的根的组合，所述变量中不包括加速度；
(4)从步骤(3)所述变量的根的组合中找到满足极值条件的最优根，并根据最优根计算非均匀结构因子和曳力系数；
(5)将步骤(4)得到的曳力系数代入步骤(2)中替换原有的曳力系数，转到步骤(2)进行迭代计算，所述迭代计算的终止条件为：所述速度分布和浓度分布满足质量非线性方程组和动量非线性方程组的收敛标准。
本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法将极值多尺度结构模型与多相流模型相结合，在现有技术的基础上基于力平衡构建了稳态下的非线性方程组，步骤(3)中变量为(ugf，upf，ugc，upc，εgc，εgf，f)，不再引入加速度等变量，发明人发现，在稳态条件下加速度的影响可以忽略，同时现有加速度的计算与实际情况相差较大，本发明不再引入加速度等变量，不仅计算过程更简便，而且提高了计算的精度和普适性，更有利于反应器的设计与放大。
本发明中极值多尺度结构模型是指：模型通过在多个特征时间尺度和特征空间尺度上的物理量来计算结构参数和流动参数，物理量以质量和力平衡方程为约束，通过极值条件来获取最优解。
本发明对多相流模型没有特殊限制，可采用本领域技术人员熟知的任何可用于计算颗粒流体两相流的模型，例如可以是双流体模型、离散颗粒模型或多相流质点网格模型等。
优选地，步骤(2)中所述计算的工具包括多相流模型计算代码或商业化计算流体力学软件。
优选地，所述商业化计算流体力学软件采用ansyscfx、ansysfluent、phoenics、star-cd、barracuda、openfoam或mfix中的任意一种。
优选地，所述质量守恒方程为计算工具中的质量守恒方程。
优选地，所述动量守恒方程为计算工具中的动量守恒方程。
优选地，步骤(2)中所述空间微元内的速度分布和浓度分布用参数ug，up和εg表示。
其中，ug表示空间微元内的流体真实速度，up表示空间微元内的固体真实速度，εg表示微元内平均空隙率。
优选地，步骤(3)中所述稳态下的质量非线性方程组包括方程组(1)：
其中，f表示密相体积分数，ugc表示团聚物密相流体速度，ugf表示稀相流体速度，upc表示团聚物密相颗粒速度，upf表示稀相颗粒速度，εgc表示团聚物密相空隙率，εgf表示稀相空隙率。
所述稳态下的动量非线性方程组包括方程组(2)：
其中，mc表示密相颗粒数密度，fdc表示密相曳力，mi表示聚团数密度，fdi表示相间曳力，ρp表示颗粒密度，ρg表示气体密度，mf表示稀相颗粒数密度，fdf表示稀相曳力，g表示重力加速度。
所述密相曳力fdc的定义式为：
所述相间曳力fdi的定义式为：
所述稀相曳力fdf的定义式为：
其中，cdc表示密相曳力系数，cdi表示相间曳力系数，cdf表示稀相曳力系数，usc表示密相表观滑移速度，usi表示相间表观滑移速度，usf表示稀相表观滑移速度，dcl表示团聚物直径，dp表示颗粒直径。
颗粒直径dp是指反应器中颗粒本身的直径，团聚物直径是指颗粒在反应过程中发生团聚形成的团聚物的直径。
所述密相曳力系数cdc的定义式为：
所述相间曳力系数cdi的定义式为：
所述稀相曳力系数cdf的定义式为：
所述密相表观滑移速度usc的定义式为：
所述相间表观滑移速度usi的定义式为：
所述稀相表观滑移速度usf的定义式为：
所述团聚物直径dcl的定义式为：
在一般的流态化系统中，团聚物的大小与输入体系的能量成反比。当床层处于广义最小流化状态时，近似认为团聚物直径为无穷大，而在床层处于最大空隙率时，团聚物消失，团聚物的当量直径等于颗粒直径，由此构建团聚物直径方程。
其中，umf表示初始流化速度，εmf表示初始流化空隙率，εmax表示最大空隙率，rec表示密相雷诺数，rei表示相间雷诺数，ref表示稀相雷诺数。
优选地，所述最大空隙率为εmax＝0.9997。
所述密相雷诺数rec的定义式为：
所述相间雷诺数rei的定义式为：
所述稀相雷诺数ref的定义式为：
其中，μg表示气相黏度。
优选地，步骤(4)中所述非均匀结构因子的定义式为：
其中，ar表示阿基米德数，其定义式为：
re表示雷诺数，其定义式为：re＝ρgdp(ug-up/(1-εg))/μg。
优选地，步骤(4)中所述曳力系数β的计算公式为：β＝hdβ0。
其中，β0表示均匀曳力系数，其计算公式为：
优选地，步骤(4)中所述极值条件为单位质量的颗粒悬浮输送能量最小。
所述单位质量的颗粒悬浮输送能量nst的定义式为式(3)：
本发明通过对颗粒和流体两相运动机制的协调分析，发现流体运动的趋势是选择阻力最小的路径，而颗粒的运动是在边界条件控制的范围内趋于势能最小。两种运动机制相互协调产生了极值多尺度结构模型的极值条件，即：单位质量的颗粒悬浮输送能量趋于最小(nst→min)。
这样，极值多尺度结构模型中7个流动及结构参数(ugf，upf，ugc，upc，εgc，εgf，f)可以通过6个守恒方程，即方程组(1)和方程组(2)，和1个极值条件方程(3)组成的极值问题得以数学求解。即：在给定的物性和流动参数(ug，up)等条件下，求解满足动力学方程组(1)和方程组(2)的所有可能根中同时满足nst→min的根，数学表达为：
minnst(x)s.t.fi(x)＝0，i＝1，2，3，4，5，6。
其中，fi(x)＝0，i＝1，2，3，4，5，6为方程组(1)和方程组(2)中的6个方程；x为变量矩阵，表达式为x＝[ugf，upf，ugc，upc，εgc，εgf，f]’。
优选地，步骤(5)中所述收敛标准为质量非线性方程组和动量非线性方程组的计算残差小于收敛值。
优选地，所述收敛值为0～0.1，例如可以是0、0.00001、0.0001、0.001、0.01或0.1，优选为0～0.001。
优选地，在步骤(5)终止迭代计算之后，转到步骤(2)重新读取空间微元内的速度分布和浓度分布，再次进行循环计算。
优选地，所述循环计算的结束标准为达到计算总时长。
本发明中由于反应器中由初始状态运行至稳定需要一段时间，因此在迭代计算的基础上随着时间的推移不断计算各空间微元内的质量守恒方程和动量守恒方程，读取所述空间微元内的速度分布和浓度分布，循环计算至反应器中达到稳态为止。
优选地，所述计算总时长为0.1～500s，例如可以是0.1s、0.2s、1s、2s、3s、4s、5s、6s、7s、8s、9s、10s、12s、15s、20s、30s、40s、50s、60s、70s、80s、90s、100s、200s、300s、400s或500s，优选为5～100s。
优选地，所述颗粒流体两相流反应器包括流态化煅烧炉、煤粉炉、循环流化床反应器或流固并流向上提升管反应器中的任意一种。
本发明中相同符号代表的物理量及物理意义相同。
第二方面，本发明提供一种计算机可读存储介质，所述存储介质上储存有根据第一方面所述的颗粒流体两相流反应器内流动参数分布的测量方法编写的程序
本发明提供的计算机可读存储介质中储存有根据第一方面所述的颗粒流体两相流反应器内流动参数分布的测量方法编写的程序，能够将该计算机存储介质与现有商业化计算流体力学的软件对接，从而实现对颗粒流体两相流反应器内流动参数分布的预测。
与现有技术相比，本发明至少具有以下有益效果：
(1)本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法考虑了空间微元中的结构对相间动量交换系数或曳力系数的影响，纠正了现有方法中曳力系数偏高的情形；
(2)本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法适用的流域范围更广、具有操作条件的普适性，计算过程更简便且计算精度更高；
(3)本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法能够合理计算颗粒流体两相流反应器内的流动过程，准确实现两相流反应器的流动计算，合理预测反应器的颗粒浓度分布，从而指导相关过程控制、优化设计乃至工艺放大等过程。
附图说明
图1是本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法的流程示意图。
图2是本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法中非均匀结构因子随微元内平均空隙率εg和气相表观速度ug变化的示意图。
图3是本发明应用例2中计算得到的反应器颗粒浓度分布图。
图4是本发明应用例3中计算得到的流化床气化反应器颗粒浓度分布图。
图5是本发明应用例1、应用例对比例1和应用对比例2中湍动床的颗粒分布预测结果与实验值的对比图。
具体实施方式
下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。
下面对本发明进一步详细说明。但下述的实例仅仅是本发明的简易例子，并不代表或限制本发明的权利保护范围，本发明的保护范围以权利要求书为准。
本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法的流程示意图如图1所示，所述方法包括如下步骤：
(1)根据流体的物性参数和操作条件，初始化流场和边界条件；
(2)计算各空间微元内的质量守恒方程和动量守恒方程，读取所述空间微元内的速度分布和浓度分布ug，up和εg；
(3)根据反应器内的操作条件即气相表观速度ug和颗粒相表观速度up，计算得到满足稳态下质量非线性方程组和动量非线性方程组的变量的根的组合，其中，所述变量为(ugf，upf，ugc，upc，εgc，εgf，f)，所述变量中不包括加速度；
(4)从步骤(3)所述变量的根的组合中找到满足nst最小的最优根，并根据最优根计算非均匀结构因子和曳力系数；
(5)将步骤(4)得到的曳力系数代入步骤(2)中替换原有的曳力系数，转到步骤(2)进行迭代计算，所述迭代计算的终止条件为：所述速度分布和浓度分布满足质量非线性方程组和动量非线性方程组的收敛标准；
(6)在步骤(5)终止迭代计算之后，转到步骤(2)重新读取空间微元内的速度分布和浓度分布，再次进行循环计算至达到计算总时长时判断计算结束。
本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法非均匀结构因子随微元内平均空隙率εg和气相表观速度ug变化的示意图如图2所示，从图2可以看出，本发明中包含了结构因素的曳力系数与商业化软件中均匀条件下的曳力系数不同，非均匀结构在大多数情况下都远低于均匀流动时的曳力。
一、实施例
实施例1
本实施例提供一种颗粒流体两相流反应器内流动参数分布的测量方法，所述方法包括如下步骤：
(1)根据流体的物性参数和操作条件，初始化流场和边界条件；所述的初始化流场及边界条件就是要根据实际的物性和操作条件给定颗粒流体两相流反应器的结构尺寸、颗粒粒径、颗粒密度、流体密度、流体粘度、表观流体速度和颗粒流动通量；
(2)计算各空间微元内的质量守恒方程和动量守恒方程，读取所述空间微元内的速度分布和浓度分布，所述空间微元内的速度分布和浓度分布用参数ug，up和εg表示，其中，ug表示空间微元内的流体真实速度，up表示空间微元内的固体真实速度，εg表示微元内平均空隙率；
(3)根据反应器内的气相表观速度ug和颗粒相表观速度up，计算得到满足稳态下质量非线性方程组和动量非线性方程组的变量的根的组合，所述变量中不包括加速度；所述稳态下的质量非线性方程组包括方程组(1)：
其中，f表示密相体积分数，ugc表示团聚物密相流体速度，ugf表示稀相流体速度，upc表示团聚物密相颗粒速度，upf表示稀相颗粒速度，εgc表示团聚物密相空隙率，εgf表示稀相空隙率；
所述稳态下的动量非线性方程组包括方程组(2)：
其中，mc表示密相颗粒数密度，fdc表示密相曳力，mi表示聚团数密度，fdi表示相间曳力，ρp表示颗粒密度，ρg表示气体密度，mf表示稀相颗粒数密度，fdf表示稀相曳力，g表示重力加速度；
所述密相曳力fdc的定义式为：
所述相间曳力fdi的定义式为：
所述稀相曳力fdf的定义式为：
其中，cdc表示密相曳力系数，cdi表示相间曳力系数，cdf表示稀相曳力系数，usc表示密相表观滑移速度，usi表示相间表观滑移速度，usf表示稀相表观滑移速度，dcl表示团聚物直径，dp表示颗粒直径；
所述密相曳力系数cdc的定义式为：
所述相间曳力系数cdi的定义式为：
所述稀相曳力系数cdf的定义式为：
所述密相表观滑移速度usc的定义式为：
所述相间表观滑移速度usi的定义式为：
所述稀相表观滑移速度usf的定义式为：
所述团聚物直径dcl的定义式为：
其中，umf表示初始流化速度，εmf表示初始流化空隙率，εmax表示最大空隙率(εmax＝0.9997)，rec表示密相雷诺数，rei表示相间雷诺数，ref表示稀相雷诺数；
所述密相雷诺数rec的定义式为：
所述相间雷诺数rei的定义式为：
所述稀相雷诺数ref的定义式为：
其中，μg表示气相黏度；
(4)从步骤(3)所述变量的根的组合中找到满足极值条件的最优根，并根据最优根计算非均匀结构因子和曳力系数；其中，所述非均匀结构因子的定义式为：
其中，ar表示阿基米德数，其定义式为：
re表示雷诺数，其定义式为：re＝ρgdp(ug-up/(1-εg))/μg。
优选地，所述曳力系数β的计算公式为：β＝hdβ0；
其中，β0表示均匀曳力系数，其计算公式为：
所述极值条件为单位质量的颗粒悬浮输送能量最小；所述单位质量的颗粒悬浮输送能量nst的定义式为式(3)：
(5)将步骤(4)得到的曳力系数代入步骤(2)中替换原有的曳力系数，转到步骤(2)进行迭代计算，所述迭代计算的终止条件为：所述速度分布和浓度分布满足质量非线性方程组和动量非线性方程组的收敛标准0.1，并保存计算得到的非均匀结构因子和曳力系数。
二、对比例
对比例1
本对比例提供一种颗粒流体两相流反应器内流动参数分布的测量方法，所述方法直接采用均匀曳力测量颗粒流体两相流反应器内流动参数分布。
对比例2
本对比例提供一种颗粒流体两相流反应器内流动参数分布的测量方法，所述方法采用cn1979111a中提供的颗粒流体两相流反应器内流动参数分布的测量方法。
三、应用例
应用例1
本应用例将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法应用于气固两相的湍动床中。
计算的工具：商业化cfd软件ansysfluent提供的通用平台
流体的物性参数和操作条件：颗粒的平均粒径0.09mm，颗粒密度1374kg/m3，气体的表观速度为0.6m/s。
反应器尺寸和结构：湍动床下部主体为一个圆管结构，内径0.05m，高0.75m，上部为一段管径拓宽区域，内径0.1m，高0.3m。
所述预测方法包括如下步骤：
(1)将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法编织成接口程序以代替软件内部的方法；
(2)对反应器进行1:1的几何构体，选择计算模型为双欧拉坐标下的双流体模型；其次设定温度并按照实际物性体系设定颗粒、气体性质；湍动床内颗粒初始设置为自然堆积状态，堆积高度0.2m，堆积浓度为0.56；气体从湍动床底部入口进入，颗粒被带入到上部管径拓宽区域，发生颗粒反混；
(3)将步骤(1)中计算曳力系数的接口文件调入和编译运行，计算按照与时间相关的动态模拟方式进行，循环计算的总时长为20s。
应用例2
本应用例将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法应用于中石油催化剂的再生反应器中。
计算的工具：商业化cfd软件ansysfluent提供的通用平台
流体的物性参数和操作条件：反应器的温度设置在660℃，操作压力为210kpa。颗粒的平均粒径0.056mm，颗粒密度1465kg/m3，气体的密度为0.782kg/m3，气体的粘度为3.719×10-5pa·s。
反应器尺寸和结构：反应器内包含一个具有复杂分布板结构，它是再生器的一个关键部件，以实现均匀布气，达到均匀流化；反应器的下部一区直径6.2m，高度15.1m，上部二区直径11.3m，高度21.7m；回料管直径1.42m，高度16.1m，进气管直径1.2。
所述应用方法包括如下步骤：
(1)将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法编织成接口程序以代替软件内部的方法；
(2)对反应器进行1:1的几何构体，选择计算模型为双欧拉坐标下的双流体模型；其次设定温度并按照实际物性体系设定颗粒、气体性质；湍动床内颗粒初始设置为自然堆积状态，堆积空隙率为0.5；
(3)将步骤(1)中计算曳力系数的接口文件调入和编译运行，计算按照与时间相关的动态模拟方式进行，循环计算的总时长为80s。
本应用例计算得到的反应器颗粒浓度分布图如图3所示，其结果显示流动模拟合理地捕捉到了床层内的颗粒聚团等介尺度结构。颗粒在床层轴向上呈上稀下浓的“s”型分布，下部浓相区与底部稀相区有明显的界线。同时颗粒在径向上呈现边壁浓中间稀的，合理捕捉到了的颗粒的壁面返混现象。计算得到的颗粒分布接近真实的流化行为。
应用例3
本应用例将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法应用于中科院山西煤化所开发的灰融聚多段分级转化流化床气化反应器。
计算的工具：商业化cfd软件barracuda提供的通用平台
流体的物性参数和操作条件：反应器温度为1000℃，颗粒密度为1200kg/m3，中心管气体流量为0.428kg/s，煤处理量为0.514kg/s。
反应器尺寸和结构：该反应器通过多段结构来组织气化反应过程，实现最优的气化效率；反应器的底部鼓泡床直径为0.8m，高度为4.4m，上部快速床直径为0.4m，高度为13m，过渡段高度为0.8m，一级旋风分离器的直径为0.52m；低床层初始堆料高度为1.02m，高床层初始堆料高度为0.68m，堆料浓度为0.5。
所述应用方法包括如下步骤：
(1)将实施例1提供的颗粒流体两相流反应器内流动参数分布的测量方法编织成接口程序以代替软件内部的方法；
(2)对反应器进行1:1的几何构体，选择计算模型为mp-pic模型；其次设定温度并按照实际物性体系设定颗粒、气体性质；底部表观气速和固体速度按照设计值给定，颗粒物料由顶部出口排出；
(3)将步骤(1)中计算曳力系数的接口文件调入和编译运行，计算按照与时间相关的动态模拟方式进行，循环计算的总时长为200s。
本应用例计算得到的流化床气化反应器颗粒浓度分布图如图4所示，该图表明，本应用例可准确预测不同反应段所特有的流动结构以及颗粒分布状态，其中底部鼓泡床颗粒浓度较高，颗粒的返混严重；在过渡段颗粒浓度迅速下降；到达上部快速床区域，颗粒的浓度较低。
应用对比例1
本应用例将对比例1提供的颗粒流体两相流反应器内流动参数分布的测量方法应用于气固两相的湍动床中，所述应用方法除不进行步骤(1)，直接调入均匀曳力系数外，其余均与应用例1相同。
应用对比例2
本应用例将对比例2提供的颗粒流体两相流反应器内流动参数分布的测量方法应用于气固两相的湍动床中，所述应用方法除将颗粒流体两相流反应器内流动参数分布的测量方法替换为对比例2提供的方法外，其余均与应用例1相同。
本发明应用例1、应用对比例1和应用对比例2中湍动床的颗粒分布预测结果与实验值的对比图如图5所示，从图可以看出，应用例1提供的方法通过采用实施例1提供的方法对湍动床的颗粒分布进行预测，较应用对比例1和应用对比例2分别采用对比例1和对比例2提供的方法进行预测而言，应用例1计算得到的颗粒浓度分布具有动态非均匀结构，更符合实际操作状态；而应用对比例1和应用对比例2得到的流场颗粒浓度轴向分布与实际情况不符，由此说明本发明提供的方法通过将极值多尺度结构模型与多相流模型相结合，且不考虑加速度的影响，计算得到的湍动床内颗粒浓度的分布与实际情况更相符。
综上，本发明提供的颗粒流体两相流反应器内流动参数分布的测量方法通过不再引入加速度等变量，只需通过使用一次极值条件即可获取最优解，不仅计算过程更简便，而且计算精度高，与实验结果更相符，能够准确预测两相流反应器内的流动分布，并掌握反应器内部规律，从根本上保障设计、运行、控制与放大。
申请人声明，本发明通过上述实施例来说明本发明的详细结构特征，但本发明并不局限于上述详细结构特征，即不意味着本发明必须依赖上述详细结构特征才能实施。所属技术领域的技术人员应该明了，对本发明的任何改进，对本发明所选用部件的等效替换以及辅助部件的增加、具体方式的选择等，均落在本发明的保护范围和公开范围之内。
1.一种颗粒流体两相流反应器内流动参数分布的测量方法，其特征在于，所述方法包括如下步骤：
(1)根据流体的物性参数和操作条件，初始化流场和边界条件；
(2)计算各空间微元内的质量守恒方程和动量守恒方程，读取所述空间微元内的速度分布和浓度分布；
(3)根据反应器内的气相表观速度ug和颗粒相表观速度up，计算得到满足稳态下质量非线性方程组和动量非线性方程组的变量的根的组合，所述变量中不包括加速度；
(4)从步骤(3)所述变量的根的组合中找到满足极值条件的最优根，并根据最优根计算非均匀结构因子和曳力系数；
(5)将步骤(4)得到的曳力系数代入步骤(2)中替换原有的曳力系数，转到步骤(2)进行迭代计算，所述迭代计算的终止条件为：所述速度分布和浓度分布满足质量非线性方程组和动量非线性方程组的收敛标准。
2.根据权利要求1所述的方法，其特征在于，步骤(2)中所述计算的工具包括多相流模型计算代码或商业化计算流体力学软件；
优选地，所述质量守恒方程为计算工具中的质量守恒方程；
优选地，所述动量守恒方程为计算工具中的动量守恒方程。
3.根据权利要求1或2所述的方法，其特征在于，步骤(2)中所述空间微元内的速度分布和浓度分布用参数ug，up和εg表示；
其中，ug表示空间微元内的流体真实速度，up表示空间微元内的固体真实速度，εg表示微元内平均空隙率。
4.根据权利要求1～3任一项所述的方法，其特征在于，步骤(3)中所述稳态下的质量非线性方程组包括方程组(1)：
其中，f表示密相体积分数，ugc表示团聚物密相流体速度，ugf表示稀相流体速度，upc表示团聚物密相颗粒速度，upf表示稀相颗粒速度，εgc表示团聚物密相空隙率，εgf表示稀相空隙率；
所述稳态下的动量非线性方程组包括方程组(2)：
其中，mc表示密相颗粒数密度，fdc表示密相曳力，mi表示聚团数密度，fdi表示相间曳力，ρp表示颗粒密度，ρg表示气体密度，mf表示稀相颗粒数密度，fdf表示稀相曳力，g表示重力加速度；
所述密相曳力fdc的定义式为：
所述相间曳力fdi的定义式为：
所述稀相曳力fdf的定义式为：
其中，cdc表示密相曳力系数，cdi表示相间曳力系数，cdf表示稀相曳力系数，usc表示密相表观滑移速度，usi表示相间表观滑移速度，usf表示稀相表观滑移速度，dcl表示团聚物直径，dp表示颗粒直径；
所述密相曳力系数cdc的定义式为：
所述相间曳力系数cdi的定义式为：
所述稀相曳力系数cdf的定义式为：
所述密相表观滑移速度usc的定义式为：
所述相间表观滑移速度usi的定义式为：
所述稀相表观滑移速度usf的定义式为：
所述团聚物直径dcl的定义式为：
其中，umf表示初始流化速度，εmf表示初始流化空隙率，εmax表示最大空隙率，rec表示密相雷诺数，rei表示相间雷诺数，ref表示稀相雷诺数。
5.根据权利要求1～4任一项所述的方法，其特征在于，步骤(4)中所述非均匀结构因子的定义式为：
其中，ar表示阿基米德数，其定义式为：
re表示雷诺数，其定义式为：re＝ρgdp(ug-up/(1-εg))/μg。
优选地，所述曳力系数β的计算公式为：β＝hdβ0；
其中，β0表示均匀曳力系数，其计算公式为：
6.根据权利要求1～5任一项所述的方法，其特征在于，步骤(4)中所述极值条件为单位质量的颗粒悬浮输送能量最小；
所述单位质量的颗粒悬浮输送能量nst的定义式为式(3)：
7.根据权利要求1～6任一项所述的方法，其特征在于，步骤(5)中所述收敛标准为质量非线性方程组和动量非线性方程组的计算残差小于收敛值；
优选地，所述收敛值为0～0.1，优选为0～0.001。
8.根据权利要求1～7任一项所述的方法，其特征在于，在步骤(5)终止迭代计算之后，转到步骤(2)重新读取空间微元内的速度分布和浓度分布，再次进行循环计算；
优选地，所述循环计算的结束标准为达到计算总时长；
优选地，所述计算总时长为0.1～500s，优选为5～100s。
9.根据权利要求1～8任一项所述的方法，其特征在于，所述颗粒流体两相流反应器包括流态化煅烧炉、煤粉炉、循环流化床反应器或流固并流向上提升管反应器中的任意一种。
10.一种计算机可读存储介质，其特征在于，所述存储介质上储存有根据权利要求1～9任一项中所述的颗粒流体两相流反应器内流动参数分布的测量方法编写的程序。
技术总结