CN107677761A

CN107677761A - 一种建立模拟卷烟燃烧的计算流体力学模型的方法和系统

Info

Publication number: CN107677761A
Application number: CN201710891565.7A
Authority: CN
Inventors: 李巧灵; 李跃锋; 刘泽春; 邓其馨; 黄华发; 黄朝章; 张建平; 黄惠贞; 林艳; 谢卫
Original assignee: China Tobacco Fujian Industrial Co Ltd
Current assignee: China Tobacco Fujian Industrial Co Ltd
Priority date: 2017-09-27
Filing date: 2017-09-27
Publication date: 2018-02-09
Anticipated expiration: 2037-09-27
Also published as: CN107677761B

Abstract

本发明涉及一种建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的方法，包括：a)建立模拟卷烟燃烧场景的几何模型；b)建立模拟卷烟燃烧场景内物理和/或化学反应的方程，包括：建立烟丝热解反应动力学方程；建立烟丝热解产物燃烧反应动力学方程；建立卷烟纸燃烧反应动力学方程；建立卷烟纸渗透率随温度变化的方程；c)采用计算流体力学软件的求解软件(例如Fluent、CFX、Phoenics、Flow3d、Flowmaster)加载步骤a)的几何模型和步骤b)的方程，建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型。该方法建立的模型能够准确模拟卷烟燃烧时温度场和/或物质浓度场。

Description

一种建立模拟卷烟燃烧的计算流体力学模型的方法和系统

技术领域

本发明属于卷烟领域，具体涉及一种建立模拟卷烟燃烧的计算流体力学模型的方法和系统。

背景技术

在卷烟燃烧模型研究进展一文(中国烟草学报，2013(2):115-122.)中，李巧灵等对卷烟燃烧模型的研究进展进行了综述。该文介绍了目前文献中所报道的阴然和吸燃模型。对于阴燃过程，气体的流速是确定的，所以此过程发生的反应在化学计量学和动力学的控制范围，因此卷烟的阴燃模型研究得相对较多。

在卷烟阴燃过程的数值模拟一文(烟草科技，2014(6))中，李巧灵等对文献中已经报道的卷烟阴燃模型进行整合，并利用Fluent软件建立了一个相对完善的卷烟阴燃模型，该模型可模拟卷烟燃烧过程中不同时刻的温度场分布，卷烟烟气中氧气、一氧化碳、二氧化碳和水蒸气的浓度场分布，并将阴燃线性燃烧速度和卷烟内部最高温度模拟值与实验值进行对比。

发明内容

发明人发现，在卷烟燃烧过程中，气体区域的空气会通过卷烟纸进入烟丝区域，卷烟纸渗透率高，则气体区域的气体进入烟丝区域的阻力小，卷烟纸渗透率低则气体区域的空气进入烟丝区域的阻力大。

发明人进一步发现，卷烟纸的渗透率随卷烟纸的燃烧状态发生改变，已燃烧的卷烟纸和未燃烧的卷烟纸的渗透率是不同的，因此，在卷烟燃烧模型中需要模拟卷烟纸的燃烧过程，进而确定哪一部分卷烟纸为已燃烧卷烟纸，哪一部分卷烟纸为未燃烧卷烟纸，进而准确地设定不同区域的卷烟纸的渗透率。

发明人进一步发现，卷烟纸除了已燃烧状态和未燃烧状态，还有一部分靠近燃烧端的卷烟纸处于受热变(焦)黄的状态，发明人创造性地发现，对于受热变(焦)黄的卷烟纸，其渗透率既不同于已燃烧卷烟纸，也不同于未燃烧卷烟纸，准确设定受热变(焦)黄的卷烟纸的渗透率，能够进一步提高模拟结果的准确性。

本发明第一方面提供一种建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的方法：

a)建立模拟卷烟燃烧场景的几何模型；

b)建立模拟卷烟燃烧场景内物理和/或化学反应的方程，包括：

建立烟丝热解反应动力学方程；

建立烟丝热解产物燃烧反应动力学方程；

建立卷烟纸燃烧反应动力学方程；

建立卷烟纸渗透率随温度变化的方程；

c)采用计算流体力学软件的求解软件(例如Fluent、CFX、Phoenics、Flow3d、Flowmaster)加载步骤a)的几何模型和步骤b)的方程，建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型。

本发明又一方面提供一种用于建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的系统，包括：

几何模型建立模块(例如Gambit软件)，其用于建立模拟烟草燃烧场景的几何模型；

物理和/或化学反应方程建立模块，该模块包括

烟丝热解反应动力学方程建立模块；

烟丝热解产物燃烧反应动力学方程建立模块；

卷烟纸燃烧反应动力学方程建立模块；

卷烟纸渗透率随温度变化的方程建立模块；

计算流体力学模型建立模块，其采用计算流体力学软件的求解软件加载几何模型建立模块建立的几何模型和物理和/或化学反应方程建立模块建立的方程，建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型。

在一个实施方案中，烟丝热解反应动力学方程建立模块用于建立烟丝热解反应动力学方程。

在一个实施方案中，烟丝热解产物燃烧反应动力学方程建立模块用于建立烟丝热解产物燃烧反应动力学方程。

在一个实施方案中，卷烟纸燃烧反应动力学方程建立模块用于建立卷烟纸燃烧反应动力学方程。

在一个实施方案中，卷烟纸渗透率随温度变化的方程建立模块用于建立卷烟纸渗透率随温度变化的方程。

在一个实施方案中，建立卷烟纸渗透率随温度变化的方程的方法包括：

设定第一温度的卷烟纸具有第一渗透率；

设定第二温度的卷烟纸具有第二渗透率；

设定第三温度的卷烟纸具有第三渗透率；

优选地，

第一温度为a K以下，第一渗透率为0.5×10^-15～5×10^-15m²；

第二温度为a～b K，第二渗透率为1.5×10^-9～6×10^-9m²；

第三温度为b K以上，第三渗透率为0.5×10⁵～5×10⁵m²；

a＝450～500K(例如473K)，b＝600～650K(例如623K)；

优选地，设定卷烟纸渗透率为单向变化，只增大不减小。

在一个实施方案中，第一温度的卷烟纸为未燃烧卷烟纸。

在一个实施方案中，第二温度的卷烟纸为受热变黄的卷烟纸。

在一个实施方案中，第三温度的卷烟纸为已燃烧卷烟纸。

在一个实施方案中，建立卷烟纸燃烧反应动力学方程的方法包括：

-在含氧气氛下加热卷烟纸样品，采用热重分析法检测卷烟纸样品的质量变化，测得卷烟纸受热时的微分热重曲线；

-对所述微分热重曲线进行分峰拟合，将曲线分为m个单峰曲线；

-将m个(例如2、3、4或5个)单峰曲线拟合为m个反应j的动力学

方程：

各反应j的动力学方程如下：

j＝1～m，j和m为正整数；

a_z，j为反应j的转化率、为反应j的转化率对时间的导数，T_z为加热卷烟纸的温度，R为理想气体常数，A_z，j为反应j的指前因子，E_z，j为反应j的活化能，n_z，j为反应j的反应级数；

a_z，j、T_z由卷烟纸受热时的微分热重曲线获得，A_z，j、E_z，j和n_z，j通过拟合获得。

在一个实施方案中，建立烟丝热解产物燃烧反应动力学方程的方法包括：

建立对应第一氧气浓度的第一烟丝热解产物燃烧反应动力学方程；

建立对应第二氧气浓度的第二烟丝热解产物燃烧反应动力学方程；以及

可选地，建立对应第二氧气浓度的第二烟丝热解产物燃烧反应动力学方程。

-在s种(例如2、3、4、5、6或7种)氧气浓度的含氧气氛下加热烟丝热解产物，采用热重分析法检测烟丝热解产物的质量变化，分别测得s个烟丝热解产物受热时的微分热重曲线，s为正整数(优选k≥2)；

