电子烟多口抽吸过程的传热传质数值模拟

黄忠辉，黄洁洁，高峄涵，李 典，高乃平，陆 漓，王萍娟，巫 毅

1. 广西中烟工业有限责任公司技术中心，广西南宁市北湖南路28 号 530001

2. 同济大学机械与能源工程学院，上海市四平路1239 号 200092

3. 上海新型烟草制品研究院有限公司，上海市浦东新区秀浦路3733 号 201315

4. 广西中烟工业有限责任公司柳州卷烟厂，广西柳州市静兰路1 号 545026

电子烟作为一种新型烟草制品，在全球范围内受到广泛关注。常见的电子烟主要由雾化芯、加热元件、储液仓、微电子传感器和电池等元件构成。多孔雾化芯是电子烟的关键部件，主要起到承载和输运烟液的作用，市面上最常见的雾化芯材料为有机棉或多孔陶瓷［1-2］。在电子烟工作过程中，雾化芯孔隙中的烟液吸收电热丝散发出来的热量而气化逸散，同时受气液界面不同相压力差的驱动，存在烟液自补充的过程以供继续受热气化。电子烟的烟液是由丙二醇（PG）、甘油（VG）、香精以及烟碱组成的多元非共沸混合溶液，其中，丙二醇和甘油在常温下不混溶，其二者在烟液质量中占比90%以上，为电子烟烟液的主要组成部分［3］。电子烟的雾化工作过程本质上是双组分或多组分的混合溶液在多孔介质内的受热相变及毛细作用下的液体传输的过程，涉及了液气相变和多孔介质毛细作用下的传热传质耦合问题，当雾化芯孔隙内的烟液被气化消耗时，由于孔隙内由液体单相变为液气双相，液气界面产生压力差，孔隙内液气界面上的毛细作用会使得与雾化芯末端相连接的储液仓里的新烟液被不断输运补充至雾化芯孔隙内。

近年来，电子烟的研究大多是测量分析方面［4-10］，如谢国勇等［5］测量了电子烟雾化剂雾化过程的温度变化特性；Lee 等［8］测量并分析了电子烟的烟液成分、烟液消耗量和消耗率；Chen等［9］测量了电阻丝加热温度、加热功率以及烟液组成等因素对电子烟雾化效果的影响；李典等［10］依据相似原理搭建了便于布置多个温度测点的电子烟扩尺试验台，用于测量电子烟内部不同位置的温度变化特性。而在理论研究方面，Shihadeh等［11］采用烟碱通量（Nicotine Flux）来描述使用者吸入的烟碱剂量和吸入速率。Talih等［12］建立了非稳态零维集总模型来研究加热功率、抽吸模式、烟液组成、烟碱含量等对烟碱递送的影响。Gao等［13］用数值模拟方法研究了电子烟多孔芯毛细管蒸发行为。目前针对电子烟的多孔介质液气相变传热传质机理的理论研究较少。因此，本研究拟通过建立计算机模型，研究电子烟多口抽吸情景下，电子烟多孔雾化芯内烟液气化相变过程的传热传质特性、烟液组成变化及系统的热效率，为电子烟的雾化机理研究提供方法参考。

1 数值模型

1.1 模型介绍

为了研究烟液在电子烟多孔雾化芯内发生的受热雾化过程，选取外绕电加热丝的柱形多孔雾化芯为研究对象，多孔雾化芯两端浸入储液仓中，另有柱形气流通道与雾化芯中央垂直相交。基于所处区域的外围是否有电加热丝缠绕引入的面热源，将雾化芯分为两段：一段是有电加热丝外绕的加热段，另一段是没有电加热丝外绕的输运段。整根雾化芯均受毛细作用影响，但由于气流通道仅和加热段区域重合，因此抽吸气流仅影响加热段。如图1所示，将雾化芯简化视为柱形多孔体，其中红色区域代表加热段，Lh为多孔雾化芯加热段的长度；蓝色区域代表输运段，Lt为多孔雾化芯输运段一端的长度；D 为柱形多孔雾化芯的直径；气流通道视为内径与加热段长度相同的柱形通道。

图1 多孔介质雾化芯结构示意图Fig.1 Structure of a porous atomizer

