CN107132156B - 一种颗粒密度和粒径动态变化流化床的模拟方法 - Google Patents

Abstract

一种颗粒密度和粒径动态变化流化床的模拟方法，步骤一、流化床内基本流动反应模型的建立；步骤二、建立描述颗粒相密度和粒径变化规律的数学模型；步骤三、建立颗粒类型分段曳力模型；步骤四、流化床内流动反应特性的预测；该方法采用计算流体力学方法对颗粒密度和粒径动态变化的流化床进行模拟研究，结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，通过颗粒类型分段曳力模型描述多种颗粒类型并存体系的气‑颗粒相间作用力，从而准确预测颗粒密度和粒径动态变化流化床内的流动反应特性。该方法无需进行大量复杂昂贵的实验研究，可以节约大量的人力、物力及时间成本。

Description

一种颗粒密度和粒径动态变化流化床的模拟方法

技术领域

本发明涉及流体力学领域，具体涉及一种基于计算流体力学预测颗粒密度和粒径动态变化流化床的模拟方法。

背景技术

流化床内的固体颗粒在流体的作用下具有流体某些表观特性，颗粒表面全部暴露于周围剧烈湍动的流体中，从而强化了传热、传质和化学反应，因此广泛用于工业领域如生物质、煤等固体燃料的热解、燃烧或气化等过程。流化床内生物质或煤等颗粒的密度和粒径随着热解、燃烧或气化等非均相化学反应的进行往往发生动态变化，表现出一定的密度和粒径变化规律，影响着流化床内的流动反应特性。计算流体力学是对研究对象建立流动反应模型并采用计算机和离散化的数值方法进行数值模拟和分析的一种方法。该方法不受现有实验技术的限制，可以全面、高效且低成本地揭示流化床内的流动反应特性，因此广泛应用于流化床的研究当中。

准确描述反应器内气相和颗粒的性质以及各种相互作用力是精确预测流化床反应器流动反应特性的重要基础。对于颗粒密度和粒径动态变化的流化床，在数值模拟时需要注意以下两点：

(1)颗粒的密度和粒径分布：目前大部分模拟工作均假设颗粒在非均相化学反应过程中粒径保持恒定不变，这与流化床反应器内的真实情况相差甚远，从而严重影响着模拟结果的准确性。虽然申请者之前的专利提出了一种基于计算流体力学流化床内颗粒密度和粒径分布预测方法(201510079554.X)，但是由于该发明仅能用于单组分一步反应机理，无法适用于更为复杂准确的多组分多步平行-顺序反应机理。

(2)颗粒的类型：Geldart依据颗粒的密度和粒径等性质将颗粒分成具有不同流态化特征的A、B、C及D类颗粒。在颗粒密度和粒径动态变化的流化床，极有可能同时存在着不同类型的颗粒。对于不同类型颗粒，其气-颗粒相间作用力(曳力)的计算方法存在较大的差异。而目前的模拟工作尚未考虑到颗粒类型可能发生的变化，对于指定的流化床反应器仅采用单一类型颗粒的曳力模型，难以准确描述多种颗粒类型并存体系的气-颗粒相间作用力。

因此，现有的模拟方法难以准确预测颗粒密度和粒径动态变化流化床内的流动反应特性，无法为此类反应器的精准调控、优化设计及过程放大提供理论基础。

发明内容

为了解决现有技术存在的问题，本发明的目的在于提供一种颗粒密度和粒径动态变化流化床的模拟方法，该方法采用计算流体力学方法对流化床进行模拟研究，结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，通过颗粒类型分段曳力模型描述多种颗粒类型并存体系的气-颗粒相间作用力，从而准确预测颗粒密度和粒径动态变化流化床内的流动反应特性。该方法不受现有实验技术的限制，可以全面、高效且低成本地揭示流化床内的流动反应特性，为此类反应器的精准调控、优化设计及过程放大提供理论基础。

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

一种颗粒密度和粒径动态变化流化床的模拟方法，包括以下步骤：

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

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

气相连续性方程：

颗粒相连续性方程：

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

气相动量方程：

颗粒相动量方程：

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

气相能量方程：

颗粒相能量方程：

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

起的能量变化源项：

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

其中c_p为比热容，μ为粘度，d为颗粒直径；

气相组分输送方程：

其中Y_gi为气相中组分i的质量分数，D_i,m为气相中组分i的扩散系数，R_gi为气相中组分i的均相反应速率，S_rgi为气相中组分i的非均相反应速率；

颗粒相组分输送方程：

