阅读说明：本技术 基于单峰热解曲线的炭化可燃物热解动力学参数计算方法 (carbonized combustible pyrolysis kinetic parameter calculation method based on unimodal pyrolysis curve ) 是由 丁彦铭 张文龙 毛少华 黄必清 于 2019-07-17 设计创作，主要内容包括：本发明公开了一种基于单峰热解曲线的炭化可燃物热解动力学参数计算方法,包括：等转化率法FWO法求取可燃物单峰热解的活化能值E；假设多种反应机理,通过CR法求取活化能值E与FWO法得到的值对比,确定可燃物单峰热解的反应机理函数f(α),求取指前因子A；确定求解参数的取值范围；粒子群人工智能优化算法(PSO)对指前因子A、活化能值E以及炭生成率v进行寻优,并在最优解时输出计算结果。本发明降低了求可燃物单峰热解高阶偏微分方程的难度；计算精度好；优化效率高；使用范围广,可以同时应用到多组升温速率和多种可燃物的单峰热解求解中,有效预测实验结果。(the invention discloses a carbonized combustible pyrolysis kinetic parameter calculation method based on a unimodal pyrolysis curve, which comprises the following steps: obtaining the activation energy value E of unimodal pyrolysis of the combustible by an equal conversion method FWO method; assuming a plurality of reaction mechanisms, comparing the activation energy value E obtained by a CR method with a value obtained by a FWO method, determining a reaction mechanism function f (alpha) of unimodal pyrolysis of the combustible, and obtaining a pre-exponential factor A; determining the value range of the solving parameter; and optimizing the pre-pointing factor A, the activation energy value E and the carbon generation rate v by a particle swarm artificial intelligence optimization algorithm (PSO), and outputting a calculation result in the optimal solution. The method reduces the difficulty of solving the high-order partial differential equation of the unimodal pyrolysis of the combustible; the calculation precision is good; the optimization efficiency is high; the application range is wide, and the method can be simultaneously applied to solving unimodal pyrolysis of multiple groups of heating rates and multiple combustibles, so that the experimental result can be effectively predicted.)

基于单峰热解曲线的炭化可燃物热解动力学参数计算方法

技术领域

本发明涉及一种基于单峰热解曲线的炭化可燃物热解动力学参数计算方法，具体涉及炭化可燃物单峰热解动力学参数的求解和优化，通过编写粒子群人工智能优化算法(PSO)脚本主程序并设置算法参数，运行PSO脚本主程序，迭代计算，搜索各参数最优解和偏差度，属于炭化可燃物单峰热解动力学参数求解及优化领域。

背景技术

在世界能源危机的大背景下，不可再生能源的消耗亟需开发新能源技术，通过炭化可燃物热解可以得到炭化可燃物、液体和气体产物，实现能源的可再生利用。比如热解可以将木质纤维素类生物质生成燃料气、生物油和生物炭等。例如，DOI为10.1016/j.enconman.2016.05.007的论文“Thermal degradation of beech wood withthermogravimetry/Fourier transform infrared analysis”指出炭化可燃物热解是能源转化的重要一环，具有重要的研究意义。其次，我们的生活离不开合成材料，生物质，建筑物等炭化可燃物物质，这些炭化可燃物在给我们带来方便的同时，也带来一定的危险隐患，火灾是最为常见的事故。热解作为火灾发生的第一步，引起越来越多的学者对其进行研究。最后，随着社会的发展，大量的炭化废弃可燃物得不到有效、合理的处理，目前的废弃物的处理方式包括：填埋、焚烧、热解，考虑到土地和环境保护，热解被认为是最有前景的废弃物处理方法。目前，炭化可燃物热解动力学的求解多存在于理论的数学计算，但是，大部分的炭化可燃物材料的热解，需要复杂的动力学机理来模拟它的热解过程，构造模型、模拟和计算难度大，准确性较低。而且，这些机理需要多个参数，而这些参数并不能通过简单的传统热解方法计算得到。