在电子烟工作过程中，雾化芯孔隙内不断发生液气相变，整个过程中烟液气化温度、气化速率、烟液组分等参数不断变化。因此，为了保证分析计算的可操作性和便利性，在数值模型建立过程中采用如下假设：①多孔介质视为各向同性的刚性介质，内部孔隙均与外界连通，孔隙分布均匀，连续介质假设和局部热平衡假设适用；②烟液作为不可压缩流体，烟液的气化产物为理想气体，烟液是丙二醇和甘油的混合物，未考虑烟碱等占比较少的组分；③电加热丝的加热功率全部传递给多孔雾化芯及烟液，不考虑其他损耗；④烟液各组分为常物性参数，混合烟液的物性参数随组分占比变化而变化；⑤在数值模拟中，整个过程经过的时间划分为多个时间步，单个时间步内各个物理量保持不变，单个时间步内烟液的流动为稳定的轴向一维流动，烟液的传输仅受毛细压力驱动，忽略接触角的影响；⑥不考虑多孔介质和烟液的热量传递时间，忽略多孔介质和储液仓接触端的烟液入口效应。

1.2 控制方程

从能量守恒来看，将含液多孔雾化芯视为系统，外界输入系统的能量为电热丝输入的热量，在抽吸过程中，这些能量将分别用于含液多孔介质的升温、烟液的气化、含液多孔介质的轴向导热、含液多孔介质与抽吸气流的对流换热以及对外辐射。在抽吸间隙，气流通道中空气流速较小，无法触发气流开关，与气流开关相连的电加热丝停止工作，电加热丝不提供能量，而其他能量消耗形式依旧存在。每个能量项都可以通过非稳态能量守恒方程来描述和计算。雾化芯加热段的能量守恒方程如下：

式中：hD，i是双组分烟液中组分 i的对流传质系数，m/s；Δρg，i是气相中组分i 在多孔介质表面富集的密度与在抽吸气流中的密度之差，kg/m3；Tb是双组分烟液的沸点温度，℃。

在蒸发浓度较低、扩散通量较小的情况下，对流传质的准则关联式和对流传热具有同一函数形式，即可以用圆柱绕流下的对流传热关联式Nu=f（Re，Pr）以相同形式类比到蒸发传质关联式Sh=f（Re，Sc）上［14］。圆柱绕流下的对流传热关联式采用Churchill-Bernstein经验式（适用于Rea·Pra＞0.2 的情形）［14］，故有蒸发传质关联式：

式中：Di是双组分烟液中组分i 在空气中的扩散系数，m2/s；Shi是混合溶液中组分 i 蒸发传质的舍伍德数，表征对流传质与扩散传质的强度相对大小；Rea=uaD/va是气流通道内的空气雷诺数；ua是抽吸气流流速，m/s；va是空气运动黏度，m2/s；Sci=va/Di是蒸发传质的施密特准则数，表征速度边界层和浓度边界层的厚度相对关系。

当抽吸气流视为纯净空气时，式（7）中Δρg，i即等于气相中组分 i 的密度，Δρg，i＝ρg，i。由克拉贝龙方程ρg，i＝PiMi（/RT），Mi是组分i 的相对分子质量，kg/mol；R 是通用气体常数，J（/mol·℃）；Pi=xiPi*是双组分烟液中组分i的气相分压，Pa；其中xi是组分i在烟液中的摩尔分数；Pi*是烟液中组分i 的饱和蒸气压，Pa，可以由安托万方程计算得到：

式中：Ai、Bi、Ci分别是烟液中组分i的安托万常数，取值参考美国国家标准与技术局化学数据库（NIST Chemistry WebBook）［15-16］，如表1所示。

表1 所研究双组分烟液的安托万方程系数Tab.1 Antoine equation coefficients of the studied e-liquid with two components

抽吸间歇电加热丝不提供能量，但因含液多孔雾化芯温度仍高于烟液蒸发温度，表面烟液仍在发生气化。由于抽吸间歇气流通道中无抽吸气流，从多孔雾化芯表面气化的丙二醇、甘油蒸汽会逐渐富集在雾化芯表面上方的空间中，使得表面附近空间内的蒸汽压增大，从而抑制多孔介质表面的烟液气化。为了定量描述这一部分影响，此处引入间歇气化系数clapse来对抽吸间歇期间的烟液气化传质进行修正：