其中Y_pi,j为第i种颗粒相中组分j的质量分数，D_j,m为第i种颗粒相中组分j的扩散系数，R_pi,j为第i种颗粒相中组分j的均相反应速率，S_rpi,j为第i种颗粒相中组分j的非均相反应速率；

颗粒拟温度方程：

其中Θ为颗粒拟温度；

气相剪应力：

颗粒相剪应力：

其中I为单位张量，σ为颗粒体相黏度；

颗粒相压力：

其中e为碰撞恢复系数；

径向分布函数：

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

当

否则

其中

颗粒体相黏度：

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

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

颗粒拟温度的输送系数：

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

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

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

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

R_i＝kΠ(C_i)ⁿ (29)

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

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

当颗粒在流化床内发生非均相化学反应时，通常多组分多步平行-顺序反应，可由以下反应式表示：

式中A、B、C、D、E、F、G、H、I、J、K、L为反应物和产物，下标n表示第n种组分，括号内s和g分别表示颗粒相组分和气相组分，k1_n、k2_n、k3_n、k4_n、k5_n为相应的速率常数；

对于密度和粒径动态变化的颗粒，其颗粒相密度变化数学模型可由下式表示：

式中Y_i和ρ_i分别为颗粒相中i组分的质量分数和密度；

而颗粒相粒径变化数学模型，则可根据颗粒尺度的质量守恒，由颗粒相的质量和密度变化数学模型获得，因此，需确定颗粒相的质量。

依据化学反应工程理论对反应式进行分析，对各颗粒相组分的生成\消耗速率建立一系列偏微分方程组，如下所示：

通过求解以上的偏微分方程组，即可获得颗粒相各组分的质量，从而可确定颗粒相的总质量和各组分的质量分数，最终根据综合密度变化数学模型方程(31)和颗粒尺度的质量守恒获得颗粒相粒径变化数学模型；

在求解偏微分方程时可采取两种方法，一种是求解析解，适用于相对较为简单的热解机理；另一种是求数值解，适用于极其复杂的热解机理。

(1)解析解法

该方法首先将各颗粒相组分的生成/消耗速率相比，即可消除反应时间项，获得各颗粒相组分之间的关系式。再根据颗粒相组分质量分数和密度变化数学模型，结合颗粒尺度的质量守恒分析获得粒径变化数学模型。

(2)数值解法

对于极其复杂的热解机理，由于组分数目多且反应步骤复杂，难以获得反应速率偏微分方程组的解析解。但是可运用数学软件对偏微分方程进行数值求解，获得不同反应时间下各颗粒相组分的质量，进而求得不同反应时间下颗粒质量及各颗粒相组分的质量分数，在数学软件中通过数据拟合的方法获得颗粒质量与各颗粒相组分质量分数的关联式，最后与密度变化数学模型结合进行颗粒尺度质量守恒分析获得粒径变化数学模型。

以上两种方法所建立粒径变化数学模型的通用形式为：

d_p＝f(Y_i,ρ_i)d_p0 (38)

式中d_p0为颗粒的初始粒径；

步骤三、建立颗粒类型分段曳力模型

对于不同类型的颗粒，特别是A类和B类颗粒，所采用的气-颗粒相曳力模型差异较大；因此，如果反应器内同时存在A类和B类颗粒，则需建立颗粒类型分段曳力模型，对不同类型的颗粒分别采用相应的曳力模型，并且需在一定程度上避免数值的跳跃：

(1)计算无因次量

(2)计算A类和B类颗粒的分界值：

(3)建立颗粒类型分段曳力模型

对于A类和B类颗粒分别选用相应的曳力模型，其曳力系数以β_A和β_B表示，则颗粒类型分段曳力模型为：

其中：

步骤四、流化床内流动反应特性的预测

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

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

本发明所述基于计算流体力学预测颗粒密度和粒径动态变化流化床流动反应特性的模拟方法，可广泛用于预测流化床生物质、煤等热解、燃烧或气化等反应器的流动反应特性。该方法通过计算流体力学方法获得流化床反应器的流动反应特性，无需进行大量复杂昂贵的实验研究，可以节约大量的人力、物力及时间成本。

本发明的有益效果是：

本发明采用计算流体力学方法对颗粒密度和粒径动态变化的流化床进行模拟研究，结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，通过颗粒类型分段曳力模型描述多种颗粒类型并存体系的气-颗粒相间作用力，从而准确预测颗粒密度和粒径动态变化流化床内的流动反应特性。该方法无需进行大量复杂昂贵的实验研究，可以节约大量的人力、物力及时间成本，为流化床的性能预测、优化控制及设计放大提供理论基础。

