CN104634708A - 基于计算流体力学流化床内颗粒密度和粒径分布预测方法 - Google Patents

Abstract

基于计算流体力学流化床内颗粒密度和粒径分布预测方法，步骤一、流化床内基本流动反应模型的建立；步骤二、建立描述颗粒相密度和粒径变化规律的数学模型；步骤三、流化床内密度和粒径分布状态的预测；该方法采用计算流体力学方法对流化床进行模拟研究，并结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，从而准确预测流化床内颗粒密度和粒径的分布状态，为流化床的性能预测、优化控制及设计放大提供理论基础。该方法通过计算流体力学方法获得流化床内颗粒密度和粒径分布状态，无需对实际运行流化床进行复杂的采样分析，节约大量的人力、物力及时间成本，是一种获得流化床内颗粒密度和粒径分布状态的新方法。

Description

基于计算流体力学流化床内颗粒密度和粒径分布预测方法

技术领域

本发明涉及一种获得流化床内颗粒密度和粒径分布状态的新方法，具体涉及一种基于计算流体力学流化床内颗粒密度和粒径分布预测方法。

背景技术

流化床内的固体颗粒在流体的作用下具有流体某些表观特性，颗粒表面能全部暴露于周围剧烈湍动的流体中，从而强化了传热、传质和化学反应，因此广泛用于工业领域如生物质、煤等固体燃料的热解、燃烧或气化等过程。生物质或煤等颗粒在流化床内发生热解、燃烧或气化等非均相化学反应时往往遵循一定的密度和粒径变化规律，从而表现出一定的密度和粒径分布状态。而流化床内的流动反应特性与颗粒的密度和粒径分布有着密切的关系，因此准确获得流化床内颗粒密度和粒径的分布状态可为流化床的性能预测、优化控制及设计放大提供理论基础。通常可以通过实验取样分析获得流化床内颗粒密度和粒径分布状态，但是该方法需耗费大量人力、物力及时间成本，并且现阶段难以实现对高温流化床进行取样分析试验，因此流化床特别是实际运行的高温流化床内的颗粒密度和粒径分布状态往往难以获得。

而计算流体力学是对研究对象建立流动反应模型并采用计算机和离散化的数值方法进行数值模拟和分析的一种方法。该方法不受现有实验技术的限制，可以全面、高效且低成本地揭示流化床内的流动反应特性，因此广泛应用于流化床的研究当中。然而，在现有的流化床计算流体力学研究中往往忽略颗粒发生非均相化学反应时的密度和粒径变化，虽然最近有研究考虑了颗粒密度的变化规律，但是由于未能考虑颗粒相应的粒径变化而得不到合理的结果，因此，现有计算流体力学方法仍无法准确获得流化床内颗粒密度和粒径的分布状态。

发明内容

为了解决现有技术存在的问题，本发明的目的在于提供一种基于计算流体力学流化床内颗粒密度和粒径分布预测方法，该方法采用计算流体力学方法对流化床进行模拟研究，并结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，从而获得流化床内颗粒密度和粒径的分布状态。该方法无需对实际运行流化床进行复杂的采样分析，节约大量的人力、物力及时间成本，是一种获得流化床内颗粒密度和粒径分布状态的新方法。

为达到以上目的，本发明采用如下技术方案：

基于计算流体力学流化床内颗粒密度和粒径分布预测方法，包括以下步骤：

步骤一、流化床内基本流动反应模型的建立

基于欧拉-欧拉方法，将气相和颗粒相均看作连续相，采用颗粒动力学理论描述颗粒相性质，颗粒相可以为一种或多种，根据实际情况确定，

气相连续性方程：

颗粒相连续性方程：

${\mathrm{\α}}_g+\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}=1---\left(3\right)$

其中下标g表示气相，下标pi表示第i种颗粒相；α为体积分数，ρ为密度，v为速度，S_m为非均相反应导致的质量变化源项。

气相动量方程：

颗粒相动量方程：

其中p为压力，τ为黏性应力张量，g为重力加速度，β为气相和颗粒相的曳力系数，ζ为不同颗粒相间的曳力系数，S_v为非均相反应导致的动量变化源项；

气相能量方程：

颗粒相能量方程：

其中H为焓值，λ为导热系数，T为温度，h_gpi为气相和颗粒相之间的对流传热系数，由于在流化床反应器中新加入的物料通常仅占炉内床料总量的5％以下，因此忽略颗粒相和颗粒相之间的热量传递，S_h为由于非均相反应引起的能量变化源项：

$h_{\mathrm{gpi}}=\frac{6{\mathrm{\λ}}_g{\mathrm{\α}}_g{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{Nu}}_{\mathrm{pi}}}{d_{\mathrm{pi}}^2}---\left(8\right)$

其中Nu_pi为无因次Nusselt准数，对于气固两相流体系采用Gunn的经验关联式计算；

${\mathrm{Nu}}_{\mathrm{pi}}=(7-10{\mathrm{\α}}_g+{5\mathrm{\α}}_g^2)(1+0.7{\mathrm{Re}}_{\mathrm{pi}}^{0.2}{\mathrm{Pr}}_g^{1/3})+(1.33-2.4{\mathrm{\α}}_g+1.2{\mathrm{\α}}_g^2){\mathrm{Re}}_{\mathrm{pi}}^{0.7}{\mathrm{Pr}}_g^{1/3}---\left(9\right)$

${\mathrm{Pr}}_g=\frac{c_{p_g}{\mathrm{\μ}}_g}{{\mathrm{\λ}}_g}---\left(10\right)$

其中c_p为比热容，μ为粘度；

气相组分输送方程：

颗粒相组分输送方程：

其中Y为组分i的质量分数，D为扩散系数，R为均相反应速率，S_r为非均相反应速率；

颗粒拟温度方程：

其中Θ为颗粒拟温度；

气相剪应力：

颗粒相剪应力：

其中I为单位张量；

颗粒相压力：

$p_{\mathrm{pi}}=[1+2\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{pi}}+d_{\mathrm{pk}}}{{2d}_{\mathrm{pi}}}\right)}^3(1+e_{\mathrm{ik}}){\mathrm{\α}}_{\mathrm{pk}}{g}_{\mathrm{ik}}]{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\Θ}}_{\mathrm{pi}}---\left(17\right)$

其中e为碰撞恢复系数；

径向分布函数：

$g_{\mathrm{ik}}=\frac{d_{\mathrm{pi}}{g}_{\mathrm{pk}}+d_{\mathrm{pk}}{g}_{\mathrm{pi}}}{d_{\mathrm{pi}}+d_{\mathrm{pk}}}---\left(18\right)$