抽吸期间，烟液的流动及雾化实质是液体在多孔介质内在毛细作用的输运驱动下受热雾化的过程，该输运过程的控制方程包括质量守恒方程、动量方程、毛细压力—饱和度方程。取右侧一端输运段轴向位置x 处一段长度为dx 的微单元进行分析，根据多孔介质基本特性参数的定义，孔隙率ε为孔隙体积与多孔介质总体积的比值，饱和度s为孔隙内液体体积与孔隙体积的比值，因此饱和度可表示为：

式中：h（x）为轴向位置x处液相连续体的液面高度，m。

对微单元dx 进行连续性分析，则质量守恒方程为：

式中：h1为微单元dx 出口处的液相连续体的液面高度，m；ρl为液体密度，kg/m3；u1为微单元dx出口处的烟液在毛细作用下发生流动的速度，m/s；h2为微单元dx入口处的液相连续体的液面高度，h2=h（x+dx），m；u2为微单元dx入口处的烟液在毛细作用下发生流动的速度，m/s。

非饱和多孔介质中液体的毛细流动满足达西定律［17］即流动驱动力为液体流动方向的前后压差。根据毛细压力的定义，毛细压力Pc实质是气液界面两端不同相的压力差，即Pc=Pg-Pl，其Pg为气相压力，Pa；Pl是孔隙内的液体压力，Pa。由于气流通道与外界连通，气相压力等同视为外界气压Pa。毛细压力是含有液气相变多孔介质中流体输运的主要驱动力，因此烟液在多孔介质内流动的动量方程可化为：

式中：K是多孔介质的绝对渗透率，m2；μ1是液体动力黏度，Pa·s。

渗流领域广泛应用Kozeny-Carman 方程作为预测多孔介质渗透率的半经验公式。基于邓永锋等［18］的汇总和Carman［19］的研究，多孔介质绝对渗透率可表示为：

式中：rave为平均孔隙半径，μm；τ为孔道弯曲因子。

采用Leverett［20］提出的毛细压力通过无因次函数J（s）与饱和度相关联，则毛细压力—饱和度方程表示为：

式中：σ是气相与液相界面处的气液表面张力系数，N/m。

忽略残余饱和度，并考虑数值求解的可行性，无因次函数采用 Udell［21］提出并得到广泛使用［17，22-23］的关联经验式：

2 模拟条件及求解方法

2.1 边界条件

模型的边界条件为第一类边界条件，即给定起始时刻的各物理量初值，随后逐步求解各个时刻的物理量。对于雾化芯内的烟液传输，有下列边界条件需要满足：

（1）由于雾化芯两端一直浸没在储液仓中，则多孔雾化芯的烟液传输入口处：

（2）多孔雾化芯输运段与加热段的交界处，毛细作用下的烟液传输质量速率mc满足：

2.2 求解方法

数值模型的求解采用MATLAB 作为计算工具，其中液气相变过程的时间求解采用一阶前向差分格式，烟液传输的空间分布求解采用二阶中心差分格式。基于给定的初始条件（初始温度、初始烟液组成、加热功率等），确定初始的物性参数后，对控制方程采用有限差分法进行离散，依次求解各个物理量，数值模拟过程的各组分参数如表2所示。

表2 烟液组分参数Tab.2 Parameters of e-liquid components

3 结果与讨论

3.1 模型校验

为了验证模型的准确性，将数值模拟条件设置为与文献［6］所述实验工况一致，以有机棉作为多孔雾化芯材料，雾化芯直径3.8 mm，加热段长度8 mm，输运段一端长度4 mm，多孔介质孔隙率44%，环境温度22 ℃，初始烟液组分比例为PG∶VG=100∶0［体积分数占比（%），下同］，加热功率为15 W，每口抽吸时长为3 s，两口抽吸之间间隔30 s，抽吸口数为7 口。引入的间歇气化系数分别取值为0、0.2、0.4、0.6、0.8、1.0进行模拟，将不同间歇气化系数取值下烟液温度随时间变化的模拟结果与文献［6］的结果进行对比，结果见图2。由图2 可以看出，不同间歇气化系数下的模拟结果均与文献数据结果一致，即在每一口抽吸期间，烟液温度迅速上升，而在电加热丝停止工作的抽吸间歇期间，烟液温度逐渐由快速下降变为缓慢下降。对比6种不同的间歇气化系数取值的模拟结果，当间隙气化系数为0.2 时，模拟结果与文献［6］的测量结果最接近，其相对平均误差为10.2%，认为模拟结果较准确，将以此间隙气化系数取值开展后续研究。