-将s个微分热重曲线分别拟合为以下反应k的动力学方程：

k为正整数；

α_c，k为反应k的转化率，为反应k的转化率α_c，k对时间t的导数，T_c为加热温度，R为理想气体常数，E_c，k为反应k的活化能，n_c，k为反应k的反应级数，ρo₂为氧气密度；

a_c，k和T根据微分热重曲线获得，A_c、E_c，k和n_c，k通过拟合获得。

在一个实施方案中，s种氧气浓度在1～25vol％范围内；

在一个实施方案中，相邻两种氧气浓度相差至少2vol％；

在一个实施方案中，s≥k≥2，s和k为正整数；

在一个实施方案中，k＝2、3、4、5、6、7或8；

在一个实施方案中，s＝2、3、4、5、6、7或8；

在一个实施方案中，s种氧气浓度包括以下7种氧气浓度：0.8～1.2vol％、1.5～2.5vol％、2.8～3.5vol％、4～6vol％、8～12vol％、16～18vol％、19～22vol％；

在一个实施方案中，建立烟丝热解反应动力学方程的方法包括：

-在非氧化气体气氛(例如氮气氛)下加热烟丝样品，采用热重分析法检测烟丝样品的质量变化，绘制烟丝样品的微分热重曲线；

-对所述微分热重曲线进行分峰拟合，将曲线分为n个单峰曲线；

-将n个(例如2、3、4、或5个)单峰曲线拟合为n个反应i的动力学

方程：

i＝1、2…n，n和i均为整数；

a_v，i为反应i的转化率、为反应i的转化率对时间的导数，T_v为加热烟丝样品的温度，β为加热烟丝样品的升温速率，R为理想气体常数，A_v，i为反应i的指前因子，m_v，i为反应i的升温速率修正参数，E_v，i为反应i的活化能，n_v，i为反应i的反应级数；

T_v、β由烟丝样品的微分热重曲线获得，A_v，i、m_v，i、E_v，i和n_v，i通过拟合获得。

在一个实施方案中，步骤b)还包括，建立对应不同氧气浓度的有害成分生成量与温度的关系方程。

在一个实施方案中，建立对应不同氧气浓度的有害成分释放量与温度的关系方程的方法包括：

在不同氧气浓度的气氛下加热烟丝样品；

检测烟丝样品受热时有害成分的释放量与加热温度的关系；

根据上一步的检测结果拟合对应不同氧气浓度的有害成分释放量与温度的关系方程；

优选地，有害成分为CO或焦油。

在一个实施方案中，步骤b)还包括建立滤棒对焦油的截留率方程的步骤：

在一个实施方案中，，滤棒对焦油的截留率方程包括：

e_s＝E_IN+E_IM+E_D+E_ID；

其中，E_IN、E_IM、E_D、E_ID是中间参数；

J＝(29.6-28ω^0.62)R_L ^2.8；

e为滤棒对焦油的截留率，e_s为滤棒的单纤维丝束对焦油的截留率，L_f是滤棒的长度，d_f是滤棒纤维丝束的直径，ω是滤棒内纤维丝束的体积分数，D_t是滤棒的总旦数，D_s是滤棒的单旦数，C_fiber是滤棒中卷曲纤维的比率，S_filter是滤棒的横截面积，ρ_g为气体密度，μ_g气体粘度，v_i是i方向上的气体速度，d_c为卷烟烟气气溶胶颗粒的直径，d_f为滤棒纤维丝束的直径，D_k为颗粒的扩散系数，K_B是波兹曼常数，T_filter是滤棒的温度。

在一个实施方案中，步骤d)还包括，采用计算流体力学软件的求解软件加载卷烟抽吸程序的步骤；

在一个实施方案中，卷烟抽吸程序设定如下：设定卷烟区域被抽吸的位置设置为速度进口，其它位置设置为压力出口，每60s进行一次抽吸周期，卷烟抽吸气体流速抽吸周期包括以下a)和b)：

a)阴燃58s，v＝0；

b)抽吸2s，v的单位为m/s。

在一个实施方案中，建立卷烟燃烧几何模型的的方法包括：

采用计算流体力学软件的前处理软件(例如Gambit、ICEM CFD)建立卷烟燃烧场景几何模型，所述几何模型包括气体区域和卷烟区域，卷烟区域位于气体区域内，卷烟区域包括烟丝区域和卷烟纸区域；

优选地，卷烟区域还包括滤棒区域；

优选地，卷烟燃烧场景几何模型是二维模型。

在一个实施方案中，本发明方法还包括，设定并用计算流体力学软件的求解软件加载以下参数和方程的步骤：

-烟丝的物理或化学参数；

-卷烟纸的物理或化学参数；

-气体的物理或化学参数；

-卷烟区域和气体区域的质量传递方程；

-卷烟区域和气体区域的动量传递方程；

-卷烟区域和气体区域的能量传递方程。

在一个实施方案中，所述物质浓度场包括以下一项或多项：

氧气浓度场；

烟丝热解产物浓度场；

卷烟纸燃烧产物浓度场；

有害成分浓度场。

本发明又一方面提供一种预测卷烟燃烧时温度和/或物质浓度的方法，根据权利要求1～10任一项所述的方法建立所述计算流体力学模型，运行所述模型，根据模型模拟结果获得卷烟燃烧时的温度和/或物质浓度。

几何模型建立模块，其用于建立模拟烟草燃烧场景的几何模型；

物理和/或化学反应方程建立模块，该模块包括

烟丝热解反应动力学方程建立模块；

烟丝热解产物燃烧反应动力学方程建立模块；

卷烟纸燃烧反应动力学方程建立模块；

卷烟纸渗透率随温度变化的方程建立模块；

在一个实施方案中，所述系统还包括烟丝热解反应热重分析模块，其用于在不含氧气氛中对烟丝进行热重分析，并将数据输出到烟丝热解反应动力学方程建立模块。

在一个实施方案中，所述系统还包括烟丝热解产物燃烧反应热重分析模块，其用于在含氧气氛中对烟丝热解产物进行热重分析，并加数据输出到烟丝热解产物燃烧反应动力学方程建立模块。

在一个实施方案中，所述系统还包括卷烟纸燃烧反应热重分析模块，器用于在含氧气氛中对卷烟纸燃烧反应进行热重分析，并将数据输出到卷烟纸燃烧反应动力学方程建立模块。

在一个实施方案中，建立卷烟纸燃烧反应热重分析模块通过在含氧气氛下加热卷烟纸样品，采用热重分析法检测卷烟纸样品的质量变化，测得卷烟纸受热时的微分热重曲线；

卷烟纸燃烧反应动力学方程建立模块通过对所述微分热重曲线进行分峰拟合，将曲线分为m个单峰曲线；

-将m个(例如2、3、4或5个)单峰曲线拟合为m个反应j的动力学

方程：

各反应j的动力学方程如下：

j＝1～m，j和m为正整数；

a_z，为反应j的转化率、为反应j的转化率对时间的导数，T_z为加热卷烟纸的温度，R为理想气体常数，A_z，j为反应j的指前因子，E_z，j为反应j的活化能，n_z，j为反应j的反应级数；

一个用于建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的系统，物理和/或化学反应方程建立模块还包括烟丝有害成分释放量方程建立模块，其用于建立对应不同氧气浓度的烟丝有害成分释放量与温度的关系方程。

在一个实施方案中，所述系统还包括烟丝有害成分释放量检测模块，其用于：

在不同氧气浓度的气氛下加热烟丝样品；

检测烟丝样品受热时有害成分的释放量与加热温度的关系；

并将检测结果输出到烟丝有害成分释放量方程建立模块。