$g_{\mathrm{pi}}=\frac{d_{\mathrm{pi}}}{2}\underset{k=1}{\overset{2}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_{\mathrm{pk}}}{d_{\mathrm{pk}}}+{[1-{\left(\frac{{\mathrm{\α}}_p}{{\mathrm{\α}}_{p,\max }}\right)}^{\frac{1}{3}}]}^{-1}---\left(19\right)$

${\mathrm{\α}}_p=\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}---\left(20\right)$

对于包含多种颗粒相的体系，最大堆积极限并不是一个固定的值，以下为计算颗粒最大堆积极限的关联式：

当 $X_i\≤\frac{{\mathrm{\α}}_{\mathrm{pi},\max }}{{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }}$

$\begin{array}{c}{\mathrm{\α}}_{p,\max }=[{\mathrm{\α}}_{\mathrm{pi},\max }-{\mathrm{\α}}_{\mathrm{pk},\max }+(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\\ \×[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\frac{X_i}{{\mathrm{\α}}_{\mathrm{pi},\max }}+{\mathrm{\α}}_{\mathrm{pk},\max }\end{array}---\left(21\right)$

否则

${\mathrm{\α}}_{p,\max }=(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }](1-X_i)+{\mathrm{\α}}_{\mathrm{pi},\max }---\left(22\right)$

其中

$X_i=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_p}---\left(23\right)$

颗粒体相黏度：

${\mathrm{\λ}}_{\mathrm{pi}}=\frac{4}{3}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(24\right)$

颗粒相剪切黏度由碰撞项、动力项及摩擦项三项组成：

${\mathrm{\μ}}_{\mathrm{pi}}=\frac{4}{5}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}+\frac{10{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{96{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{4}{5}{g}_{\mathrm{ik}}{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})]}^2+\frac{p_{\mathrm{pi}}\sin {\mathrm{\θ}}_i}{2\sqrt{I_{2D}}}---\left(25\right)$

其中I_2D为偏应力张量的第二不变量；

颗粒拟温度的输送系数：

$k_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{150{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{384(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{6}{5}{\mathrm{\α}}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})]}^2+2{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{d}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(26\right)$

颗粒碰撞导致的颗粒拟温度耗散项：

${\mathrm{\γ}}_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{12(1-e_{\mathrm{ik}}^2)g_{\mathrm{ik}}}{d_{\mathrm{pi}}\sqrt{\mathrm{\π}}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{\mathrm{\Θ}}_{\mathrm{pi}}^{3/2}---\left(27\right)$

气相与颗粒相之间的曳力系数采用Gidaspow曳力模型：

当α_g>0.8，

其中C_D为单颗粒曳力系数，可由下式获得：

$C_D=\left\{\begin{array}{c}\frac{24}{{\mathrm{Re}}_{\mathrm{pi}}}[1+0.15{\left({\mathrm{Re}}_{\mathrm{pi}}\right)}^{0.687}],{\mathrm{Re}}_{\mathrm{pi}}\≤1000\\ 0.44,{\mathrm{Re}}_{\mathrm{pi}}>1000\end{array}\right.---\left(29\right)$

当α_g≤0.8，

颗粒相与颗粒相之间的曳力系数采用Syamlal曳力模型：

其中C_fr为颗粒间的摩擦系数；

均相和非均相化学反应速率均可采用阿伦尼乌斯方程描述：

$R_i=\mathrm{kexp}(-\frac{E}{\mathrm{RT}})\mathrm{\Π}{\left(C_i\right)}^n---\left(32\right)$

其中R_i为i组分的化学反应速率，E为反应活化能，R为理想气体常数，C_i为i组分摩尔浓度，n为反应级数；

步骤二、建立描述颗粒相密度和粒径变化规律的数学模型

当颗粒在流化床内发生非均相化学反应时，通常由以下反应式表示：

A(s)+B(g)→C(s)+D(g)

颗粒中化合物A与气相中化合物B发生化学反应生成颗粒中化合物C和气相中化合物D，其中气相化合物B和D并不是必须存在的；颗粒中纯化合物A和C的密度分别为ρ_A和ρ_C，质量分数分别为Y_A和Y_C。化学反应式中化合物C和化合物A的质量比为a；

颗粒在流化床内发生非均相化学反应时，密度和粒径将会同时或单独发生变化，由颗粒性质和化学反应类型所决定；而颗粒的密度和粒径变化规律与所发生的化学反应以及颗粒中各化合物的质量分数和密度相关，并且必须相互匹配，遵循质量守恒定律；

对于尚未发生化学反应的颗粒，颗粒中仅含化合物A，假设其质量为m_A0，密度为ρ_A，粒径为d_p0，则此时颗粒体积为：

$V_0=\frac{m_{A0}}{{\mathrm{\ρ}}_A}---\left(33\right)$

假设m_A化合物A发生反应生成m_C化合物C，颗粒体积为V，粒径为d_p。

对于密度和粒径同时发生变化的颗粒，其颗粒相密度变化数学模型：

$\mathrm{\ρ}=\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C}}---\left(34\right)$

根据化学反应式和质量守恒定律及方程33可得：

m_C＝am_A  (35)

$Y_A=\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}---\left(36\right)$

$Y_C=1-Y_A=\frac{m_C}{m_{A0}-m_A+m_C}---\left(37\right)$

$V=\frac{m_{A0}-m_A+m_C}{\mathrm{\ρ}}=(m_{A0}-m_A+m_C)(\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C})=\frac{m_{A0}-m_A}{{\mathrm{\ρ}}_A}+\frac{m_C}{{\mathrm{\ρ}}_c}=V_0-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_A{\mathrm{\ρ}}_C}{m}_A---\left(38\right)$

由方程33、35和36可得：

Y_Am_A0-Y_Am_A+Y_Aam_A＝m_A0-m_A  (39)

m_A(1-Y_A+Y_Aa)＝m_A0(1-Y_A)  (40)

$m_A=\frac{m_{A0}(1-Y_A)}{1-Y_A(1-a)}=\frac{{\mathrm{\ρ}}_A{V}_0(1-Y_A)}{1-Y_A(1-a)}---\left(41\right)$

由方程38和41可得：

$V=V_0(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})---\left(42\right)$

则颗粒相粒径变化数学模型：