全局优化算法是计算、优化热解参数的有效方法，PSO作为改进的模拟法，首先编写PSO脚本主程序并设置算法参数，设置模拟输出值。然后运行该算法，在设置的范围内可以自动搜索在实验数据基础上获得的指前因子A，活化能E，炭生成率v，搜索各参数最优解和偏差度。算法提供的偏差度表示预测结果和实验值的偏差程度，该偏差度越小越好，接近于0最佳。例如，DOI为10.1016/j.energy.2019.04.030的论文“The accuracy andefficiency of GA and PSO optimization schemes on estimating reaction kineticparameters of biomass pyrolysis”在实验值的基础上构造计算模型，对炭化可燃物热解动力学参数进行计算和优化，计算时间短，精度高，预测值与实验值高度一致。因此本发明提出一种用于炭化可燃物单峰热解的化学反应动力学参数求解及优化方法，有效降低了求解炭化可燃物单峰热解高阶偏微分方程的难度和复杂程度；计算精度好；优化效率高；使用范围广，可以同时应用到多组升温速率和多种炭化可燃物的单峰热解求解中，有效预测实验结果。

发明内容

本发明针对现有炭化可燃物单峰热解动力学参数求解存在的问题，提供一种基于单峰热解曲线的炭化可燃物热解动力学参数计算方法，该方法降低了求解炭化可燃物单峰热解高阶偏微分方程的难度和复杂程度；计算精度好；优化效率高；使用范围广，可以同时应用到多组升温速率和多种炭化可燃物的单峰热解求解中，有效计算出结果。

本发明解决其技术问题，所采用的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法，包含如下步骤：

S1、基于单峰热解曲线的炭化可燃物热解单步反应公式：固体→v炭+(1-v)挥发物，v为炭生成率，采用等转化率法FWO在不同的转化率α下求取炭化可燃物热解的活化能值E，最后求取不同转化率下活化能值E的平均值；

S2、在多种反应机理下，采用CR法求取在不同机理下炭化可燃物热解的活化能平均值E，与FWO法得到的值进行对比，找到FWO法计算的活化能值最接近CR法计算的活化能值时CR法对应的机理就是炭化可燃物单峰热解的反应机理函数f(α)，继而求得指前因子A；

S3、将指前因子A以及最接近的活化能值E分别上下浮动一定百分比，分别得到指前因子A和最接近的活化能值E的寻优范围，炭生成率v的寻优范围设置为(0,1)；

S4、在步骤S3确定的寻优范围的基础上，采用粒子群人工智能优化算法PSO程序对指前因子A、最接近的活化能值E以及炭生成率v进行寻优，并在最优解时输出转化率α、转化速率dα/dt、质量损失m/m₀以及质量损失速率d(m/m₀)/dt，从而完成炭化可燃物热解动力学参数的计算；寻优过程中，参数的迭代更新公式如下：

α_i＝α_i-1+(dα/dt)_i-1×(t_i-t_i-1)，

(m/m₀)_i＝1-α_i×(1-v)，

(d(m/m₀)/dt)_i＝((m/m₀)_i-(m/m₀)_i-1)/(t_i-t_i-1)，

dα/dt＝Af(α)exp(-E/RT)；

其中t表示时间，m、m₀表示炭化可燃物热解瞬时质量和初始质量，R为普适气体常数，T为绝对温度，i表示更新迭代次数。

进一步地，在本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法中，步骤S1具体为：炭化可燃物单峰热解单步反应公式：固体→v炭+(1-v)挥发物，v为炭生成率，在炭化可燃物单峰热解反应机理未知的情况下，通过等转化率法FWO法，在不同的升温速率下，测量与实验结果对应的固定转化率的温度，通过公式：lnβ＝ln(AE/Rg(α))-5.331-1.052(E/RT)，其中R表示普适气体常数，β表示升温速率，g(α)表示微分反应机理函数，T表示绝对温度，根据公式在多种不同的转化率下，得到lnβ与1/T的直线斜率-1.052(E/R)，得到炭化可燃物单峰热解的活化能值E，最后求取各种转化率下活化能E的平均值。