一个用于建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的系统，，物理和/或化学反应方程建立模块还包括滤棒对焦油截留率方程建立模块，其被配制为按以下参数和方程建立滤棒对焦油截留率方程：

e_s＝E_IN+E_IM+E_D+E_ID；

其中，E_IN、E_IM、E_D、E_ID是中间参数；

J＝(29.6-28ω^0.62)R_L ^2.8；

在一个实施方案中，所述方法还包括设定以下参数的步骤并用计算流体力学软件加载的步骤：

-设定烟丝区域的烟丝的物理化学参数，例如烟丝为多孔介质；

-设定卷烟纸区域的卷烟纸的物理化学参数，例如设定卷烟纸为多孔介质；

-设定气体区域的气体的物理化学参数，例如设定气体为气体流场。

术语解释

Fluent软件例如是ANSYS公司的计算流体力学(CFD)软件。

除非特别说明，温度单位为K。

有益技术效果：

1.采用300～800K/min的升温速率进行烟丝热解的热重分析实验；

2.在烟丝热解反应动力学方程中引入了参数：其中β为升温速率，m_v，i为升温速率修正参数；

3.在不同氧气浓度下，进行了烟丝热解产物燃烧的热重分析实验，根据所得实验数据拟合获得烟丝热解产物燃烧反应动力学方程；

4.建立卷烟纸热解燃烧动力学，引入了卷烟纸在不同温度下的渗透率参数；

5.根据烟丝受热时CO和焦油的释放量仅与温度和氧含量相关，建立了CO和焦油的释放量；

6.建立了滤棒对焦油的截留模型，预测了滤棒对焦油的截留率。

7.模型对温度的模拟较准确；

8.模型对物质浓度的模拟较准确。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的限定。在附图中：

图1为卷烟燃烧模型中四个区域的划分；

图2为牌号1烟丝粉末样品的DTG曲线(升温速率300K·min^-1)；

图3为牌号1烟丝在300Kmin^-1～800K·min^-1升温速率下的热解DTG曲线实验值和拟合值对比图；

图4为不同氧浓度下烟丝热解产物燃烧的DTG曲线实验值与计算值的对比图；

图5为卷烟纸热解燃烧DTG曲线实验值与计算值的对比图；

图6烟丝热解燃烧平台示意图；

图7为不同氧浓度下烟丝的CO释放量随温度变化曲线；

图8为不同氧浓度下烟丝焦油释放量随温度变化曲线

图9为不同氧浓度下烟丝的焦油释放实验的滤片照片；

图10为牌号1卷烟在60s和120s时，卷烟燃烧模型拟合的烟草热解产物浓度图(上)和卷烟纸燃烧线图(下)；

图11为牌号1卷烟在180s和240s时，卷烟燃烧模型拟合的烟草热解产物浓度图(上)和卷烟纸燃烧线图(下)；

图12为牌号1卷烟在300s和360s时，卷烟燃烧模型拟合的烟草热解产物浓度图(上)和卷烟纸燃烧线图(下)；

图13为牌号1卷烟在一次抽吸时，烟草热解产物浓度变化图和卷烟纸燃烧线变化图；

图14为牌号1卷烟在22mm、24mm、26mm、28mm、30mm、32mm、34mm、36mm处的温度曲线实验值和预测值

图15为牌号1卷烟滤棒在一次抽吸时(180s～182s)进出口的焦油释放量以及焦油截留率

图16为牌号1卷烟在180s时卷烟的气流速度场图；

图17为牌号1卷烟在181s时卷烟的气流速度场图；

图18为牌号1卷烟在182s时卷烟的气流速度场图。

具体实施方式

材料与仪器

卷烟烟丝(福建中烟工业有限责任公司，牌号1、牌号2、牌号3)。取适量卷烟烟丝，置于温度(22±1)℃、相对湿度(60±2)％的恒温恒湿箱中48h待用。

STA 449F3TG-DTA/DSC同步热分析仪(德国Netzsch公司)；TF-M100烟支燃吸气相温度场分析仪(北京紫东科技有限公司)；三孔道吸烟机；XP205电子天平(感量0.00001g，瑞士Mettler-Toledo公司)；GM200研磨仪(德国Retsch公司)；KBF720恒温恒湿箱(德国Binder公司)。J2KN多功能烟气分析仪(德国rbr公司)

一、基于Fluent建立卷烟燃烧(吸燃)模型

1.在Fluent软件中建立二维几何体模型

1.1划分网格区域

在二维模型中划分四个计算域，如图1所示：烟丝区域101；卷烟纸区域102；滤棒区域104；气体区域104。卷烟区域包括烟丝区域、卷烟纸区域和滤棒区域。

烟丝区域101为长方形。卷烟纸区域102包括两个长方形，卷烟纸区域102与烟丝区域的长度相同。卷烟纸区域102位于烟丝区域101的外侧，且卷烟纸区域102的一条长边与烟丝区域102的长边重合。滤棒区域104为长方体，滤棒区域104的短边与烟丝区域101的短边长度相同。滤棒区域104位于烟丝区域101的外侧，且滤棒区域104的一条短边与烟丝区域101的短边重合。

其中区域101、102、103为多孔介质，区域104为气体流场。气体区域的初始成分为23vol％氧气，余量为氮气。在卷烟烟丝区域进行烟丝的热解、燃烧过程的模拟运算；在卷烟纸区域进行卷烟纸燃烧过程及卷烟纸渗透率变化的模拟运算；在滤嘴区域进行滤棒对焦油的截留率的模拟运算。

如图1所示，卷烟燃烧模型中烟丝区域长0.058m(x轴)，宽0.008m(y轴)。滤棒区域长0.028m(x轴)，宽0.008m(y轴)。烟纸区域长0.058m，宽0.05mm，空气区域长0.115m，宽0.04m。

边界条件：抽吸位置为速度入口，速度值由抽吸曲线方程控制；其余边界均为压力入口，表压为0。

1.2各区域物理化学参数设置

各区域的材料按照卷烟各项的物理化学性质进行设置，具体见表1。

表1

气体区域的气体成分设置为23％氧气，余量为氮气。假设气体的密度、比热容、导热系数在卷烟吸燃过程中保持不变。

起始温度设为300K，烟气中一氧化碳、二氧化碳、水蒸气起始质量分数均设为0。