$d_p={(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(\text{43}\right)$

对于仅粒径发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相密度数学模型：

ρ＝ρ_A＝ρ_C  (44)

颗粒相粒径变化数学模型可由方程43和44获得：

$d_p={(1-\frac{(1-a)(1-Y_A)}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(\text{45}\right)$

对于仅密度发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相粒径数学模型：

V＝V₀  (46)

d_p＝d_p0  (47)

结合方程33、35、36和37可推导出颗粒相密度变化数学模型：

$\begin{array}{c}\mathrm{\ρ}=\frac{m_{A0}-m_A+m_C}{V}=\frac{m_{A0}-m_A+m_C}{V_0}=\frac{m_{A0}-m_A+m_C}{\frac{m_{A0}}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{m_{A0}-m_A+m_A}{{\mathrm{\ρ}}_A}\frac{1}{m_{A0}-m_A+m_C}}\\ =\frac{1}{\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}+\frac{m_C/a}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{a{\mathrm{\ρ}}_A}}\end{array}---\left(48\right)$

为了与方程34的形式一致，则可得：

ρ_C＝aρ_A  (49)

步骤三、流化床内密度和粒径分布状态的预测

基于流动反应模型和颗粒相密度和粒径变化数学模型对流化床内密度和和粒径分布状态进行模拟预测，首先根据流化床结构对模拟体系进行网格划分，设置基本流动反应模型，输入各化合物、颗粒的物性及反应动力学数据，定义出入口和壁面边界条件，设置时间步长和收敛条件开始求解；求解时根据方程1、2、4、5和14先求解连续性、动量和颗粒拟温度方程，随后根据方程6和7求解能量方程，再根据方程12和13求解组分方程，根据各组分含量和相应的颗粒密度和粒径变化数学模型，如密度和粒径同时发生变化根据方程34和43；如仅粒径发生变化根据方程44和45；如仅密度发生变化根据方程47、48和49对颗粒相密度和粒径进行修正更新。若整个计算体系内连续性方程、动量方程、颗粒拟温度方程、能量方程和组分方程两侧差值的绝对值之和(即为残差)均小于0.001，则计算结果收敛，否则不收敛；如不收敛则重复迭代，如收敛则判断计算时间是否完成，如未完成刚进入下一时间步进行求解，如完成则停止计算求解，导出颗粒密度和粒径数据，获得流化床内颗粒密度和粒径的分布状态。

本发明和现有技术相比，具有如下优点：

本发明所述基于计算流体力学的流化床内颗粒密度和粒径分布预测方法可广泛用于获得流化床生物质、煤等热解、燃烧或气化等过程中颗粒密度和粒径分布状态。该方法通过计算流体力学方法获得流化床内颗粒密度和粒径分布状态，无需对实际运行流化床进行复杂的采样分析，可以节约大量的人力、物力及时间成本。

本发明的有益效果是：

本发明采用计算流体力学方法对流化床进行模拟研究，并结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，从而准确预测流化床内颗粒密度和粒径的分布状态。该方法无需对实际运行流化床进行复杂的采样分析，可以节约大量的人力、物力及时间成本，为流化床的性能预测、优化控制及设计放大提供理论基础。

附图说明

图1生物质快速热解流化床。

图2为实施例一颗粒密度和粒径同时变化工况示意图。

图3为实施例二仅颗粒粒径发生变化工况示意图。

图4为实施例三仅颗粒密度发生变化工况示意图。

图5本发明步骤三流程图。

具体实施方式

下面结合附图对本发明做详细叙述。

基于计算流体力学流化床内颗粒密度和粒径分布预测方法，包括以下步骤：

步骤一、流化床内基本流动反应模型的建立

基于欧拉-欧拉方法，将气相和颗粒相均看作连续相，采用颗粒动力学理论描述颗粒相性质，颗粒相可以为一种或多种，根据实际情况确定，

气相连续性方程：

颗粒相连续性方程：

${\mathrm{\α}}_g+\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}=1---\left(52\right)$

其中下标g表示气相，下标pi表示第i种颗粒相；α为体积分数，ρ为密度，v为速度，S_m为非均相反应导致的质量变化源项。

气相动量方程：

颗粒相动量方程：

其中p为压力，τ为黏性应力张量，g为重力加速度，β为气相和颗粒相的曳力系数，ζ为不同颗粒相间的曳力系数，S_v为非均相反应导致的动量变化源项。

气相能量方程：

颗粒相能量方程：

其中H为焓值，λ为导热系数，T为温度，h_gpi为气相和颗粒相之间的对流传热系数。由于在流化床反应器中新加入的物料通常仅占炉内床料总量的5％以下，因此忽略颗粒相和颗粒相之间的热量传递。S_h为由于非均相反应引起的能量变化源项。

$h_{\mathrm{gpi}}=\frac{6{\mathrm{\λ}}_g{\mathrm{\α}}_g{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{Nu}}_{\mathrm{pi}}}{d_{\mathrm{pi}}^2}---\left(57\right)$

其中Nu_pi为无因次Nusselt准数，对于气固两相流体系采用Gunn的经验关联式计算。

${\mathrm{Nu}}_{\mathrm{pi}}=(7-10{\mathrm{\α}}_g+{5\mathrm{\α}}_g^2)(1+0.7{\mathrm{Re}}_{\mathrm{pi}}^{0.2}{\mathrm{Pr}}_g^{1/3})+(1.33-2.4{\mathrm{\α}}_g+1.2{\mathrm{\α}}_g^2){\mathrm{Re}}_{\mathrm{pi}}^{0.7}{\mathrm{Pr}}_g^{1/3}---\left(58\right)$

${\mathrm{Pr}}_g=\frac{c_{p_g}{\mathrm{\μ}}_g}{{\mathrm{\λ}}_g}---\left(59\right)$

其中c_p为比热容，μ为粘度。

气相组分输送方程：

颗粒相组分输送方程：

其中Y为组分i的质量分数，D为扩散系数，R为均相反应速率，S_r为非均相反应速率。

颗粒拟温度方程：

其中Θ为颗粒拟温度。

气相剪应力：

颗粒相剪应力：

其中I为单位张量。

颗粒相压力：

$p_{\mathrm{pi}}=[1+2\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{pi}}+d_{\mathrm{pk}}}{{2d}_{\mathrm{pi}}}\right)}^3(1+e_{\mathrm{ik}}){\mathrm{\α}}_{\mathrm{pk}}{g}_{\mathrm{ik}}]{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\Θ}}_{\mathrm{pi}}---\left(66\right)$

其中e为碰撞恢复系数。

径向分布函数：

$g_{\mathrm{ik}}=\frac{d_{\mathrm{pi}}{g}_{\mathrm{pk}}+d_{\mathrm{pk}}{g}_{\mathrm{pi}}}{d_{\mathrm{pi}}+d_{\mathrm{pk}}}---\left(67\right)$

$g_{\mathrm{pi}}=\frac{d_{\mathrm{pi}}}{2}\underset{k=1}{\overset{2}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_{\mathrm{pk}}}{d_{\mathrm{pk}}}+{[1-{\left(\frac{{\mathrm{\α}}_p}{{\mathrm{\α}}_{p,\max }}\right)}^{\frac{1}{3}}]}^{-1}---\left(68\right)$