进一步地，在本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法中，所述步骤S2具体为：分别在炭化可燃物单峰热解的多种不同的反应机理下，采用CR法求取炭化可燃物单峰热解的活化能值E，也就是一个反应机理对应一个活化能值，通过公式：ln(g(α)/T²)＝ln(AR/βE)-(E/RT)，在各个热解反应机理下计算多个升温速率的ln(g(α)/T²)与1/T的斜率得到炭化可燃物热解单峰的活化能值E，并求取多组升温速率活化能的平均值，通过CR法计算的平均活化能值与FWO法计算的平均活化能值比较，找到FWO法计算的活化能值最接近CR法计算的活化能值时CR法对应的机理就是炭化可燃物单峰热解的反应机理，确定反应机理函数f(α)，在确定后机理下，通过公式的截距，ln(AR/βE)可以得到指前因子A。

进一步地，在本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法中，步骤S3中，上下浮动一定百分比是指上下各自浮动50％。

进一步地，在本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法中，步骤S4具体为：PSO寻优时，适应度函数值Φ为预测值和实验值的偏差程度，具体计算公式为：

其中，φ_m、φ_mlr、φ_α及φ_dα/dt分别表示质量损失、质量损失速率、转化率、转化速率的目标函数；N表示实验次数；n表示每个实验的实验数据点的个数；CML_mod、CML_exp表示累积质量损失的模拟值和实验值；MLR_mod、MLR_exp表示质量损失速率的模拟值和实验值；α_mod，α_exp表示转化率的模拟值和实验值；dα/dt_mod、dα/dt_exp表示转化速率的模拟值和实验值；w_CML、w_MLR、w_α、w_dα/dt表示质量损失、质量损失速率、转化率转化速率的权重系数。

进一步地，在本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法中，步骤S4中各个实验值是指步骤S1中进行FWO法所测量的实验值的计算公式如下：

CML_exp＝m/m₀

MLR_exp＝d(m/m₀)/dt＝

((m/m₀)_i-(m/m₀)_i-1)/(t_i-t_i-1)

dα/dt＝Af(α)exp(-E/RT)；

其中，m和m_t表示反应瞬时的质量，m₀表示样品初始质量，m_∞表示最终的质量。

相对于现有热解求解方法，本发明的有益效果如下：1)本发明提供的一种基于单峰热解曲线的炭化可燃物热解动力学参数求解方法，以炭化可燃物单峰热解实验数据为基础得到炭化可燃物单峰热解的反应机理，改进全局优化算法—粒子群优化算法(PSO)，降低了在理论分析的公式推导过程中，求解高阶偏微分方程的难度；该算法能够自动计算并绘制图表，在高次数的迭代计算中，计算速度快且误差较小；验证了实验结果的合理性；确定了理论构建模型的准确性并将之应用到更大尺寸的工况研究中；节省了人力物力财力；2)本发明提供的一种基于单峰热解曲线的炭化可燃物热解动力学参数求解方法，该算法得到模拟值的准确性较高；使用范围广，可以同时应用到多组升温速率和多种炭化可燃物的单峰热解求解中，有效预测实验结果；3)本发明提供的一种基于单峰热解曲线的炭化可燃物热解动力学参数求解方法，该算法在实验数据的基础上，可以模拟热解的复杂过程，在实验基础上贴合实际情况优化数据，有效提高实验结果的准确性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明所述的基于单峰热解曲线的炭化可燃物热解动力学参数求解方法的流程图；

图2为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的转化率α的模拟值和实验值对比图；

图3所示为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的转化速率dα/dt的模拟值和实验值对比图；

图4所示为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的质量损失m/m0的模拟值和实验值对比图；

图5所示为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的质量损失速率d(m/m0)/dt的模拟值和实验值对比图；

图6所示为本发明实施案例一中以外墙保温材料EPS单峰热解为例在氮气环境中热解得到的偏差度图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明提供一种基于单峰热解曲线的炭化可燃物热解动力学参数求解方法，其方法流程如图1所示，由图1可知，本发明的基于单峰热解曲线的炭化可燃物热解动力学参数计算方法包含如下步骤：

S1、基于单峰热解曲线的炭化可燃物热解单步反应公式：固体→v炭+(1-v)挥发物，v为炭生成率，采用等转化率法FWO(Flynn–Wall–Ozawa method)在不同的转化率α下求取炭化可燃物热解的活化能值E，最后求取不同转化率下活化能值E的平均值。