下文的步骤2～9的方程均运用UDF自定义编程(UDF即用户自定义函数，是用户自编的程序，可以对Fluent进行个性化设置。Fluent求解器可以动态地加载UDF。用Fluent6.3.26求解以上方程，采用分离求解器。压强-速度耦合算法采用了SIMPLE格式。用隐式格式对方程进行线化和求解。离散化采用一阶迎风格式计算通量。

2.建立烟丝热解、燃烧方程

假设烟丝区域烟丝的燃烧过程分为两个阶段，即烟丝热解阶段和烟丝热解产物燃烧阶段。下面通过对烟丝(福建中烟，牌号1)进行热重分析实验，根据实验数据，建立烟丝热解方程和烟丝热解产物燃烧方程。

2.1建立烟丝热解方程

2.1.1烟丝热解热重分析

取适量烟丝，置于温度(22±1)℃、相对湿度(60±2)％的恒温恒湿箱中48h。将烟丝研磨成粉，过筛选取粒径为198～165μm(80～100目)的烟粉作为烟丝样品。

称取9.50mg所述烟丝样品，平铺于同步热分析仪的托盘中，将热分析仪从室温升至873K，之后冷却至室温，得到烟丝热解产物。实验载气为氮气(N₂)，载气流速为50mL·min^-1。分别使用不同的升温速率(300K·min^-1、400K·min^-1、500K·min^-1、600K·min^-1、700K·min^-1和800K·min^-1)进行多次热分析实验。

2.1.2烟丝热解方程拟合

通过热重实验记录时间(min)、温度(K)、质量(mg)等物理量，得到样品质量损失率m(即损失质量占起始质量的质量百分数)对温度T的导数(wt％·K^-1)随温度T(K)的变化曲线，即微分失重曲线(微分热重曲线)。

图2为烟丝粉末样品受热时的微分热重曲线(升温速率300K·min^-1)。用Origin软件将微分热重曲线分成五个谷(R1，R2，R3，R4和R5)。假设五个谷分别对应五类主要成分i以及它们各自的独立平行的热解反应i(i＝1～5)。第一个谷(R1)例如代表水分的蒸发；第二到第五个谷例如分别代表不同热解前驱体的热解。第二个谷(R2)例如代表糖类，尼古丁，果胶和一些其他挥发物种的热解。第三个谷(R3)和第四个谷(R4)例如代表半纤维素和纤维素的热解。第五个谷(R5)例如代表以木质素作为主要成分的热解。

测算图2的微分热重曲线的五个谷R1、R2、R3、R4和R5的峰面积所占百分比，分别为9.52％、17.71％、18.04％、13.58％和41.16％。以该比值为依据对其他升温速率β(400K·min^-1、500K·min^-1、600K·min^-1、700K·min^-1和800K·min^-1)获得的热解微分热重曲线进行分峰。

假设五个谷分别对应五类组分i(i＝1～5)的五个独立平行反应i(i-1～5)，烟丝样品的可由下式表示：

方程(2.1)中，(i＝1～5)表示每类组分i的质量损失百分数对烟丝温度T的导数。

每类组分i的热解反应的转化率α_i(％)由下式(2.2)表示：

方程(2.2)中，m_i0和m_i∞是指每类组分i在热解反应起始和终止时的质量。T₀和T_i∞指在每类组分i在热解反应起始和终止时烟丝的温度。

对方程(2.2)进行整理可得：

对于反应i，其动力学方程f(a_v，i)可定义为：

式(2.4)中，α_v，i是反应i的转化率，A_v，i(min^-1)是反应i的指前因子，β(K/min)是加热烟丝的升温速率，m_v，i是反应i的升温速率修正参数，E_v，i(J/mol)是组分i的活化能，R(J/mol/K)是理想气体常数，T_v(K)是加热烟丝的温度。

假定f(α_v，i)＝(1-α_v，)^n，i，5个烟丝热解反应的拟合方程分别如式(2.5)所示：

对式(2.5)，a_v，i、T_v可从微分热重曲线中获得，用Origin软件非线性拟合功能中的整体拟合方法，将5个谷的数据同时进行拟合以寻求参数的最优解。该拟合过程中，Origin采用chi-square minimization(卡方最小化)的方法来确定最优动力学参数：热解反应活化能E_v，i、烟丝热解反应指前因子A_v，i、反应级数n_v，i和升温速率修正参数m_v，i，即最小化理论曲线与实验曲线的偏差。

所得的热解反应动力学参数如表2所示。表中f_v，i代表各组分的含量，由高斯分峰的峰面积比得到。

表2

图3示出红狼烟丝在300Kmin^-1～800K·min^-1升温速率下的热解DTG的实验曲线和拟合曲线对比图。如图3所示，在300Kmin^-1、400Kmin^-1、500Kmin^-1、600Kmin^-1、700Kmin^-1、800Kmin^-1，实验曲线与拟合曲线的吻合度都较高。

2.2建立烟丝热解产物的燃烧方程

2.2.1烟丝热解产物燃烧的热重实验

烟丝热解后，烟丝热解产物在氧气中发生氧化反应，即烟丝热解产物的氧化/燃烧。

为了建立烟丝热解产物的燃烧方程，进行如下实验。取适量红狼烟丝，置于温度(22±1)℃、相对湿度(60±2)％的恒温恒湿箱中48h。将烟丝研磨成粉，过筛选取粒径为198～165μm(80～100目)的烟粉作为实验样品。准确称取9.50mg样品，平铺于同步热分析仪的托盘中，热分析仪从室温升至873K，升温速率为10K·min^-1，实验载气为氮气，流速50mL·min^-1，此时得到烟草热解产物，待温度降到常温，然后根据需求改变载气中氧气的浓度(1％、2％、3％、5％、10％、15％、20％O₂，余量为N₂)，流速50mL·min^-1，热分析仪从室温升至873K，升温速率为10K·min^-1，得到烟丝热解产物在不同氧浓度下的质量损失。

图4示出了牌号1卷烟烟丝的热解产物在不同氧浓度下的燃烧DTG的实验值曲线和拟合值曲线。1％、2％、3％、5％、10％、15％、20％O₂氧气浓度下，实验值和拟合值曲线十分接近。

2.2.2烟丝热解产物燃烧方程拟合

参考上述的烟丝热解方程拟合过程，将热解产物的燃烧微分热重曲线拟合为方程。

假定烟丝热解产物氧化反应为1级反应。用独立平行反应模型描述烟丝热解产物燃烧的行为，烟丝热解产物燃烧方程如式(2.6)所示：

式(2.6)中，变量(kg·m^-3)代表O₂密度(0.28kg·m^-3)，T为加热温度，a_c为燃烧反应转化率，为燃烧反应转化率对时间t的导数，变量T、a_c、可由DTG实验得到，E_c为氧化活化能，A_c为氧化反应指前因子，n_c为对应不同氧气浓度的相关系数。

对于式2.6，用Origin软件非线性拟合功能中的整体拟合方法，该拟合过程中，Origin采用chi-square minimization(卡方最小化)的方法来确定最优动力学参数：氧化活化能E_c，氧化反应指前因子A_c，对应不同氧气浓度的相关系数n_c，i，即最小化理论曲线与实验点的偏差，结果见表3。

表3

图4示出了牌号1卷烟烟丝的热解产物在氧浓度1％～20％下的燃烧DTG的实验值曲线。

在编写UDF参数时，A_c、E_c和n_c的取值随氧气浓度变化而变化，当氧气浓度为0～2.5％时，A_c＝1.48E+07min^-1，E_c＝91.04KJ/mol，n_c＝1.09；当氧气浓度为2.5～7.5％时，A_c＝4.26E+07min^-1，E_c＝111.20KJ/mol，n_c＝0.957；当氧气浓度大于7.5％时，A_c＝8.30E+07min^-1，E_c＝116.31KJ/mol，n_c＝0.36。

建立卷烟纸热解燃烧方程：

将卷烟纸从牌号1卷烟中剥离，准确称取20.0mg样品，置于于同步热分析仪的托盘中，从室温升至773K，升温速率为10K·min-1。实验载气为空气(23v％氧气，余量为氮气)，流速50mL·min^-1。

通过利用同步热分析仪记录时间(min)、温度(K)、质量(mg)等物理量，得到样品质量损失率(wt％)随温度(K)的变化曲线，即失重曲线TG。对TG曲线求导，得到失重速率(％·min^-1)及微分失重曲线DTG。图5为卷烟纸样品在空气气氛下热解的TG和微分热重曲线。用Origin软件将微分热重曲线分成3个谷(R1、R2、R3)，即代表3个独立平行的热解反应。

2.3卷烟纸热解燃烧方程拟合

参考烟丝热解方程拟合过程，利用3个独立平行反应模型描述卷烟纸各组分的热解燃烧行为，卷烟纸各组分的热解燃烧的动力学方程可定义为：

方程(3.1)变量a_z，j、T_z、可由卷烟纸热解燃烧微分热重曲线获得，A_z，j(min^-1)是卷烟纸热解燃烧反应指前因子，E_z，j(kJ/mol)是卷烟纸热解燃烧反应活化能，R(J/mol/K)是理想气体常数，n_z，j是卷烟纸热解燃烧反应的反应级数。

用Origin软件非线性拟合功能中的整体拟合方法，将3个谷的数据同时进行拟合以寻求参数的最优解。该拟合过程中，Origin采用chi-square minimization(卡方最小化)的方法来确定最优动力学参数：卷烟纸热解燃烧反应活化能E_z，i，卷烟纸热解燃烧反应指前因子A_z，i和反应级数n_z，i，即最小化理论曲线与实验点的偏差，通过拟合确定上述3个反应各自的动力学参数，结果见表4。

表4

卷烟纸热解燃烧	R1	R2	R3
				f_z，j(％)	79.8	11.5	8.7
A_z，j(min^-1)	2.01E+08	1.031E+26	1.35E+45
				E_z，j(kJ/mol)	98.19	347.67	614.06
n_z，j	1.42	1.52	1.14
				R²	0.997	0.999	0.990

3.质量传递方程

3.1氧气在气相中的扩散系数

氧气在气相中的扩散系数D与气相温度有关：

D＝D₀(T_g/273)^1.75 (4.1)

D₀是O₂在273K、1个大气压下多孔介质中扩散系数的参考值，T_g是气相温度(由能量方程计算获得)。

D₀的值和烟丝孔隙率有关，烟丝孔隙率

其中，D_g是氧气的无限制扩散系数。在氮气气氛下，氧气的扩散系数

D_g＝2×10^-5m²s^-1。

3.2氧气质量传递方程源项

假定氧化反应发生在固相表面，氧气和气体产物进入气相。忽略气体从固相进入气相的边界层阻力。氧气的源项(kg·m^‐3·s^‐1)为：

式(4.3)中，是氧化反应中氧气的化学计量系数1.65，其它参数含义和取值参见表3。

3.3固相质量传递方程源项

3.3.1烟丝热解模型

将烟丝假设为五种热解前驱体(2.1.2中的组分i＝1～5)分别发生热解反应，热解反应方程可表示为：

其中，ρ_v，i是热解前驱体组分i的密度(kg/m³)，ρ_v是五种热解前驱体组分总密度，ρ_v，o五种热解前驱体组分起始总密度，即烟丝的起始密度740kg·m^-3；T_s(K)是固相温度(根据下文的能量方程计算获得)，A_v，i(min^-1)是热解前驱体i热解反应的指前因子，E_v，i(J/mol)是热解前驱体i热解反应的活化能，R(J/mol/K)是理想气体常数，f_v，i(％)是热解组分i的质量分数，n_vi是反应级数，取值参见表2。

3.3.2热解产物燃烧模型：

根据热解产物是热解前驱体组分i(i＝3、4、5)的热解产物，热解产物的燃烧反应方程如下：

此方程中右侧第一部分代表热解产生热解产物，第二部分代表热解产物氧化燃烧；其中f_c为热解前驱体转化为热解产物的质量转化比率系数，f_c＝0.3083；ρ_c是热解产物的密度，kg·m^-3；是氧气密度。

A_c(min^-1)是热解产物燃烧反应的指前因子，E_c(J/mol)是热解产物燃烧反应活化能，n_c氧气浓度相关系数，参见表3。

3.3.3烟灰模型

根据热解产物燃烧后残留的灰份质量分数为13.01％，灰份的生成方程表示为：

方程(4.9)中，f_ash为热解产物转化为烟灰的质量转化比率系数，f_ash＝0.1301。

3.3.4卷烟纸热解燃烧模型

卷烟纸热解燃烧反应方程可表示为：

其中ρ_z，(kg·m^-3)是卷烟纸组分j的密度；ρ_z(kg·m^-3)是卷烟纸总密度；ρ_z0卷烟纸起始密度，为548.2kg·m^-3。A_z，j，E_z，j，n_z，j取值参见表4。

4.能量方程

4.1有效热导率

在烟丝燃烧过程中，当温度超过1000K时，辐射的影响很大。在固相中，辐射对温度方程的影响用Rosseland近似的方法模拟。在气相中，忽略辐射效应。

气相有效热导率(k_g，eff)为：

固相有效热导率(k_s，eff)为：

k_g是空气的导热系数，k_g＝0.0242w·m^-1·K^-1，是烟丝孔隙率，ε是烟丝辐射系数，k_s是烟丝导热系数，d_p是烟丝孔径，取值见表1。σ是斯忒藩–波尔兹曼常数，σ＝5.67×10^-8W/(m²·K⁴)。

能量方程源项

多孔介质中固相能量方程：

多孔介质中气相能量方程：

方程(5.3)和(5.4)中，

C_p，s(J·kg^-1·K^-1)是烟丝比热容，1043J·kg^-1·K^-1；

C_p，g(J·kg^-1·K^-1)是气体比热容，1004J·kg^-1·K^-1；

ρ_s(kg/m³)是所有固体的总密度，包括热解前驱体、热解产物和烟灰ρ_s＝ρ_v+ρ_c+ρ_z+ρ_ash，它一个变量；

ρ_g(kg/m³)是气体的密度，1.225kg·m^-3；

T_s是烟丝温度，K；

T_g是气相温度，K。

多孔介质比表面方程：

固体能量源项方程：

其中，ΔH_i(kJ/kg)分别为烟丝热解产物燃烧释放的氧化热(ΔΗ_c，17570 kJ·kg^-1)和烟丝中水分的蒸发反应热(ΔΗ_w，-2257kJ·kg^-1)，取值见表1。假定热解反应热很小，可以忽略。

气固间的传热系数h(W.m^-2.K^-1)，可通过下列方程计算得到：

其中μ_g是气体粘度，kg·s^-1·m^-1；v_i是i方向上的速度，m·s^-1；Re是雷诺数；Pr是普朗特数；Nu是努谢尔数。

本模型中卷烟圆柱体表面辐射简化到圆柱固体能量源项计算中，计算公式为：

其中，S_r(J/m³s)为能量源项，σ是斯忒藩–波尔兹曼常数；R'(m)是卷烟烟支圆柱体半径，0.004m；L(m)是卷烟烟支圆柱体长度，0.058m。

5.动量传递方程

6.1当气流从卷烟燃烧锥前端进入

此时，动量方程源项由粘性阻力和惯性阻力组成：

其中，v_i是气体在i方向上的速度，m·s^-1；K是烟丝渗透率，m²；C是惯性项的经验常数。

K与C使用欧根方程求解，

其中，d_p是烟丝孔径，是烟丝孔隙率。

K为烟丝渗透率

在烟丝燃烧的过程中，烟丝渗透率K会发生变化。假定烟丝渗透率随未燃烧烟丝的密度呈线性变化：

K＝K_u(1-g)+K_bg (5.4)

K_u是未燃烧烟丝渗透率，5×10^-10m²，K_b是已燃烧烟丝渗透率10⁵m²。g是一个插入因子。ρ_s是固体的总质量浓度，固体包括热解前驱体、热解产物和烟灰，它是一个变量。ρ_v，0是烟丝的起始密度(见表1)。

6.2当气流从卷烟燃烧线后端进入

此时，动量方程源项组成：

其中K_z是卷烟纸的渗透率，m²。发明人发现，卷烟纸渗透率与温度有关。因此，设定卷烟纸在473-623K温度范围时的渗透率K_z，m＝1.5×10^-9m²。未燃烧卷烟纸渗透率K_z，u＝10^- ¹⁵m²，已燃烧区域的卷烟纸渗透率K_z，b＝10⁵m²。

6.CO和焦油释放量方程

称取1.0g卷烟烟丝样品装入CN104267140A的热解燃烧实验平台。图6为热解燃烧实验平台的示意图。其中：1-气体源，2-质量流量控制器，3-计算机，4-温控系统，5-烟气分析装置，6-剑桥滤片，7-红外灯管，8-烟草样品，9-热电偶，10-保温层，11-腔体，12-石英玻璃管。此处引入CN104267140A全部内容。

烟丝样品装入石英玻璃管中部，并置于快速管式升温炉中，左侧通入不同的气氛，气体流量设定为2.1L/min，通气3min后，控制升温速率为20K/s，升温至目标温度下平衡10min，用剑桥滤片捕集烟气中的焦油，透过剑桥滤片的气相物质中的CO利用烟气分析仪在线监测。每个实验条件重复3次，取平均值。

7.1烟气中CO释放量随温度和氧浓度变化的关系

图7是烟丝分别在2％O₂+98％N₂混合气、10％O₂+90％N₂混合气以及在空气(23％O₂+67％N₂)中快速热解生成的CO总量随温度的变化趋势图，纵坐标为平均每克烟丝热解生成CO的质量。

7.2烟气中焦油的释放量随温度和氧浓度变化的关系

图8示出焦油在不同温度以及不同氧气浓度条件下的释放量。当温度低513K时，焦油释放量随温度的上升迅速增加至90～100mg，然后温度继续升高，焦油释放量变化并不大。

图9展示了不同条件热解燃烧剑桥滤片所收集的焦油状态。总体来看，在三种气氛中剑桥滤片的颜色相似，423K条件下的剑桥滤片是白色，表明该温度下产生的焦油极少；483K条件下的剑桥滤片呈淡黄色，表明焦油开始慢慢增多；而在温度高于513K以后，剑桥滤片呈棕黄色，焦油已经大量生成，升高温度至663K，剑桥滤片呈棕色，继续升高温度，剑桥滤片的颜色差别不大。

7.3CO和焦油释放量方程

根据CO和焦油释放量只与氧气浓度和温度相关。建立牌号1，对单位烟丝所释放的CO和焦油在不同的氧浓度下随温度变化的方程，见表5。将这些数学关系式以UDF的形式输入到Fluent的卷烟燃烧数学模型中，用于预测焦油和CO在不同氧浓度以及不同温度下的生成情况。

表5

因此对应不同温度和氧气浓度下焦油和CO的源项(kg·m^-3·s^-1)可表示为：

其中ρ_s是指所有固体的总密度，包括五种热解前驱体、焦炭和烟灰的总密度，它是一个变量。

7.滤棒对焦油的过滤模型

焦油属于卷烟烟气中的粒相物，表现为气溶胶颗粒。卷烟烟气通过滤棒时，一部分气溶胶颗粒透过滤棒，另一部分气溶胶颗粒被滤棒截留。通过用截留效率e表示滤棒对卷烟烟气气溶胶的收集能力，定义为进入滤棒并被滤棒所截留颗粒的分数，表示如下：

式中m_in和m_out分别表示进入和离开焦油的质量。

气溶胶颗粒沉积于滤棒纤维丝束表面，主要作用机制包括拦截效应E_IN，惯性效应E_IM，扩散效应E_D以及扩散-拦截相互作用E_ID。滤棒单纤维丝束的截留效率e_s可以表示为：

e_s＝E_IN+E_IM+E_D+E_ID (8.2)

卷烟烟气气溶胶在滤棒的截留率e通过风扇模型过滤理论来表达：

其中L_f是滤棒的长度，0.028m；d_f是滤棒纤维丝束的直径，2.51×10^-5m；ω是滤棒内纤维丝束的体积分数，ω的表达式为：

其中D_t是滤棒的总旦数，35000，d_f是滤棒纤维丝束的直径，2.51×10^-5m，D_s是滤棒的单旦数，3.0；C_fiber是滤棒中卷曲纤维的比率，0.17；S_filter是滤棒的横截面积，5.024×10^-5m²。

焦油在卷烟燃烧过程中的源项S_tar(kg·m^-3·s^-1)为：

根据此源项方程，可在Fluent模型中得到滤棒入口面所有网格的焦油流量(kg/s)，根据动量守恒方程，可在Fluent模型中得到滤棒入口面所有网格的气流流量(kg/s)。将二者的比值作为滤棒入口面烟气中焦油的质量分数f_tar(％)。

滤棒不同位置的截留率可表达为：

其中L_x表示滤棒不同位置距离滤棒前端的x轴向长度。

当烟气经过滤嘴后，部分的焦油沉积于滤棒纤维丝束表面，则焦油的质量减少，因此焦油在滤棒中的源项S′_tar(kg·m^-3·s^-1)此时为负值，可表示为：

根据以上方程，在Fluent模型中可计算抽吸前0s滤棒入口面所有网格的焦油释放量(mg)以及抽吸2s后滤棒出口面所有网格的焦油释放量(mg)，并根据方程(8.1)计算滤棒对焦油的截留率e。

8.1拦截效应E_IN

在气流弯曲时，颗粒随着气流的流线恰好进入到距离滤棒纤维表面小于颗粒半径的时候，粒子有可能被收集，这种作用就是拦截作用，拦截效率为：

截留效率取决于无量纲的拦截参数R_L：

其中d_c为卷烟烟气气溶胶颗粒的直径，假设气溶胶颗粒的直径均为4.4×10^-7m，d_f是滤棒纤维丝束的直径，2.51×10^-5m。

K′_u为桑原动力学因子，表示由于其他滤棒纤维丝束的存在造成的气流弯曲的影响，表达式为：

8.2惯性效应E_IM

卷烟烟气通过滤棒时会发生弯曲，具有一定质量的气溶胶颗粒可能不会完全随着弯曲而运动，而由于惯性的作用，气溶胶颗粒可能还来不及转向而撞击到滤棒纤维上而被截留，这种效益称为惯性效应。惯性效应导致的截留率公式：

J＝(29.6-28ω^0.62)R_L ^2.8 (8.10)

Stk是斯托克斯数，惯性效应依赖Stk数，表达式为：

其中ρ_g为气体密度，μ_g气体粘度，v_i是i方向上的气体速度，d_c为卷烟烟气气溶胶颗粒的直径，d_f为滤棒纤维丝束的直径。

8.3扩散效应E_D

卷烟烟气中小粒径的粒子不会完全随气流一起运动，通常情况下还会由于布朗运动而脱离气流。粒子撞击到滤棒纤维表面，使得表面附近粒子浓度减为零。表面附近的浓度差别将促使粒子在表面的扩散沉积，这种效应称为扩散效应。由于扩散作用产生则收集效率为：

其中Pe是贝克莱数，表示单丝纤维截留效率，表达式为：

其中，D_k为颗粒的扩散系数，m²·s^-1；K_B是波兹曼常数，1.38×10^-23J·K^-1，其中T_filter是滤棒的温度，313K。

8.4扩散—拦截相互效应E_ID

扩散和拦截相互作用而产生的截留效率表达式为：

8.卷烟点燃和吸燃程序

将点烟端5毫米区域设为1000K，点燃8秒。

根据GB/T 19609-2004，按ISO模式抽吸卷烟，抽吸容量35ml，每60s抽吸一口，即抽吸2s，阴燃58s、抽吸2s、阴燃58s…如此循环直至烟支燃尽熄灭的过程。

在抽吸的0～2s内，

v表示每口抽吸容量。

将1～8项的参数编写为用户自定义函数(UDF)，加载到Fluent软件中，建立卷烟燃烧(吸燃)的计算流体力学模型，通过运行该模型能够对卷烟燃烧过程的温度场、物质浓度场进行准确地模拟，进而能够获得卷烟燃烧时有害成分释放量。

通过有限体积法对控制方程进行离散，使用分离求解器求解，并选用隐式模式对控制方程进行线化和求解。用一阶迎风格式将控制方程中所有的项进行离散化操作。SIMPLE算法用于压力--速度耦合。收敛性判据设定为10^-3。迭代间隔时间为0.001s。

二、卷烟吸燃模型温度预测结果与实验值的对比

1.热解产物密度场与卷烟纸燃烧线(已燃烧卷烟纸和未燃烧卷烟纸界面)

图10、11、12展示了不同时刻(60s、120s、180s、240s、300s和360s)卷烟燃烧时的热解产物的密度场变化。60s时燃烧锥已形成，随着燃烧的进行，燃烧中心形状慢慢呈锥形，且向后移动，说明卷烟燃烧持续进行。燃烧锥末端与卷烟纸燃烧线的位置相平齐。这说明本实施例的卷烟吸燃计算流体力学模型预测结果准确。

2.抽吸2s内燃烧线移动距离

图13示出在抽吸的2s(180～182s)燃烧锥末端与卷烟纸燃烧线位置变化图。观察实验中卷烟纸燃烧线向后迁移接近3～4mm，通过卷烟燃烧模型中卷烟纸燃烧的位置变化，可以发现模型预测的燃烧线迁移也接近3～4mm，这说明卷烟抽吸2s内模型与实际参与热解燃烧的烟草量是接近的。

从点燃至滤棒前端3mm处(以模型中卷烟纸燃烧至3mm为参考)，卷烟吸燃模型共抽吸了7口，实际卷烟抽吸6.2口，说明卷烟模型的燃烧速度与实际卷烟的燃烧速度基本保持一致。

3.卷烟燃烧温度预测值与实际值

为了检验模型预测的卷烟燃烧温度与实际温度的准确性，本文利用TF-M100烟支燃吸气相温度场分析仪和三孔道吸烟机，将微细热电偶插入烟支中特定位置，并将烟支置于吸烟机持烟器，运行温度采集软件采集烟支阴燃以及吸燃的温度数据。

定位卷烟中心深度距离卷烟前端22mm，24mm，26mm，28mm，30mm，32mm，34mm和36mm卷烟中心的实时温度，在实验过程中，当观测卷烟纸燃烧线在26mm处进行抽吸，只抽吸1口，之前保持卷烟的阴燃状态。通过实验得到26mm位置在抽吸时刻的温度为933K。在卷烟燃烧模型中，按照实验定点监测距离卷烟前端22mm，24mm，26mm，28mm，30mm，32mm，34mm和36mm卷烟中心的实时温度，设置26mm中心位置温度达到933K时进行抽吸，只抽吸1口，之前保持卷烟的阴燃状态。图10展示了0～350s各位置点温度的实验值和模型拟合值。从图14可以看出，实验值和拟合值的温度曲线的趋势保持一致，整体温度偏差较小。

卷烟燃烧时不同位置点的实验温度与模型预测温度的标准均方根误差如下表所示：

位置(mm)	22	24	26	28	30	32	34	36
									标准均方根误差(％)	16.0	17.5	11.0	8.5	13.9	16.0	15.8	11.9

各测量点误差在18％以下，说明预测值与实验值基本相符合。

三、CO和焦油释放量预测结果与实验值

参考上文建立牌号1卷烟吸燃计算流体力学模型的步骤。通过同样的方式建立牌号2和牌号3的卷烟吸燃计算流体力学模型。

计算抽吸每支卷烟后主流烟气中CO和焦油的释放量，并于实验所测的主流烟气中CO和焦油的释放量进行对比，结果见表6。从表6可以看出，模型预测的3个牌号卷烟的抽吸口数与实验值非常接近。

通过计算每口的焦油和CO的释放量，得到总体抽吸7口卷烟模型预测的焦油的释放量为11.4mg/支，卷烟实际焦油释放量为11.2mg/支，相对偏差为1.8％；卷烟模型预测的CO的释放量为14.1mg/支，卷烟实际CO释放量为13.2mg/支，相对偏差为6.4％。

还对牌号2和3的卷烟也按照牌号1卷烟的方法进行了预测，结果如表6所示，相对偏差也很小。

表6

牌号1	口数	CO(mg)	焦油(mg)
				实验值	6.2	13.2	11.4
预测值	6	14.1	11.2
				相对偏差		6.4％	1.8％
牌号2	口数	CO(mg)	焦油(mg)
				实验值	5.7	11.9	11.2
预测值	6	12.8	11.6
				相对偏差		7.5％	3.6％
牌号3	口数	CO(mg)	焦油(mg)
				实验值	6.4	9.8	8.5
预测值	6	9.3	8.4
				相对偏差		5.1％	1.1％
平均相对偏差		8.5％	2.2％

四、焦油截留率预测结果与实验值

计算牌号1卷烟抽吸6口的总体焦油在滤棒进口释放量为21.1mg/支，焦油的滤棒出口释放量为11.4mg/支，焦油的截留率为46.0％。

因为烟碱主要存在与焦油内，因此实验中通过检测烟碱的截留率用于代表焦油的截留率。实验测定的滤棒对烟碱的截留率为44.5％。因此，预测值与实验值的相对偏差为3.4％。

图15为牌号1卷烟滤棒在一次抽吸时(180s～182s)进出口的焦油释放量以及焦油截留率。计算180～182s单口内，焦油在滤棒进口释放量为3.31mg/口，焦油在滤棒出口释放量为1.82mg/口，焦油的截留率为45.0％。计算卷烟抽吸7口的总体焦油在滤棒进口释放量为21.1mg/支，焦油的滤棒出口释放量为11.4mg/支，焦油的截留率为46.0％。

由于实验测定焦油在滤棒的截留率采用称重法，实验误差较大，因为烟碱主要存在与焦油内，本实施例采用烟碱的截留率测定值作为焦油截留率。实验测定的滤棒对烟碱的截留率为44.5％，相对偏差为3.4％

五、卷烟纸渗透率对预测结果的影响

对于卷烟纸在473-623K温度范围内的渗透率，按下表进行了对比实验，并运行模型获得了结果

表7

通过以上对比实验可知，473-623K的卷烟纸渗透率设定为1.5×10^-9m²，建立的烟草燃烧计算流体力学模型更准确，获得的有害成分释放量结果更接近实验测定值。

六、气流速度场

图16、17、18分别示出卷烟在抽吸180s、181s和182s(即一次抽吸的2s内)的气流速度场变化。

在180s，卷烟阴燃时气流速度较小，最大气流速度只有0.00195m·s^-1，当抽吸180～181s时，气流速度增大，且燃烧锥中心区域的的气流速度较大，最大气流速度为1.032m·s^-1，当抽吸181～182s时，气流速度下降，但是靠近卷烟纸区域的气流速度大于燃烧锥中心区域的气流速度，这是由于此时燃烧线后端的卷烟纸变黄区域的渗透率增大，因此气流将从靠近卷烟纸区域进入。

Claims

1.一种建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的方法，包括：

a)建立模拟卷烟燃烧场景的几何模型；