${\mathrm{\α}}_p=\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}---\left(69\right)$

对于包含多种颗粒相的体系，最大堆积极限并不是一个固定的值，以下为计算颗粒最大堆积极限的关联式：

当 $X_i\≤\frac{{\mathrm{\α}}_{\mathrm{pi},\max }}{{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }}$

$\begin{array}{c}{\mathrm{\α}}_{p,\max }=[{\mathrm{\α}}_{\mathrm{pi},\max }-{\mathrm{\α}}_{\mathrm{pk},\max }+(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\\ \×[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\frac{X_i}{{\mathrm{\α}}_{\mathrm{pi},\max }}+{\mathrm{\α}}_{\mathrm{pk},\max }\end{array}---\left(70\right)$

否则

${\mathrm{\α}}_{p,\max }=(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }](1-X_i)+{\mathrm{\α}}_{\mathrm{pi},\max }---\left(71\right)$

其中

$X_i=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_p}---\left(72\right)$

颗粒体相黏度：

${\mathrm{\λ}}_{\mathrm{pi}}=\frac{4}{3}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(73\right)$

颗粒相剪切黏度由碰撞项、动力项及摩擦项三项组成：

${\mathrm{\μ}}_{\mathrm{pi}}=\frac{4}{5}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}+\frac{10{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{96{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{4}{5}{g}_{\mathrm{ik}}{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})]}^2+\frac{p_{\mathrm{pi}}\sin {\mathrm{\θ}}_i}{2\sqrt{I_{2D}}}---\left(74\right)$

其中I_2D为偏应力张量的第二不变量。

颗粒拟温度的输送系数：

$k_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{150{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{384(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{6}{5}{\mathrm{\α}}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})]}^2+2{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{d}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(75\right)$

颗粒碰撞导致的颗粒拟温度耗散项：

${\mathrm{\γ}}_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{12(1-e_{\mathrm{ik}}^2)g_{\mathrm{ik}}}{d_{\mathrm{pi}}\sqrt{\mathrm{\π}}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{\mathrm{\Θ}}_{\mathrm{pi}}^{3/2}---\left(76\right)$

气相与颗粒相之间的曳力系数采用Gidaspow曳力模型：

当α_g>0.8，

其中C_D为单颗粒曳力系数，可由下式获得：

$C_D=\left\{\begin{array}{c}\frac{24}{{\mathrm{Re}}_{\mathrm{pi}}}[1+0.15{\left({\mathrm{Re}}_{\mathrm{pi}}\right)}^{0.687}],{\mathrm{Re}}_{\mathrm{pi}}\≤1000\\ 0.44,{\mathrm{Re}}_{\mathrm{pi}}>1000\end{array}\right.---\left(78\right)$

当α_g≤0.8，

颗粒相与颗粒相之间的曳力系数采用Syamlal曳力模型：

其中C_fr为颗粒间的摩擦系数。

均相和非均相化学反应速率均可采用阿伦尼乌斯方程描述：

$R_i=\mathrm{kexp}(-\frac{E}{\mathrm{RT}})\mathrm{\Π}{\left(C_i\right)}^n---\left(81\right)$

其中R_i为i组分的化学反应速率，E为反应活化能，R为理想气体常数，C_i为i组分摩尔浓度，n为反应级数；

步骤二、建立描述颗粒相密度和粒径变化规律的数学模型

当颗粒在流化床内发生非均相化学反应时，通常由以下反应式表示：

A(s)+B(g)→C(s)+D(g)

颗粒中化合物A与气相中化合物B发生化学反应生成颗粒中化合物C和气相中化合物D，其中气相化合物B和D并不是必须存在的。颗粒中纯化合物A和C的密度分别为ρ_A和ρ_C，质量分数分别为Y_A和Y_C。化学反应式中化合物C和化合物A的质量比为a。

颗粒在流化床内发生非均相化学反应时，密度和粒径将会同时或单独发生变化，由颗粒性质和化学反应类型所决定。而颗粒的密度和粒径变化规律与所发生的化学反应以及颗粒中各化合物的质量分数和密度相关，并且必须相互匹配，遵循质量守恒定律。

对于尚未发生化学反应的颗粒，颗粒中仅含化合物A，假设其质量为m_A0，密度为ρ_A，粒径为d_p0，则此时颗粒体积为：

$V_0=\frac{m_{A0}}{{\mathrm{\ρ}}_A}---\left(82\right)$

假设m_A化合物A发生反应生成m_C化合物C，颗粒体积为V，粒径为d_p。

对于密度和粒径同时发生变化的颗粒，其颗粒相密度变化数学模型：

$\mathrm{\ρ}=\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C}}---\left(83\right)$

根据化学反应式和质量守恒定律及方程33可得：

m_C＝am_A  (84)

$Y_A=\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}---\left(85\right)$

$Y_C=1-Y_A=\frac{m_C}{m_{A0}-m_A+m_C}---\left(86\right)$

$V=\frac{m_{A0}-m_A+m_C}{\mathrm{\ρ}}=(m_{A0}-m_A+m_C)(\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C})=\frac{m_{A0}-m_A}{{\mathrm{\ρ}}_A}+\frac{m_C}{{\mathrm{\ρ}}_c}=V_0-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_A{\mathrm{\ρ}}_C}{m}_A---\left(87\right)$

由方程33、35和36可得：

Y_Am_A0-Y_Am_A+Y_Aam_A＝m_A0-m_A  (88)

m_A(1-Y_A+Y_Aa)＝m_A0(1-Y_A)  (89)

$m_A=\frac{m_{A0}(1-Y_A)}{1-Y_A(1-a)}=\frac{{\mathrm{\ρ}}_A{V}_0(1-Y_A)}{1-Y_A(1-a)}---\left(90\right)$

由方程38和41可得：

$V=V_0(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})---\left(91\right)$

则颗粒相粒径变化数学模型：