在本实施例中，步骤S1具体为：炭化可燃物单峰热解单步反应公式：固体→v炭+(1-v)挥发物，v为炭生成率，在炭化可燃物单峰热解反应机理未知的情况下，通过等转化率法FWO法，在不同的升温速率下，测量与实验结果对应的固定转化率的温度，通过公式：lnβ＝ln(AE/Rg(α))-5.331-1.052(E/RT)，其中R表示普适气体常数，β表示升温速率，g(α)表示微分反应机理函数，T表示绝对温度，根据公式在多种不同的转化率下，得到lnβ与1/T的直线斜率-1.052(E/R)，得到炭化可燃物单峰热解的活化能值E，最后求取各种转化率下活化能值E的平均值。

S2、在多种反应机理下，采用CR(Coats-Redfern method)法求取在不同机理下炭化可燃物热解的活化能平均值E，与FWO法得到的值进行对比，找到FWO法计算的活化能值最接近CR法计算的活化能值时CR法对应的反应机理函数f(α)，继而求得指前因子A。

在本实施例中，步骤S2具体为：分别在炭化可燃物单峰热解的多种不同的反应机理下，采用CR法求取炭化可燃物单峰热解的活化能值E，也就是一个反应机理对应一个活化能值通过公式：ln(g(α)/T²)＝ln(AR/βE)-(E/RT)，在假设的多个热解反应机理下计算多个升温速率的ln(g(α)/T²)与1/T的斜率得到炭化可燃物热解单峰的活化能值E，并求取多组升温速率活化能的平均值，通过CR法计算的平均活化能值与FWO法计算的平均活化能值比较，找到FWO法计算的活化能值最接近CR法计算的活化能值时CR法对应的机理就是炭化可燃物单峰热解的反应机理，确定反应机理函数f(α)，在确定机理后，通过公式的截距，ln(AR/βE)可以得到指前因子A。

S3、将指前因子A以及活化能值E(FWO求出的活化能值E)分别上下浮动一定百分比，分别得到指前因子A和活化能值E的寻优范围，根据炭化可燃物热解主要生成炭、液体和气体产物，炭化可燃物单峰热解单步反应公式：固体→v炭+(1-v)挥发物，所以，炭生成率v在0-1之间，即炭生成率v的寻优范围设置为(0,1)，例如上下浮动一定百分比是指上下各自浮动50％，以A为例：即下限50％*A，上限150％*A。

S4、在步骤S3确定的寻优范围的基础上，采用粒子群人工智能优化算法PSO程序对指前因子A、活化能值E以及炭生成率v进行寻优，并在最优解时输出转化率α、转化速率dα/dt、质量损失m/m₀以及质量损失速率d(m/m₀)/dt，从而完成炭化可燃物热解动力学参数的计算；寻优过程中，参数的迭代更新公式如下：

α_i＝α_i-1+(dα/dt)_i-1×(t_i-t_i-1)，

(m/m₀)_i＝1-α_i×(1-v)，

(d(m/m₀)/dt)_i＝((m/m₀)_i-(m/m₀)_i-1)/(t_i-t_i-1)，

dα/dt＝Af(α)exp(-E/RT)；

其中t表示时间，m、m₀表示炭化可燃物热解瞬时质量和初始质量，R为普适气体常数，T为绝对温度，i表示更新迭代次数。

PSO寻优时，适应度函数值Φ为预测值和实验值的偏差程度，具体计算公式为：

其中，适应度函数值Φ越小越好，φ_m、φ_mlr、φ_α及φ_dα/dt分别表示质量损失、质量损失速率、转化率、转化速率的目标函数；N表示实验次数；n表示每个实验的实验数据点的个数；CML_mod、CML_exp表示累积质量损失的模拟值和实验值；MLR_mod、MLR_exp表示质量损失速率的模拟值和实验值；α_mod，α_exp表示转化率的模拟值和实验值；dα/dt_mod、dα/dt_exp表示转化速率的模拟值和实验值；w_CML、w_MLR、w_α、w_dα/dt表示质量损失、质量损失速率、转化率转化速率的权重系数。这些实验值是指步骤S1中进行FWO法所计算的实验值，或者是步骤S2中进行CR法所计算得到的实验值，实验值的计算公式如下：