建立烟丝热解反应动力学方程；

建立烟丝热解产物燃烧反应动力学方程；

建立卷烟纸燃烧反应动力学方程；

建立卷烟纸渗透率随温度变化的方程；

2.根据权利要求1所述的方法，建立卷烟纸渗透率随温度变化的方程的方法包括：

设定第一温度的卷烟纸具有第一渗透率；

设定第二温度的卷烟纸具有第二渗透率；

设定第三温度的卷烟纸具有第三渗透率；

优选地，

第一温度为a K以下，第一渗透率为0.5×10^-15～5×10^-15m²；

第二温度为a～b K，第二渗透率为1.5×10^-9～6×10^-9m²；

第三温度为b K以上，第三渗透率为0.5×10⁵～5×10⁵m²；

a＝450～500K，b＝600～650K；

优选地，设定卷烟纸渗透率为单向变化，只增大不减小。

3.根据权利要求1所述的方法，建立卷烟纸燃烧反应动力学方程的方法包括：

-将m个单峰曲线拟合为m个反应j的动力学方程：

各反应j的动力学方程如下：

<mrow> <mfrac> <mrow> <msub> <mi>da</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msub> <mi>A</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>&CenterDot;</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msub> <mi>E</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <msub> <mi>RT</mi> <mi>z</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>a</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>n</mi> <mrow> <mi>z</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </msup> <mo>;</mo> </mrow>