$d_p={(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(\text{92}\right)$

对于仅粒径发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相密度数学模型：

ρ＝ρ_A＝ρ_C  (93)

颗粒相粒径变化数学模型可由方程43和44获得：

$d_p={(1-\frac{(1-a)(1-Y_A)}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(\text{94}\right)$

对于仅密度发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相粒径数学模型：

V＝V₀  (95)

d_p＝d_p0  (96)

结合方程33、35、36和37可推导出颗粒相密度变化数学模型：

$\begin{array}{c}\mathrm{\ρ}=\frac{m_{A0}-m_A+m_C}{V}=\frac{m_{A0}-m_A+m_C}{V_0}=\frac{m_{A0}-m_A+m_C}{\frac{m_{A0}}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{m_{A0}-m_A+m_A}{{\mathrm{\ρ}}_A}\frac{1}{m_{A0}-m_A+m_C}}\\ =\frac{1}{\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}+\frac{m_C/a}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{a{\mathrm{\ρ}}_A}}\end{array}---\left(97\right)$

为了与方程34的形式一致，则可得：

ρ_C＝aρ_A  (98)

参照图5，步骤三、流化床内密度和粒径分布状态的预测

基于流动反应模型和颗粒相密度和粒径变化数学模型对流化床内密度和和粒径分布状态进行模拟预测，首先根据流化床结构对模拟体系进行网格划分，设置基本流动反应模型，输入各化合物、颗粒的物性及反应动力学数据，定义出入口和壁面边界条件，设置时间步长和收敛条件开始求解；求解时根据方程1、2、4、5和14先求解连续性、动量和颗粒拟温度方程，随后根据方程6和7求解能量方程，再根据方程12和13求解组分方程，根据各组分含量和相应的颗粒密度和粒径变化数学模型，如密度和粒径同时发生变化根据方程34和43；如仅粒径发生变化根据方程44和45；如仅密度发生变化根据方程47、48和49对颗粒相密度和粒径进行修正更新。若整个计算体系内连续性方程、动量方程、颗粒拟温度方程、能量方程和组分方程两侧差值的绝对值之和(即为残差)均小于0.001，则计算结果收敛，否则不收敛；如不收敛则重复迭代，如收敛则判断计算时间是否完成，如未完成刚进入下一时间步进行求解，如完成则停止计算求解，导出颗粒密度和粒径数据，获得流化床内颗粒密度和粒径的分布状态。

以生物质流化床快速热解过程为具体实施例对本发明作进一步的详细描述。

生物质在小型流化床内进行快速热解反应。流化床内径为3.81cm，高度为34.29cm。氮气从流化床底部进入，而生物质在高度为1.7cm处与伴随氮气一起进入体系。底部氮气温度为773K，速度为0.319m/s。生物质流率为100g/h，伴随氮气气速为0.143m/s，温度均为300K。反应温度由加热恒温装置维持在773K。初始床层内填充有5.5cm高的砂子，密度为2649kg/m³，粒径为520μm。生物质的快速热解化学反应可由下式表示：

A(s)→C(s)+D(g)

其中A为表示生物质的化合物，密度ρ_A为400kg/m³，C为表示焦炭的化合物，D为热解气体产物。该化学反应式中化合物C和化合物A的质量比a为0.27。颗粒未发生化学反应时的d_p0初始粒径为325μm。热解反应速率常数为1.3×10¹⁰，活化能为1.505×10⁸J/kmol。

在商业软件Fluent中对生物质流化床快速热解过程建立流动反应数学模型，将该流化床简化为二维结构如图1所示，网格尺寸为3.81mm宽×3.65mm高，根据上述内容设置基础物性、反应动力学数据及出入口条件，气相和颗粒相的壁面条件均为无滑移壁面条件，忽略湍流作用，时间步长为0.001s，收敛条件为残差小于0.001。通过用户自定义函数根据颗粒密度和粒径变化数学模型对颗粒密度和粒径进行更新修正，模拟计算120s。以下为分别结合三种不同颗粒密度和粒径变化规律所获得的流化床内颗粒密度和粒径分布状态。

实施例一

如生物质快速热解过程中密度和粒径同时发生变化，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相密度变化数学模型：

$\mathrm{\ρ}=\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C}}$

颗粒相粒径变化数学模型：

$d_p={(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})}^{1/3}{d}_{p0}$

设置焦炭密度ρ_C为2330kg/m³，在数值模拟过程中根据颗粒相密度和粒径变化数学模型对颗粒密度和粒径进行实时修正，从而获得流化床内颗粒密度和粒径分布状态。计算120s时颗粒的体积分数、密度和粒径分布云图如图2所示。从图中可以看出颗粒主要聚集在床层的右上部，其密度和粒径在流化床内均发生了变化。颗粒密度分布范围为400-1800kg/m³，而颗粒粒径分布范围为125-325μm。流化床内大部分颗粒的密度约为1470kg/m³，粒径约为140μm。

实施例二

如生物质快速热解过程中仅粒径发生变化，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相密度数学模型：

ρ＝ρ_A＝ρ_C

颗粒相粒径变化数学模型：

$d_p={(1-\frac{(1-a)(1-Y_A)}{1-Y_A(1-a)})}^{1/3}{d}_{p0}$

设置焦炭密度ρ_C与生物质密度一致ρ_A，均为400kg/m³，在数值模拟过程中根据颗粒相粒径变化数学模型对粒径进行实时修正，从而获得流化床内粒径分布状态。计算120s时颗粒的体积分数和粒径分布云图如图3所示。从图中可以看出流化床内颗粒的粒径分布范围为225-325μm，大部分颗粒的粒径约为260μm。

实施例三

如生物质快速热解过程中仅密度发生变化，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相粒径数学模型：

d_p＝d_p0

颗粒相密度变化数学模型：

ρ_C＝aρ_A

$\mathrm{\ρ}=\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C}}$

保持颗粒相粒径不变为325μm，设置焦炭密度ρ_C＝aρ_A＝108kg/m³，在数值模拟过程中根据颗粒相密度变化数学模型对密度进行实时修正，从而获得流化床内密度分布状态。计算120s时颗粒的体积分数和密度分布云图如图4所示。从图中可以看出流化床内颗粒的密度分布范围为170-400kg/m³，大部分颗粒密度约为300kg/m³。

由以上三个实施例可知，本发明所阐述计算流体力学方法可以根据不同颗粒密度和变化规律分别获得相应的颗粒密度和粒径分布状态，为流化床的性能预测、优化控制及设计放大提供理论基础。

Claims (1)

1.基于计算流体力学流化床内颗粒密度和粒径分布预测方法，其特征在于，包括以下步骤：

步骤一、流化床内基本流动反应模型的建立

基于欧拉-欧拉方法，将气相和颗粒相均看作连续相，采用颗粒动力学理论描述颗粒相性质，颗粒相可以为一种或多种，根据实际情况确定，

气相连续性方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\right)+\▿\·\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{v_g}\right)=S_{\mathrm{mg}}---\left(1\right)$

颗粒相连续性方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\right)+\▿\·\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}\right)=S_{\mathrm{mpi}}---\left(2\right)$

${\mathrm{\α}}_g+\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}=1---\left(3\right)$

其中下标g表示气相，下标pi表示第i种颗粒相；α为体积分数，ρ为密度，v为速度，S_m为非均相反应导致的质量变化源项；