附图说明

图1是本发明步骤四流程图。

图2是数据拟合示意图。

图3是数值解与解析解两种方法获得的粒径变化数学模型对比图。

图4是流化床简化为二维结构的模拟体体系示意图。

图5是生物质的体积分数、密度和粒径分布模拟结果。

具体实施方式

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

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

其中b为生物质，密度ρ_b为660kg/m³，c为生物炭，密度ρ_c为350kg/m³，l为生物油，g为热解气。颗粒未发生化学反应时的d_p0初始粒径为600μm。各反应的反应动力学常数如表1所示。

表1反应动力学常数

反应	指前因子，s<sup>-1</sup>	活化能，J/mol	反应热，J/kg
				k<sub>1</sub>	2.66×10<sup>10</sup>	1.065×10<sup>5</sup>	4.18×10<sup>5</sup>
k<sub>2</sub>	1.48×10<sup>10</sup>	1.127×10<sup>5</sup>	4.18×10<sup>5</sup>
				k<sub>3</sub>	2.37×10<sup>10</sup>	8.4×10<sup>10</sup>	4.18×10<sup>5</sup>

根据本发明内容对该流化床内的流动反应特性进行模拟预测，具体步骤如下：

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

建立基本流动反应模型(方程1至30)。

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

建立颗粒相密度变化数学模型(方程31)。

建立颗粒相各组分反应速率偏微分方程组：

式中m为组分的质量；

(1)解析解法

两式相比获得两种固相组分之间的关系式：

积分后获得下式：

对于生物质原料颗粒m_c0＝0，令

m_c＝x_c(m_b0-m_b) (49)

结合颗粒相组分质量分数表达式进行推导：

Y_bm_b+Y_bm_c＝m_b (52)

Y_bm_b+Y_bx_cm_b0-Y_bx_cm_b＝m_b (53)

m_b0＝ρ_bV₀ (55)

式中V₀为生物质原料颗粒的体积；

对反应过程中进行颗粒尺度质量守恒分析获得下式：

则粒径变化数学模型为：

(2)数值解法

将方程44和45左右两边同时除以生物质原料质量m_b0，同时令：

则方程44和45转化为：

运用Matlab软件对方程60和61进行数值求解，并由下式计算颗粒总质量和生物质组分的质量分数：

具体数值解结果如表2所示：

表2数值解结果

注：上表中E表示科学记数法。

将上表中w_b+w_c和Y_b数据在Origin作图，参照图2：

对图2中的曲线进行数据拟合获得下式：

与密度变化数学模型结合进行颗粒尺度质量守恒分析：

则粒径变化数学模型为：

将数值解和解析解两种方法获得的粒径变化数学模型作图对比，如图3所示，可见两种方法所计算的粒径几乎一致，说明两种方法均可有效描述颗粒粒径的变化规律。

步骤三、建立颗粒类型分段曳力模型

对于A类颗粒，采用EMMS曳力模型：

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

H_D为非均匀系数，由下式获得：

H_D＝a(Re_pi+b)^c (69)

当0.4≤α_g≤0.46，

当0.46<α_g≤0.545，

当0.545<α_g≤0.99，

当0.99<α_g≤0.9997，

当0.9997≤α_g≤1，

对于B类颗粒，气相与颗粒相之间的曳力系数采用Gidaspow曳力模型：

当α_g>0.8，

当α_g≤0.8，

将方程67和方程75、76代入方程41中，即可获得颗粒类型分段曳力模型。

步骤四、流化床内流动反应特性的预测

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

模拟计算完成可获得流化床内的流动反应特性，其中生物质的体积分数、密度和粒径分布如图5所示。

由以上实施例可知，本发明采用计算流体力学方法对颗粒密度和粒径动态变化的流化床进行模拟研究，结合描述颗粒相密度和粒径变化规律的数学模型对颗粒的密度和粒径进行实时修正，通过颗粒类型分段曳力模型描述多种颗粒类型并存体系的气-颗粒相间作用力，从而预测颗粒密度和粒径动态变化流化床内的流动反应特性。该方法无需进行大量复杂昂贵的实验研究，可以节约大量的人力、物力及时间成本，为流化床的性能预测、优化控制及设计放大提供理论基础。

Claims (3)

1.一种颗粒密度和粒径动态变化流化床的模拟方法，其特征在于，包括以下步骤：

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

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