j＝1～m，j和m为正整数；

4.根据权利要求1所述的方法，建立烟丝热解产物燃烧反应动力学方程的方法包括：

建立对应第二氧气浓度的第二烟丝热解产物燃烧反应动力学方程。

5.根据权利要求1所述的方法，建立烟丝热解产物燃烧反应动力学方程的方法包括：

-在s种氧气浓度的含氧气氛下加热烟丝热解产物，采用热重分析法检测烟丝热解产物的质量变化，分别测得k个烟丝热解产物受热时的微分热重曲线；

-将s个微分热重曲线分别拟合为以下反应k的动力学方程：

k为正整数；

s≥k≥2，s和k为正整数。

6.根据权利要求1所述的方法，建立烟丝热解反应动力学方程的方法包括：

-将n个单峰曲线拟合为n个反应i的动力学方程：

<mrow> <mfrac> <mrow> <msub> <mi>da</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> </mrow> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <msub> <mi>A</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <msup> <mi>&beta;</mi> <mrow> <mn>1</mn> <mo>-</mo> <msub> <mi>m</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> </mrow> </msup> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msub> <mi>E</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mrow> <msub> <mi>RT</mi> <mi>v</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>a</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>n</mi> <mrow> <mi>v</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> </msup> <mo>;</mo> </mrow>