气相动量方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{v_g}\right)+\▿\·\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{v_g}\stackrel{\→}{v_g}\right)=-{\mathrm{\α}}_g\▿p+\▿\·\stackrel{\‾}{\stackrel{\‾}{{\mathrm{\τ}}_g}}+{\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{g}+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\β}}_i(\stackrel{\→}{v_{\mathrm{pi}}}-\stackrel{\→}{v_g})+S_{\mathrm{vg}}---\left(4\right)$

颗粒相动量方程：

$\begin{array}{c}\frac{\∂}{\∂t}\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}\right)+\▿\·\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}\stackrel{\→}{v_{\mathrm{pi}}}\right)=-{\mathrm{\α}}_p\▿p+\▿\·\stackrel{\‾}{\stackrel{\‾}{{\mathrm{\τ}}_{\mathrm{pi}}}}+{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{g}+{\mathrm{\β}}_i(\stackrel{\→}{v_g}-\stackrel{\→}{v_{\mathrm{pi}}})\\ +{\mathrm{\ζ}}_{\mathrm{ik}}(\stackrel{\→}{v_{\mathrm{pk}}}-\stackrel{\→}{v_{\mathrm{pi}}})+S_{\mathrm{vpi}}\end{array}---\left(5\right)$

其中p为压力，τ为黏性应力张量，g为重力加速度，β为气相和颗粒相的曳力系数，ζ为不同颗粒相间的曳力系数，S_v为非均相反应导致的动量变化源项；

气相能量方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g{H}_g\right)+\▿\·\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{v_g}{H}_g\right)=\▿\·({\mathrm{\λ}}_g\▿T_g)+\underset{i}{\mathrm{\Σ}}\left[h_{\mathrm{gpi}}(T_{\mathrm{pi}}-T_g)\right]+S_{\mathrm{hg}}---\left(6\right)$

颗粒相能量方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{H}_{\mathrm{pi}}\right)+\▿\·\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}{H}_{\mathrm{pi}}\right)=\▿\·({\mathrm{\λ}}_{\mathrm{pi}}\▿T_{\mathrm{pi}})+h_{\mathrm{gpi}}(T_g-T_{\mathrm{pi}})+S_{\mathrm{hpi}}---\left(7\right)$

其中H为焓值，λ为导热系数，T为温度，h_gpi为气相和颗粒相之间的对流传热系数，由于在流化床反应器中新加入的物料通常仅占炉内床料总量的5％以下，因此忽略颗粒相和颗粒相之间的热量传递，S_h为由于非均相反应引起的能量变化源项：

$h_{\mathrm{gpi}}=\frac{6{\mathrm{\λ}}_g{\mathrm{\α}}_g{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{Nu}}_{\mathrm{pi}}}{d_{\mathrm{pi}}^2}---\left(8\right)$

其中Nu_pi为无因次Nusselt准数，对于气固两相流体系采用Gunn的经验关联式计算；

${\mathrm{Nu}}_{\mathrm{pi}}=(7-10{\mathrm{\α}}_g+5{\mathrm{\α}}_g^2)(1+0.7R{e}_{\mathrm{pi}}^{0.2}{\mathrm{Pr}}_g^{1/3})+(1.33-2.4{\mathrm{\α}}_g+1.2{\mathrm{\α}}_g^2){\mathrm{Re}}_{\mathrm{pi}}^{0.7}{\mathrm{Pr}}_g^{1/3}---\left(9\right)$

${\mathrm{Pr}}_g=\frac{c_{p_g}{\mathrm{\μ}}_g}{{\mathrm{\λ}}_g}---\left(10\right)$

${\mathrm{Re}}_{\mathrm{pi}}=\frac{{\mathrm{\α}}_g{\mathrm{\ρ}}_g|\stackrel{\→}{v_g}-\stackrel{\→}{v_{\mathrm{pi}}}|d_{\mathrm{pi}}}{{\mathrm{\μ}}_g}---\left(11\right)$

其中c_p为比热容，μ为粘度；

气相组分输送方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g{Y}_{\mathrm{gi}}\right)+\▿\·\left({\mathrm{\α}}_g{\mathrm{\ρ}}_g\stackrel{\→}{v_g}{Y}_{\mathrm{gi}}\right)=\▿\·({\mathrm{\α}}_g{\mathrm{\ρ}}_g{D}_{i,m}\▿Y_i]+R_{\mathrm{gi}}+S_{\mathrm{rgi}}---\left(12\right)$

颗粒相组分输送方程：

$\frac{\∂}{\∂t}\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{Y}_{\mathrm{pi},j}\right)+\▿\·\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}{Y}_{\mathrm{pi},j}\right)=\▿\·({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{D}_{i,m}\▿Y_{\mathrm{pi},j}]+R_{\mathrm{pi},j}+S_{\mathrm{rpi},j}---\left(13\right)$

其中Y为组分i的质量分数，D为扩散系数，R为均相反应速率，S_r为非均相反应速率；

颗粒拟温度方程：

$\frac{3}{2}[\frac{\∂}{\∂t}\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\Θ}}_{\mathrm{pi}}\right)+\▿\·\left({\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}\stackrel{\→}{v_{\mathrm{pi}}}{\mathrm{\Θ}}_{\mathrm{pi}}\right)]=\stackrel{\‾}{\stackrel{\‾}{{\mathrm{\τ}}_{\mathrm{pi}}}}:\▿\stackrel{\→}{v_{\mathrm{pi}}}+\▿\·(k_{{\mathrm{\Θ}}_{\mathrm{pi}}}\▿{\mathrm{\Θ}}_{\mathrm{pi}})-{\mathrm{\γ}}_{{\mathrm{\Θ}}_{\mathrm{pi}}}-3({\mathrm{\β}}_i+{\mathrm{\ζ}}_{\mathrm{ik}}){\mathrm{\Θ}}_{\mathrm{pi}}---\left(14\right)$

其中Θ为颗粒拟温度。

气相剪应力：

$\stackrel{\‾}{\stackrel{\‾}{{\mathrm{\τ}}_g}}={\mathrm{\α}}_g{\mathrm{\μ}}_g[\▿\stackrel{\→}{v_g}+{(\▿\stackrel{\→}{v_g})}^T-\frac{2}{3}(\▿\·\stackrel{\→}{v_g})\stackrel{\‾}{\stackrel{\‾}{I}}]---\left(15\right)$

颗粒相剪应力：

$\stackrel{\‾}{\stackrel{\‾}{{\mathrm{\τ}}_{\mathrm{pi}}}}=(-p_{\mathrm{pi}}+{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\λ}}_{\mathrm{pi}}\▿\·\stackrel{\→}{v_{\mathrm{pi}}})\stackrel{\‾}{\stackrel{\‾}{I}}+{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\μ}}_{\mathrm{pi}}[\▿\stackrel{\→}{v_{\mathrm{pi}}}+{(\▿\stackrel{\→}{v_{\mathrm{pi}}})}^T-\frac{2}{3}(\▿\·\stackrel{\‾}{v_{\mathrm{pi}}})\stackrel{\‾}{\stackrel{\‾}{I}}]---\left(16\right)$