图2 不同间隙气化系数下的烟液温度变化Fig.2 Temperature variation of e-liquid at different interval vaporization coefficients

3.2 液气相变传质过程模拟结果

3.2.1 逐口抽吸开始、结束时刻烟液温度变化

为了研究抽吸口数对多孔介质液气相变传热传质过程的影响，用已验证的多口抽吸模型模拟了总抽吸口数为5口的多口抽吸过程。设定多孔介质孔隙率50%，雾化芯直径3 mm，加热功率9.09 W，加热段长度8 mm，输运段一端长度4 mm，初始烟液组成为PG∶VG = 50∶50，每口抽吸时长为3 s，两口抽吸之间间隔为30 s，环境温度为10 ℃。逐口抽吸始末烟液温度随抽吸口数变化的模拟结果如图3所示。

图3 逐口抽吸始末烟液温度Fig.3 Temperature of e-liquid at the beginning and end of each puff for 5 consecutive puffs

由图3 可以看出，随着抽吸口数增加，每口抽吸开始时刻和抽吸结束时刻的烟液温度在逐渐上升，从第2口开始，每口抽吸结束时烟液温度比上一口平均提高4.4%。从多口抽吸的角度分析，每口抽吸结束时的烟液状态除了与3 s 抽吸期间的多孔介质液气相变传热传质过程有关，还与抽吸开始时烟液的状态有关；由于前一口抽吸结束后，有余温的烟液在抽吸间隙继续进行传热传质。雾化芯孔隙中的烟液在抽吸间隙期间会发生降温，这是由于抽吸间隙电加热丝未向多孔介质输入能量，而含液多孔介质与外界环境依然存在温差，导热、辐射、对流损失仍然发生。降温的幅度除了与上述3种热损失有关以外，还与抽吸间隙的时长有关。在环境温度条件下，抽吸间隙30 s 时间内烟液的降温幅度有限，烟液无法降至环境温度，从而使后续抽吸始末的烟液温度相继提高。

3.2.2 逐口抽吸烟液PG、VG组分占比变化

为了进一步考察雾化过程中烟液组成的变化，对不同抽吸口序下的烟液气化速率以及每口抽吸结束时各组分所占体积分数进行分析，结果见图4和图5。针对PG∶VG = 50∶50的烟液在不同雾化温度下的释放实验［24］表明，当雾化温度为260 ℃时，PG、VG的释放量接近一致；当雾化温度高于260 ℃时，VG相对释放量高于PG。模拟结果与上述实验结论显示出一致的规律，当第4 口抽吸结束时，由图3 可见烟液的平均温度为259.3 ℃，同时在图4 中第4 口抽吸时VG 组分的平均气化速率略小于PG 组分；在第5 口抽吸结束时，由图3 可见烟液的平均温度为263.3 ℃，同时在图 4 中第 4 口抽吸时 VG 组分的平均气化速率已超过PG组分。此外，烟液整体的每口平均气化速率与烟气中体积分数占比较大的PG 组分变化趋势一致，呈先升高后下降的趋势，气化较少的VG 组分每口抽吸平均气化速率则是随着口序增大而递增。第2口抽吸时的平均气化速率最大，达到了2.8 mg/s，比第1 口提高了12.2%，从第3 口开始，每一口抽吸时的烟液平均气化速率比上一口平均降低4.3%。虽然每口抽吸始末的烟液温度不断上升，但第2口抽吸结束时，多孔介质空隙内烟液的PG组分体积分数已经从初始的50%降到21%，由于VG组分具有相对更高的汽化潜热值，因此随着烟液中PG组分的减少，将导致双组分混合烟液的汽化潜热值增大，即相同的能量耗费下一定时间内的气化产物质量将减少，导致烟液整体平均气化速率降低。

图4 逐口抽吸平均气化速率Fig.4 Average vaporization rate of each puff for 5 consecutive puffs