步骤三、建立颗粒类型分段曳力模型；

步骤四、流化床内流动反应特性的预测；

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

步骤二具体为：

当颗粒在流化床内发生非均相化学反应时，多组分多步平行-顺序反应，由以下反应式表示：

式中A、B、C、D、E、F、G、H、I、J、K、L为反应物和产物，下标n表示第n种组分，括号内s和g分别表示颗粒相组分和气相组分，k1_n、k2_n、k3_n、k4_n、k5_n为相应的速率常数；

对于密度和粒径动态变化的颗粒，其颗粒相密度变化数学模型可由下式表示：

式中Y_i和ρ_i分别为颗粒相中i组分的质量分数和密度；

而颗粒相粒径变化数学模型，则可根据颗粒尺度的质量守恒，由颗粒相的质量和密度变化数学模型获得，因此，需确定颗粒相的质量；

依据化学反应工程理论对反应式进行分析，对各颗粒相组分的生成\消耗速率建立一系列偏微分方程组，如下所示：

通过求解以上的偏微分方程组，即可获得颗粒相各组分的质量，从而可确定颗粒相的总质量和各组分的质量分数，最终根据综合密度变化数学模型方程(31)和颗粒尺度的质量守恒获得颗粒相粒径变化数学模型；

在求解偏微分方程时可采取两种方法，一种是求解析解，适用于相对较为简单的热解机理；另一种是求数值解，适用于极其复杂的热解机理；

(1)解析解法

该方法首先将各颗粒相组分的生成/消耗速率相比，即可消除反应时间项，获得各颗粒相组分之间的关系式；再根据颗粒相组分质量分数和密度变化数学模型，结合颗粒尺度的质量守恒分析获得粒径变化数学模型；

(2)数值解法

对于极其复杂的热解机理，由于组分数目多且反应步骤复杂，难以获得反应速率偏微分方程组的解析解；但是可运用数学软件对偏微分方程进行数值求解，获得不同反应时间下各颗粒相组分的质量，进而求得不同反应时间下颗粒质量及各颗粒相组分的质量分数，在数学软件中通过数据拟合的方法获得颗粒质量与各颗粒相组分质量分数的关联式，最后与密度变化数学模型结合进行颗粒尺度质量守恒分析获得粒径变化数学模型；

以上两种方法所建立粒径变化数学模型的通用形式为：

d_p＝f(Y_i,ρ_i)d_p0

式中d_p0为颗粒的初始粒径。

2.根据权利要求1所述的一种颗粒密度和粒径动态变化流化床的模拟方法，其特征在于，步骤一具体为：

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

气相连续性方程：

颗粒相连续性方程：

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

气相动量方程：

颗粒相动量方程：

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

气相能量方程：

颗粒相能量方程：

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

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

其中c_p为比热容，μ为粘度，d为颗粒直径；

气相组分输送方程：

其中Y_gi为气相中组分i的质量分数，D_i,m为气相中组分i的扩散系数，R_gi为气相中组分i的均相反应速率，S_rgi为气相中组分i的非均相反应速率；

颗粒相组分输送方程：

其中Y_pi,j为第i种颗粒相中组分j的质量分数，D_j,m为第i种颗粒相中组分j的扩散系数，R_pi,j为第i种颗粒相中组分j的均相反应速率，S_rpi,j为第i种颗粒相中组分j的非均相反应速率；

颗粒拟温度方程：

其中Θ为颗粒拟温度；

气相剪应力：

颗粒相剪应力：

其中I为单位张量，σ为颗粒体相黏度；

颗粒相压力：

其中e为碰撞恢复系数；

径向分布函数：

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

当

否则

其中

颗粒体相黏度：

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

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

颗粒拟温度的输送系数：

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

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

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

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

R_i＝kΠ(C_i)ⁿ (29)

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

3.根据权利要求1所述的一种颗粒密度和粒径动态变化流化床的模拟方法，其特征在于，步骤三具体为：

对于不同类型的颗粒，特别是A类和B类颗粒，所采用的气-颗粒相曳力模型差异较大；因此，如果反应器内同时存在A类和B类颗粒，则需建立颗粒类型分段曳力模型，对不同类型的颗粒分别采用相应的曳力模型，并且需在一定程度上避免数值的跳跃：

(1)计算无因次量

(2)计算A类和B类颗粒的分界值：

(3)建立颗粒类型分段曳力模型

对于A类和B类颗粒分别选用相应的曳力模型，其曳力系数以β_A和β_B表示，则颗粒类型分段曳力模型为：

其中：