其中I为单位张量；

颗粒相压力：

$p_{\mathrm{pi}}=[1+2\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\left(\frac{d_{\mathrm{pi}}+d_{\mathrm{pk}}}{2d_{\mathrm{pi}}}\right)}^3(1+e_{\mathrm{ik}}){\mathrm{\α}}_{\mathrm{pk}}{g}_{\mathrm{ik}}]{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\Θ}}_{\mathrm{pi}}---\left(17\right)$

其中e为碰撞恢复系数；

径向分布函数：

$g_{\mathrm{ik}}=\frac{d_{\mathrm{pi}}{g}_{\mathrm{pk}}+d_{\mathrm{pk}}{g}_{\mathrm{pi}}}{d_{\mathrm{pi}}+d_{\mathrm{pk}}}---\left(18\right)$

$g_{\mathrm{pi}}=\frac{d_{\mathrm{pi}}}{2}\underset{k=1}{\overset{2}{\mathrm{\Σ}}}\frac{{\mathrm{\α}}_{\mathrm{pk}}}{d_{\mathrm{pk}}}+{[1-{\left(\frac{{\mathrm{\α}}_p}{{\mathrm{\α}}_{p,\max }}\right)}^{\frac{1}{3}}]}^{-1}---\left(19\right)$

${\mathrm{\α}}_p=\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{pk}}---\left(20\right)$

对于包含多种颗粒相的体系，最大堆积极限并不是一个固定的值，以下为计算颗粒最大堆积极限的关联式：

当 $X_i\≤\frac{{\mathrm{\α}}_{\mathrm{pi},\max }}{{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }}$

$\begin{array}{c}{\mathrm{\α}}_{p,\max }=[{\mathrm{\α}}_{\mathrm{pi},\max }-{\mathrm{\α}}_{\mathrm{pk},\max }+(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\\ \×[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }]\frac{X_i}{{\mathrm{\α}}_{\mathrm{pi},\max }}+{\mathrm{\α}}_{\mathrm{pk},\max }\end{array}---\left(21\right)$

否则

${\mathrm{\α}}_{p,\max }=(1-\sqrt{\frac{d_{\mathrm{pk}}}{d_{\mathrm{pi}}}})[{\mathrm{\α}}_{\mathrm{pi},\max }+(1-{\mathrm{\α}}_{\mathrm{pi},\max }){\mathrm{\α}}_{\mathrm{pk},\max }](1-X_i)+{\mathrm{\α}}_{\mathrm{pi},\max }---\left(22\right)$

其中

$X_i=\frac{{\mathrm{\α}}_{\mathrm{pi}}}{{\mathrm{\α}}_p}---\left(23\right)$

颗粒体相黏度：

${\mathrm{\λ}}_{\mathrm{pi}}=\frac{4}{3}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(24\right)$

颗粒相剪切黏度由碰撞项、动力项及摩擦项三项组成：

${\mathrm{\μ}}_{\mathrm{pi}}=\frac{4}{5}{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}+\frac{10{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{96{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{4}{5}{g}_{\mathrm{ik}}{\mathrm{\α}}_{\mathrm{pi}}(1+e_{\mathrm{ik}})]}^2+\frac{p_{\mathrm{pi}}\sin {\mathrm{\θ}}_i}{2\sqrt{I_{2D}}}---\left(25\right)$

其中I_2D为偏应力张量的第二不变量；

颗粒拟温度的输送系数：

$k_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{150{\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}\sqrt{{\mathrm{\Θ}}_{\mathrm{pi}}\mathrm{\π}}}{384(1+e_{\mathrm{ik}})g_{\mathrm{ik}}}{[1+\frac{6}{5}{\mathrm{\α}}_{\mathrm{pi}}{g}_{\mathrm{ik}}(1+e_{\mathrm{ik}})]}^2+2{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{d}_{\mathrm{pi}}(1+e_{\mathrm{ik}})g_{\mathrm{ik}}\sqrt{\frac{{\mathrm{\Θ}}_{\mathrm{pi}}}{\mathrm{\π}}}---\left(26\right)$

颗粒碰撞导致的颗粒拟温度耗散项：

${\mathrm{\γ}}_{{\mathrm{\Θ}}_{\mathrm{pi}}}=\frac{12(1-e_{\mathrm{ik}}^2)g_{\mathrm{ik}}}{d_{\mathrm{pi}}\sqrt{\mathrm{\π}}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pi}}^2{\mathrm{\Θ}}_{\mathrm{pi}}^{3/2}---\left(27\right)$

气相与颗粒相之间的曳力系数采用Gidaspow曳力模型：

当α_g>0.8，

${\mathrm{\β}}_i=\frac{3}{4}{C}_D\frac{{\mathrm{\α}}_g{\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}|\stackrel{\→}{v_g}-\stackrel{\→}{v_{\mathrm{pi}}}|}{d_{\mathrm{pi}}}{{\mathrm{\α}}_g}^{-2.65}---\left(28\right)$

其中C_D为单颗粒曳力系数，可由下式获得：

$C_D=\left\{\begin{array}{cc}\frac{24}{{\mathrm{Re}}_{\mathrm{pi}}}[1+0.15{\left({\mathrm{Re}}_{\mathrm{pi}}\right)}^{0.687}],& {\mathrm{Re}}_{\mathrm{pi}}\≤1000\\ 0.44,& {\mathrm{Re}}_{\mathrm{pi}}>1000\end{array}\right.---\left(29\right)$

当α_g≤0.8，

${\mathrm{\β}}_i=150\frac{{\mathrm{\α}}_{\mathrm{pi}}(1-{\mathrm{\α}}_g){\mathrm{\μ}}_g}{{\mathrm{\α}}_g{d_{\mathrm{pi}}}^2}+1.75\frac{{\mathrm{\ρ}}_g{\mathrm{\α}}_{\mathrm{pi}}|\stackrel{\→}{v_g}-\stackrel{\→}{v_{\mathrm{pi}}}|}{d_{\mathrm{pi}}}---\left(30\right)$

颗粒相与颗粒相之间的曳力系数采用Syamlal曳力模型：

${\mathrm{\ζ}}_{\mathrm{ik}}=\frac{3(1+e_{\mathrm{ik}})(\frac{\mathrm{\π}}{2}+C_{\mathrm{fr},\mathrm{ik}}\frac{{\mathrm{\π}}^2}{8}){\mathrm{\α}}_{\mathrm{pi}}{\mathrm{\ρ}}_{\mathrm{pi}}{\mathrm{\α}}_{\mathrm{pk}}{\mathrm{\ρ}}_{\mathrm{pk}}{(d_{\mathrm{pi}}+d_{\mathrm{pk}})}^2{g}_{\mathrm{ik}}}{2\mathrm{\π}({\mathrm{\ρ}}_{\mathrm{pi}}{d}_{\mathrm{pi}}^3+{\mathrm{\ρ}}_{\mathrm{pk}}{d}_{\mathrm{pk}}^3)}|\stackrel{\→}{v_{\mathrm{pi}}}-\stackrel{\→}{v_{\mathrm{pk}}}|---\left(31\right)$

其中C_fr为颗粒间的摩擦系数；

均相和非均相化学反应速率均可采用阿伦尼乌斯方程描述：

$R_i=\mathrm{kexp}(-\frac{E}{\mathrm{RT}})\mathrm{\Π}{\left(C_i\right)}^n---\left(32\right)$

其中R_i为i组分的化学反应速率，E为反应活化能，R为理想气体常数，C_i为i组分摩尔浓度，n为反应级数；

步骤二、建立描述颗粒相密度和粒径变化规律的数学模型

当颗粒在流化床内发生非均相化学反应时，通常由以下反应式表示：

A(s)+B(g)→C(s)+D(g)

颗粒中化合物A与气相中化合物B发生化学反应生成颗粒中化合物C和气相中化合物D，其中气相化合物B和D并不是必须存在的；颗粒中纯化合物A和C的密度分别为ρ_A和ρ_C，质量分数分别为Y_A和Y_C；化学反应式中化合物C和化合物A的质量比为a；

颗粒在流化床内发生非均相化学反应时，密度和粒径将会同时或单独发生变化，由颗粒性质和化学反应类型所决定；而颗粒的密度和粒径变化规律与所发生的化学反应以及颗粒中各化合物的质量分数和密度相关，并且必须相互匹配，遵循质量守恒定律；

对于尚未发生化学反应的颗粒，颗粒中仅含化合物A，假设其质量为m_A0，密度为ρ_A，粒径为d_p0，则此时颗粒体积为：

$V_0=\frac{m_{A0}}{{\mathrm{\ρ}}_A}---\left(33\right)$

假设m_A化合物A发生反应生成m_C化合物C，颗粒体积为V，粒径为d_p。

对于密度和粒径同时发生变化的颗粒，其颗粒相密度变化数学模型：

$\mathrm{\ρ}=\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C}}---\left(34\right)$

根据化学反应式和质量守恒定律及方程33可得：

m_C＝am_A            (35)

$Y_A=\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}---\left(36\right)$

$Y_C=1-Y_A=\frac{m_C}{m_{A0}-m_A+m_C}---\left(37\right)$

$V=\frac{m_{A0}-m_A+m_C}{\mathrm{\ρ}}=(m_{A0}-m_A+m_C)(\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{{\mathrm{\ρ}}_C})=\frac{m_{A0}-m_A}{{\mathrm{\ρ}}_A}+\frac{m_C}{{\mathrm{\ρ}}_C}=V_0-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_A{\mathrm{\ρ}}_C}{m}_A---\left(38\right)$