i＝1、2...n，n和i均为整数；

a_v，i为反应i的转化率、的反应i的转化率对时间的导数，T_v为加热烟丝样品的温度，β为加热烟丝样品的升温速率，R为理想气体常数，A_v，i为反应i的指前因子，m_v，i为反应i的升温速率修正参数，E_v，i为反应i的活化能，n_v，i为反应i的反应级数；

7.根据权利要求1所述的方法，步骤b)还包括，建立对应不同氧气浓度的有害成分生成量与温度的关系方程，

优选地，建立对应不同氧气浓度的有害成分释放量与温度的关系方程的方法包括：

在不同氧气浓度的气氛下加热烟丝样品；

检测烟丝样品受热时有害成分的释放量与加热温度的关系；

优选地，有害成分为CO或焦油。

8.根据权利要求1所述的方法，步骤b)还包括建立滤棒对焦油的截留率方程的步骤：

优选地，滤棒对焦油的截留率方程包括：

<mrow> <mi>e</mi> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mn>4</mn> <msub> <mi>&omega;e</mi> <mi>s</mi> </msub> <msub> <mi>L</mi> <mi>f</mi> </msub> </mrow> <mrow> <mi>&pi;</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <msub> <mi>d</mi> <mi>f</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

<mrow> <mi>&omega;</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>D</mi> <mi>t</mi> </msub> <msup> <mi>&pi;df</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <msub> <mi>C</mi> <mrow> <mi>f</mi> <mi>i</mi> <mi>b</mi> <mi>e</mi> <mi>r</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mn>4</mn> <msub> <mi>D</mi> <mi>s</mi> </msub> <msub> <mi>S</mi> <mrow> <mi>f</mi> <mi>i</mi> <mi>l</mi> <mi>t</mi> <mi>e</mi> <mi>r</mi> </mrow> </msub> </mrow> </mfrac> <mo>;</mo> </mrow>