CML_exp＝m/m₀

MLR_exp＝d(m/m₀)/dt＝

((m/m₀)_i-(m/m₀)_i-1)/(t_i-t_i-1)

dα/dt＝Af(α)exp(-E/RT)；

其中，m和m_t表示反应瞬时的质量，m₀表示样品初始质量，m_∞表示最终的质量。

具体实施例：

本实施例提供基于单峰热解曲线的炭化可燃物热解动力学参数求解及优化方法的实施案例中以外墙保温材料EPS单峰为例在氮气环境中热解。

步骤S1中外墙保温材料EPS单峰在氮气环境中热解，该热解单步反应公式：

固体→v炭+(1-v)挥发物，v为炭生成率。在不掌握EPS单峰反应机理时，通过等转化率法FWO法，在不同的升温速率5k/min，20k/min下，测量在0.10，0.15，0.20，0.25，0.30，0.35，0.40，0.45，0.50，0.55，0.60，0.65，0.70，0.75，0.80共15个转化率值时的温度。根据公式：lnβ＝ln(AE/Rg(α))-5.331-1.052(E/RT)，其中R表示普适气体常数，β表示升温速率，g(α)表示微分反应机理，T表示绝对温度。得到lnβ与1/T的直线斜率-1.052(E/R)，计算固定每个转化率时的活化能，并求取所有转化率的平均活化能值E为163.58kJ/mol，该活化能平均值则是FWO法得到的外墙保温材料EPS单峰热解的活化能值E。

步骤S2中外墙保温材料EPS单峰在氮气环境中热解，通过CR法公式：

ln(g(α)/T²)＝ln(AR/βE)-(E/RT)，假设13个反应机理，分别为：g(α)＝α；g(α)＝-ln(1-α)；g(α)＝[1-(1-α)1-n]/(1-n)；g(α)＝1-(1-α)^1/2；g(α)＝1-(1-α)^1/3；g(α)＝α^1/3；g(α)＝α^1/4；g(α)＝1-(1-α)^1/4；g(α)＝1-(1-α)²；g(α)＝1-(1-α)³；g(α)＝3/2[(1-α)-^1/3-1]-¹；g(α)＝[-ln(1-α)]-¹；g(α)＝1/2α。通过ln(g(α)/T²)与1/T的斜率-(E/R)分别计算13个机理下的2组升温速率(5k/min，20k/min)的活化能，并求取2组升温速率的平均活化能值E，通过CR法计算的平均活化能为158.17kJ/mol与FWO法计算的平均活化能比较(FWO法可以根据实验数据直接计算活化能值E，不需要知道反应机理，但该法不能往下继续求解反应机理等其它参数。但为了找到反应机理，所以才和CR法计算的活化能值进行比较。因为CR法是假设一些反应机理，在假设的多个反应机理下进行求解活化能值，每个活化能对应一个反应机理，所以两种方法比较活化能，FWO法与CR法的活化能接近时，因为CR法已经假设了反应机理，所以就认为实验的反应机理是活化能最接近时对应的反应机理)，找到两种方法计算的活化能最接近时其对应的机理函数f(α)，外墙保温材料EPS在氮气环境中的反应机理函数微分方程g(α)＝1-(1-α)³。在该机理下，通过CR法公式的截距ln(AR/βE)计算2组升温速率的指前因子A，并求得平均值，由于A的值较大，实施例1求取指前因子平均值A的对数，lnA为23.77[ln(s-¹)]。然后继续进行步骤S3、S4，完成数据的计算，具体可参见图2-图6，其中图2为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的转化率α的模拟值和实验值对比图，图3为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的转化速率dα/dt的模拟值和实验值对比图，图4为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的质量损失m/m₀的模拟值和实验值对比图，图5为本发明实施案例一中算法得到的以外墙保温材料EPS单峰热解为例在氮气环境中热解的质量损失速率d(m/m₀)/dt的模拟值和实验值对比图，图6为本发明实施案例一中以外墙保温材料EPS单峰热解为例在氮气环境中热解得到的偏差度图。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。