图5 逐口抽吸烟液组分变化量Fig.5 Component variation of e-liquid released by each puff for 5 consecutive puffs

从烟气组分来看，图5中每口烟气中VG组分的占比不断增大，PG 组分占比均≥50%。从烟液组分来看，雾化芯孔隙中液相PG、VG的逐口变化规律呈现与气相组分相反的趋势。Talih等［12］的同类模拟结果也指出：在单口抽吸的早期阶段，挥发性更强的PG 以远高于挥发性较弱的VG 的速率从烟液中排出，因此气相组分几乎完全由PG组成。在小功率下（2 W），几乎仅有PG进行蒸发，抽吸3 s后气相中PG质量分数达95.6%；在大功率下（11 W），随着抽吸的进行，烟液温度继续升高，蒸发速率加快，雾化芯中残余烟液中VG组分比例的不断升高使得更多的VG向气相转移，导致蒸气成分更快地接近母液，约1.25 s后气化达到稳态。本研究中的多口抽吸过程可视为单口抽吸过程的组合，但进一步考虑了抽吸间隙对下一口抽吸的烟液初始状态的影响，本研究中的加热功率条件（9.09 W）介于Talih 模拟的两种功率之间，由此可判断在逐口抽吸开始时，PG 组分先进行快速蒸发，且达到稳定时间晚于1.25 s，导致气相中PG 比例大于VG；但随着每口抽吸结束时刻烟液温度不断上升，达到稳定气化的时间会缩短，使得气相中PG、VG 的体积比更接近储油仓中两种组分的比例。对于液相而言，由于在前4口抽吸中烟液平均温度未达到260 ℃，VG组分相对于PG组分气化更慢，VG 逐渐富集在孔隙内，表现为VG 组分在液相中的体积分数不断上升。

3.2.3 逐口抽吸热能耗及热效率

不同抽吸口序下的热消耗及热效率如图6所示，从图中可以看出，烟液在多孔介质孔隙内的液气相变过程中，多孔雾化芯和烟液组成的系统，其每口抽吸的系统升温耗能随抽吸口序的增大呈现先减小后增大的趋势，对于总抽吸口数为5 口的情况，第1～2口抽吸系统升温耗能从19.02 W降低至17.78 W，第3～5口升温耗能从17.94 W升高至18.25 W。每口抽吸的烟液气化耗能则随抽吸口数变化先增大后减小，与每口抽吸平均气化速率的变化趋势一致。轴向导热损失、空气对流损失和对外辐射损失均随抽吸口序的增大而增大，这是由于导热、对流、辐射3项热损失与该口抽吸结束时的烟液温度正相关。每口抽吸的轴向导热损失、空气对流损失和对外辐射损失分别平均相对提高8.4%、6.8%、11.6%。不同抽吸口数下，每口抽吸热效率的范围为90.9%～93.3%。每增加一口抽吸口数，每口抽吸的热效率平均减小0.6%。

图6 逐口抽吸热消耗及热效率变化量Fig.6 Variations in heat consumption and heat efficiency of each puff for 5 consecutive puffs

4 结论

利用计算机软件MATLAB 建立了电子烟烟液在多孔雾化芯内液气相变和毛细作用下的液体传输的数值模型，模拟了多口抽吸下的烟油型电子烟烟液相变的传热传质过程，研究了多口抽吸条件下，电子烟烟液相变过程的传热传质特性、烟液组成变化以及热利用情况。结果表明：①烟液在抽吸间隙发生降温但无法降至环境温度，随着抽吸口数的增加，逐口抽吸开始时刻和结束时刻的烟液温度较上一口相比平均提高4.4%。②烟液的平均气化速率随抽吸口序的增大而呈现先增后减的趋势，第2口抽吸时达到最大值2.80 mg/s，比第1 口抽吸时相对提高12.2%。③多口抽吸时，随抽吸口序的增加，逐口的系统升温耗能呈现先减小后增大的趋势，逐口的烟液气化耗能则先增大后减小；逐口的导热、对流、辐射3项热损失逐级增大。每口抽吸的热效率随抽吸口数增加而减小，当抽吸口数从1 口增加至5 口，每口抽吸的热效率从93.3%降至90.9%。