e_s＝E_IN+E_IM+E_D+E_ID；

其中，E_IN、E_IM、E_D、E_ID是中间参数；

<mrow> <msubsup> <mi>K</mi> <mi>u</mi> <mo>&prime;</mo> </msubsup> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>l</mi> <mi>n</mi> <mi>&omega;</mi> </mrow> <mn>2</mn> </mfrac> <mo>-</mo> <mfrac> <mn>3</mn> <mn>4</mn> </mfrac> <mo>+</mo> <mi>a</mi> <mo>-</mo> <mfrac> <msup> <mi>&omega;</mi> <mn>2</mn> </msup> <mn>4</mn> </mfrac> <mo>;</mo> </mrow>

<mrow> <msub> <mi>E</mi> <mrow> <mi>I</mi> <mi>M</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>(</mo> <mi>S</mi> <mi>t</mi> <mi>k</mi> <mo>)</mo> <mi>J</mi> </mrow> <mrow> <mn>2</mn> <msup> <msubsup> <mi>K</mi> <mi>u</mi> <mo>&prime;</mo> </msubsup> <mn>2</mn> </msup> </mrow> </mfrac> <mo>;</mo> </mrow>

J＝(29.6-28ω^0.62)R_L ^2.8；

9.根据权利要求1所述的方法，步骤d)还包括，采用计算流体力学软件的求解软件加载卷烟抽吸程序的步骤。

10.根据权利要求1～6任一项所述的方法，步骤a)建立卷烟燃烧几何模型的的方法包括：

优选地，卷烟区域还包括滤棒区域；

优选地，卷烟燃烧场景几何模型是二维模型。

11.根据权利要求1所述的方法，所述物质浓度场包括以下一项或多项：

氧气浓度场；

烟丝热解产物浓度场；

卷烟纸燃烧产物浓度场；

有害成分浓度场。

12.一种预测卷烟燃烧时温度和/或物质浓度的方法，根据权利要求1～11任一项所述的方法建立所述计算流体力学模型，运行所述模型，根据模型模拟结果获得卷烟燃烧时的温度和/或物质浓度。

13.一种用于建立模拟卷烟燃烧时温度场和/或物质浓度场的计算流体力学模型的系统，包括：

物理和/或化学反应方程建立模块，该模块包括

烟丝热解反应动力学方程建立模块；

烟丝热解产物燃烧反应动力学方程建立模块；

卷烟纸燃烧反应动力学方程建立模块；

卷烟纸渗透率随温度变化的方程建立模块；

优选地，所述系统还包括烟丝热解反应热重分析模块，其用于在不含氧气氛中对烟丝进行热重分析，并将数据输出到烟丝热解反应动力学方程建立模块；

优选地，所述系统还包括烟丝热解产物燃烧反应热重分析模块，其用于在含氧气氛中对烟丝热解产物进行热重分析，并加数据输出到烟丝热解产物燃烧反应动力学方程建立模块；

优选地，所述系统还包括卷烟纸燃烧反应热重分析模块，器用于在含氧气氛中对卷烟纸燃烧反应进行热重分析，并将数据输出到卷烟纸燃烧反应动力学方程建立模块；

优选地，所述系统还包括烟丝有害成分释放量检测模块，其用于：在不同氧气浓度的气氛下加热烟丝样品；检测烟丝样品受热时有害成分的释放量与加热温度的关系；并将检测结果输出到烟丝有害成分释放量方程建立模块；

优选地，所述系统包括滤棒对焦油截留率方程建立模块，其被配制为按以下参数和方程建立滤棒对焦油截留率方程：

e_s＝E_IN+E_IM+E_D+E_ID；

其中，E_IN、E_IM、E_D、E_ID是中间参数；

J＝(29.6-28ω^0.62)R_L ^2.8；

<mrow> <msub> <mi>E</mi> <mrow> <mi>I</mi> <mi>D</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mn>1.24</mn> <msup> <msub> <mi>R</mi> <mi>L</mi> </msub> <mfrac> <mn>2</mn> <mn>3</mn> </mfrac> </msup> </mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>K</mi> <mi>u</mi> <mo>&prime;</mo> </msubsup> <mi>P</mi> <mi>e</mi> <mo>)</mo> </mrow> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </msup> </mfrac> <mo>;</mo> </mrow>