由方程33、35和36可得：

Y_Am_A0-Y_Am_A+Y_Aam_A＝m_A0-m_A           (39)

m_A(1-Y_A+Y_Aa)＝m_A0(1-Y_A)             (40)

$m_A=\frac{m_{A0}(1-Y_A)}{1-Y_A(1-a)}=\frac{{\mathrm{\ρ}}_A{V}_0(1-Y_A)}{1-Y_A(1-a)}---\left(41\right)$

由方程38和41可得：

$V=V_0(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})---\left(42\right)$

则颗粒相粒径变化数学模型：

$d_p={(1-\frac{{\mathrm{\ρ}}_C-{\mathrm{\ρ}}_{A}a}{{\mathrm{\ρ}}_C}\frac{1-Y_A}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(43\right)$

对于仅粒径发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相密度数学模型：

ρ＝ρ_A＝ρ_C           (44)

颗粒相粒径变化数学模型可由方程43和44获得：

$d_p={(1-\frac{(1-a)(1-Y_A)}{1-Y_A(1-a)})}^{1/3}{d}_{p0}---\left(45\right)$

对于仅密度发生变化的颗粒，可用下列数学模型描述颗粒相粒径和密度变化规律：

颗粒相粒径数学模型：

V＝V₀             (46)

d_p＝d_p0             (47)

结合方程33、35、36和37可推导出颗粒相密度变化数学模型：

$\begin{array}{c}\mathrm{\ρ}=\frac{m_{A0}-m_A+m_C}{V}=\frac{m_{A0}-m_A+m_C}{V_0}=\frac{m_{A0}-m_A+m_C}{\frac{m_{A0}}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{m_{A0}-m_A+m_A}{{\mathrm{\ρ}}_A}\frac{1}{m_{A0}-m_A+m_C}}\\ =\frac{1}{\frac{m_{A0}-m_A}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}+\frac{m_C/a}{m_{A0}-m_A+m_C}\frac{1}{{\mathrm{\ρ}}_A}}\\ =\frac{1}{\frac{Y_A}{{\mathrm{\ρ}}_A}+\frac{Y_C}{a{\mathrm{\ρ}}_A}}\end{array}---\left(48\right)$

为了与方程34的形式一致，则可得：

ρ_C＝aρ_A              (49)

步骤三、流化床内密度和粒径分布状态的预测

基于流动反应模型和颗粒相密度和粒径变化数学模型对流化床内密度和和粒径分布状态进行模拟预测，首先根据流化床结构对模拟体系进行网格划分，设置基本流动反应模型，输入各化合物、颗粒的物性及反应动力学数据，定义出入口和壁面边界条件，设置时间步长和收敛条件开始求解；求解时根据方程1、2、4、5和14先求解连续性、动量和颗粒拟温度方程，随后根据方程6和7求解能量方程，再根据方程12和13求解组分方程，根据各组分含量和相应的颗粒密度和粒径变化数学模型，如密度和粒径同时发生变化根据方程34和43；如仅粒径发生变化根据方程44和45；如仅密度发生变化根据方程47、48和49对颗粒相密度和粒径进行修正更新；若整个计算体系内连续性方程、动量方程、颗粒拟温度方程、能量方程和组分方程两侧差值的绝对值之和(即为残差)均小于0.001，则计算结果收敛，否则不收敛；如不收敛则重复迭代，如收敛则判断计算时间是否完成，如未完成刚进入下一时间步进行求解，如完成则停止计算求解，导出颗粒密度和粒径数据，获得流化床内颗粒密度和粒径的分布状态。