Numerical Simulation of the Burning Process in a King-Size Cigarette Based on Experimentally Derived Reaction Kinetics

The cigarette samples were obtained from China Tobacco Fujian Industrial Corporation (Xiamen, Fujian, P.R. China). For 48 h prior to analysis, the samples were conditioned in a chamber at 295 K and relative humidity of 60%.

The cut tobacco from the cigarette sample was pulverized into powder and then sifted through a 100-mesh sieve. 10 mg of tobacco powder were loaded evenly into an open ceramic pan. The pan was then placed on the sample holder of a SAT 449 F3 (Netzsch Feinmahltechnik GmbH, Selb, Germany). In the first stage, the tobacco powder was pyrolyzed under nitrogen atmosphere. The temperature was increased from room temperature to 873 K in pure nitrogen with a steady flow rate of 100 mL·min⁻¹. The heating rates were 300 K·min⁻¹, 400 K·min⁻¹, 500 K·min⁻¹, 600 K·min⁻¹, 700 K·min⁻¹ and 800 K·min⁻¹, respectively. In the second stage, the obtained residual char was cooled to room temperature first, and then the mass loss of char at different oxygen concentrations was obtained by changing the concentration of oxygen in the carrier gas from 1% to 20%. The flow rate was 100 mL·min⁻¹, and the temperature increased from room temperature to 873 K with a heating rate of 10 K·min⁻¹.

The schematic diagram of the pyrolysis and combustion experiment is shown in Figure 1. Approximately 1.0 g of tobacco sample was used. The quartz sample tube was flushed out by the carrier gas for 5 min. The applied carrier gas was 2% O₂ + 98% N₂, 10% O₂ + 90% N₂ and air, respectively. The tobacco sample was heated up from room temperature to the target temperature at the heating rate of 20 K·s⁻¹. When the target temperature was reached, the tobacco sample was held for 5 min. The gas flow rate was maintained at 2.0 L·min⁻¹, which proved high enough to eliminate the smoke from the tube quickly. The “tar” was trapped in a condensing system using a Cambridge filter pad, which was placed on the outlet side of the furnace, while CO was detected online by smoke analyzer (J2KN, rbr, Iserlohn, Germany). Each experiment was performed twice with good repeatability in order to obtain an average of “tar” and CO data.

Figure 1

The gas temperature measurement system is shown in Figure 2. Eight thermocouples were inserted into the center of the cigarette, their locations started at 22 mm from the lighting end of the cigarette. The temperature data at specific locations during the smoking process were collected online by software.

Figure 2

The implementation of the model was carried out in a custom version of FLUENT's structured mesh solver, FLUENT 14.0 (28). The cigarette model was simplified as a two-dimensional model due to the cylindrical axis symmetric structure. The simulation domain was a central plane through the cigarette containing the longitudinal axis of the cigarette, as shown in Figure 3. There are four computational domains. The tobacco rod, cigarette paper and filter rod domains consisted of porous media. In the external environment, the pressure, temperature and gas composition were fixed at ambient conditions. A non-uniform structured mesh of 15854 control volumes was used in the domains, and the meshes near the surfaces were refined. Domain 1) was divided into 2790 grids (99 for the x-axis and 30 for the y-axis) while domains 2) and 3) were divided into 738 and 720 grids separately. A number of important parameters and values related to each domain are shown in Table 2.

Figure 3

Domain 1) was the point source of energy and mass, which included the tobacco pyrolysis and combustion reaction kinetics. The tobacco was consumed and produced char, ash, and smoke, meanwhile releasing energy. Domain 2) affected the resistance to gas flow entering from the cigarette paper. Domain 3) covered the filtration model of the aerosols across the filter. Domain 4) is the external environment. All the exchange of heat, mass and momentum among the four domains were calculated based on the laminar flow model with energy equation and the species transport models with the effect of diffusion.

The initial pressure, temperature and gas composition were also set at ambient conditions within the entire computational domains. Tobacco and cigarette paper were initialized in the unburned state. The outside boundary of the Domain 4) was set as pressure outlet. The mouth end of the cigarette was set as the velocity inlet. The flow velocity was prescribed depending on the International Organization for Standardization smoking regime (ISO 3308) with 35-mL puffs of 2 s duration, and a rate of one puff every 60 s.

To ignite the cigarette, the temperature at the tip of the cigarette was raised to 1000 K, and 35 mL air was drawn for two seconds. After this ignition period, the pre-programmed smoking regime was applied.

Accurate pyrolysis reaction kinetics are necessary to reasonably predict the yield of char, the fuel for the combustion process. The differential thermal gravity (DTG) curve of tobacco pyrolysis at 300 K·min⁻¹ is presented in Figure 4. The pyrolysis DTG curve can be approximated by five Gaussian peaks (R1–R5), so the tobacco is considered as consisting of five precursors. They are moisture, volatiles, hemicellulose, cellulose and lignin, which is consistent with the results reported in (20). According to the area percentages, the mass fractions f_p,j for each pre cursor are listed in Table 3.

Figure 4

The pyrolysis of tobacco can be regarded as the parallel reactions of five precursors (21, 22). The pyrolysis kinetics of each precursor are expressed by Arrhenius equation: [1] ∂αp,j∂t=Ap,jexp(−Ep,jRTs)(1−αp,j)np,j {{\partial {\alpha _{p,j}}} \over {\partial t}} = {A_{p,j}}\exp \left( { - {{{E_{p,j}}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _{p,j}}} \right)^{{n_{p,j}}}} where α_p,j (%) is the conversion ratio of the j^th precursor, t (min) is the time, A_p,j (min⁻¹) is the j^th reaction pre-exponential factor, E_p,j (kJ·mol⁻¹) is the j^th reaction activation energy, T_s (K) is the solid temperature, R (8.314 J·mol⁻¹·K⁻¹) is the universal gas constant, n_p,j is j^th reaction order.

For reactions that are heterogeneous and non-isothermal with β = dT_s / dt, equation [1] can be rewritten as: [2] ∂αp,j∂Ts=Ap,jβexp(−Ep,jRTs)(1−αp,j)np,j {{\partial {\alpha _{p,j}}} \over {\partial {T_s}}} = {{{A_{p,j}}} \over \beta }\exp \left( { - {{{E_{p,j}}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _{p,j}}} \right)^{{n_{p,j}}}} where β (K·min⁻¹) is the heating rate. A wide range of heating rates was applied, but mostly at much lower heating rates (10 K·min⁻¹ – 100 K·min⁻¹) than are actually present in a cigarette during smoldering and puffing. Therefore, the heating rate range was extended to 800 K·min⁻¹ in this study. Figure 5 shows the approximated Gaussian peaks of five precursors at different heating rates. As seen from it, the effects of heating rate on pyrolysis reactions of five precursors are not the same. Hence, the pyrolysis kinetics of each precursor can be expressed by a modified Arrhenius which contains a calibration factor of the heating rate m_p,j (20): [3] ∂αp,j∂Ts=Ap,jβmp,jexp(−Ep,jRTs)(1−αp,j)np,j {{\partial {\alpha _{p,j}}} \over {\partial {T_s}}} = {{{A_{p,j}}} \over {{\beta ^{{m_{p,j}}}}}}\exp \left( { - {{{E_{p,j}}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _{p,j}}} \right)^{{n_{p,j}}}} The conversion ratio of the j^th precursor α_p,j (%) is defined as: [4] αp,j=Wp,j0−Wp,jfp,j(Wp0−Wpf) {\alpha _{p,j}} = {{{W_{p,j0}} - {W_{p,j}}} \over {{f_{p,j}}\left( {{W_{p0}} - {W_{pf}}} \right)}} where the initial mass percentage of the j^th precursor is W_p,j0 = 100%, the initial mass percentage of tobacco is W_p0 = 100%, the final mass percentage of tobacco is W_pf = 30.83%, measured by the TG experiment of tobacco pyrolysis, W_p,j (%) is the mass percentage of the j^th precursor.

Figure 5

Based on equation [4], the mass loss rate of each precursor can be given by: [5] ∂Wp,j∂Ts=−fp,j(Wp0−Wpf)∂αp,j∂Ts=−fp,j(Wp0−Wpf)Ap,jβmp,jexp(−Ep,jRTs)(1−αp,j)np,j {{\partial {W_{p,j}}} \over {\partial {T_s}}} = - {f_{p,j}}\left( {{W_{p0}} - {W_{pf}}} \right){{\partial {\alpha _{p,j}}} \over {\partial {T_s}}} = - {f_{p,j}}\left( {{W_{p0}} - {W_{pf}}} \right){{{A_{p,j}}} \over {{\beta ^{{m_{p,j}}}}}}\exp \left( { - {{{E_{p,j}}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _{p,j}}} \right)^{{n_{p,j}}}} The overall mass loss rate would be the sum of each precursor mass loss rate: [6] ∂Wp∂Ts=∑j=15∂Wp,j∂Ts {{\partial {W_p}} \over {\partial {T_s}}} = \sum\nolimits_{j = 1}^5 {{{\partial {W_{p,j}}} \over {\partial {T_s}}}} where W_p (%) is the mass percentage of tobacco.

Inserting equation [5] into equation [6]: [7] ∂Wp∂Ts=−∑j=15fp,j(Wp0−Wpf)Ap,jβmp,jexp(−Ep,jRTs)(1−αp,j)np,j {{\partial {W_p}} \over {\partial {T_s}}} = - \sum\nolimits_{j = 1}^5 {{f_{p,j}}\left( {{W_{p0}} - {W_{pf}}} \right){{{A_{p,j}}} \over {{\beta ^{{m_{p,j}}}}}}\exp \left( { - {{{E_{p,j}}} \over {R{T_s}}}} \right){{\left( {1 - {\alpha _{p,j}}} \right)}^{{n_{p,j}}}}} In order to determine the kinetic parameters of each precursor, denoting the experimental data by (∂Wp∂Ts)exp {\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)^{exp }} and the calculated data by (∂Wp∂Ts)calc {\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)^{calc}} which can be calculated by equations [7] and equation [8]: [8] αp,ji+1=αp,ji−(∂Wp,ji∂Ts)calc∂Tsfp,j(Wp0−Wpf)(i=1,2,….N) {\alpha _{p,ji + 1}} = {\alpha _{p,ji}} - {{{{\left( {{{\partial {W_{p,ji}}} \over {\partial {T_s}}}} \right)}^{calc}}\partial {T_s}} \over {{f_{p,j}}\left( {{W_{p0}} - {W_{pf}}} \right)}}(i = {1},{2}, \ldots .N) where N is the number of experimental points, the initial conversion ratio of the j^th precursor α_p,j1 (%) is assumed to be 0%.

The unknown kinetic parameters (A_p,j, E_p,j, n_p,j, and m_p,j) are searched at which the determination coefficient (R²) is maximal. [9] R2=1−∑1N((∂Wp∂Ts)exp−(∂Wp∂Ts)calc)2∑1N((∂Wp∂Ts)exp−(∂Wp∂Ts)exp¯)2 {R^2} = 1 - {{\sum\nolimits_1^N {{{\left( {{{\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)}^{{exp}}} - {{\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)}^{calc}}} \right)}^2}} } \over {\sum\nolimits_1^N {{{\left( {{{\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)}^{{exp}}} - \overline {{{\left( {{{\partial {W_p}} \over {\partial {T_s}}}} \right)}^{{exp}}}} } \right)}^2}} }} The corresponding kinetic parameters are listed in Table 3. The experimental DTG curves and the fitted DTG curves of tobacco pyrolysis at different heating rates are shown in Figure 6, and they are in a good agreement.

Figure 6

The char is formed through the pyrolysis of tobacco and consumed in an oxidation reaction. Reasonable char combustion reaction kinetics are critical for the prediction of the heat generation as they are the driving force for the flameless combustion process.

The char combustion kinetics can be expressed by the Arrhenius equation with β = dT_s/dt: [10] ∂αc∂Ts=Acβexp(−EcRTs)(1−αc)nc(ρO2)no=Acβexp(−EcRTs)(1−αc)nc(ρg0WO2)no {{\partial {\alpha _c}} \over {\partial {T_s}}} = {{{A_c}} \over \beta }\exp \left( { - {{{E_c}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _c}} \right)^{{n_c}}}{\left( {{\rho_{{O_2}}}} \right)^{{n_o}}} = {{{A_c}} \over \beta }\exp \left( { - {{{E_c}} \over {R{T_s}}}} \right){\left( {1 - {\alpha _c}} \right)^{{n_c}}}{\left( {{\rho _{g0}}{W_{{O_2}}}} \right)^{{n_o}}} where α_c (%) is the conversion ratio of char, ρ_O² (kg·m⁻³) is the oxygen density in gas mixture, ρ_g0 (kg·m⁻³) is the initial gas density, W_O² (%) is the mass fraction of oxygen, A_c (min⁻¹) is the char oxidation reaction pre-exponential factor, E_c (kJ·mol⁻¹) is the char oxidation reaction activation energy, n_c is the reaction order of char, which is assumed to be 1, n_o is the reaction order of oxygen.

The conversion ratio of char α_c (%) is defined as: [11] αc=Wc0−WcWc0−Wcf {\alpha _c} = {{{W_{c0}} - {W_c}} \over {{W_{c0}} - {W_{cf}}}} where the initial mass percentage of char is W_c0 = 100%, the final mass percentage of char is W_cf = 52.93%, measured by the TG experiment of char combustion, W_c (%) is the mass percentage of char.

Based on equation [11], the mass loss rate of char can be given by: [12] ∂Wc∂Tc=−(Wc0−Wcf)∂αc∂Ts=−(Wc0−Wcf)Acβexp(−EcRTs)(1−αc)(ρg0WO2)no \matrix{{{\partial {W_c}} \over {\partial {T_c}}} &= - \left( {{W_{c0}} - {W_{cf}}} \right){{\partial {\alpha _c}} \over {\partial {T_s}}} \hfill\cr &= - \left( {{W_{c0}} - {W_{cf}}} \right){{{A_c}} \over \beta }\exp \left( { - {{{E_c}} \over {R{T_s}}}} \right)\left( {1 - {\alpha _c}} \right){\left( {{\rho _{g0}}{W_{{O_2}}}} \right)^{{n_o}}}} Inside a burning cigarette, the supply rate of oxygen is the dominant factor in determining the rate of char combustion and heat generation. Therefore, the char combustion kinetics need to be set up at different oxygen mass fractions W_O² (%). In this study, W_O² = 1%, 2%, 3%, 5%, 10%, 15% and 20% were considered.

In order to determine the kinetic parameters of char, denoting the experimental data by (∂Wc∂Ts)exp {\left( {{{\partial {W_c}} \over {\partial {T_s}}}} \right)^{exp }} and the calculated data by (∂Wc∂Ts)calc {\left( {{{\partial {W_c}} \over {\partial {T_s}}}} \right)^{calc}} which can be calculated by equation [12] and equation [13]: [13] αc,i+1=αc,i−(∂Wc,i∂Ts)calc∂TsWc0−Wcf(i=1,2,….N) {\alpha _{c,i + 1}} = {\alpha _{c,i}} - {{{{\left( {{{\partial {W_{c,i}}} \over {\partial {T_s}}}} \right)}^{{\rm{calc}}}}\partial {T_s}} \over {{W_{c0}} - {W_{cf}}}}(i = {1},{2}, \ldots .N) where the initial conversion ratio of char α_c,1 (%) is assumed to be 0%.

The unknown kinetic parameters (A_c, E_c, and n_O² are searched at which the determination coefficient (R²) is maximal. However, the fitting results showing the reaction order of oxygen cannot be the same at different W_O². Therefore, the oxygen mass fraction is divided into three ranges (0% ≤ W_O² ≤ 2%, 2% < W_O² ≤ 10% and 10% < W_O² ≤ 23%) and the corresponding kinetic parameters are listed in Table 4.

The experimental DTG curves and the fitted DTG curves of char combustion at different oxygen concentrations are shown in Figure 7, and they are in a reasonable agreement.

Figure 7

For single fibers, the filtration efficiency was calculated by the filtration mechanism model of aerosols in filters proposed by Du et al. (19), which includes four filtration efficiencies, such as the interception efficiency (E_IN), the inertial impaction efficiency (E_IM), the diffusion impaction efficiency (E_D), and diffusion-inertial impaction efficiency (E_ID).

Interception efficiency (E_IN) can be calculated by: [41] EIN=f(G)2Ku {E_{IN}} = {{f(G)} \over {2{K_u}}} [42] f(G)=11+G−(1+G)+2(1+G)ln(1+G) f(G) = {1 \over {1 + G}} - \left( {1 + G} \right) + 2\left( {1 + G} \right)\ln \left( {1 + G} \right) [43] G=dpdf G = {{{d_p}} \over {{d_f}}} [44] Ku=−lnαf2−34+αf−αf24 {K_u} = - {{\ln {\alpha _f}} \over 2} - {3 \over 4} + {\alpha _f} - {{\alpha _f^2} \over 4} [45] αf=Dtπdf2(1+Cfiber)4DsSfilter {\alpha _f} = {{{D_t}\pi d_f^2\left( {1 + {C_{fiber}}} \right)} \over {4{D_s}{S_{filter}}}} where G is the dimensionless interception parameter, K_u is the Kuwabara hydrodynamic factor, α_f is the solid fraction of filter, d_p (m) is the particle diameter, d_f (m) is the single fiber diameter, D_t is the total denier of filter, D_s is the denier of per single fiber, C_fiber is the crimping ratio of fibers, S_filter (m²) is the cross-sectional area of filter rod. Inertial impaction efficiency (E_IM) can be calculated by: [46] EIM=(Stk)J2Ku2 {E_{IM}} = {{\left( {Stk} \right)J} \over {2K_u^2}} [47] Stk=ρgdp2v18μgdf Stk = {{{\rho _g}d_p^2v} \over {18{\mu _g}{d_f}}} [48] J=(29.6−28αf0.62)G2−27.5G2.8 J = \left( {29.6 - 28\alpha _f^{0.62}} \right){G^2} - 27.5{G^{2.8}} where Stk is the Stokes number, v (m·s⁻¹) is the gas velocity, μ_g (kg·s⁻¹·m⁻¹) is the gas dynamic viscosity, which can be calculated by equation [38].

Diffusion impaction efficiency (E_D) can be calculated by: [49] ED=2.7Pe23 {E_D} = 2.7P{e^{{2 \over 3}}} [50] Pe=dfvDp Pe = {{{d_f}v} \over {{D_p}}} [51] Dp=KBTfilter3πμgdp {D_p} = {{{K_B}{T_{filter}}} \over {3\pi {\mu _g}{d_p}}} where Pe is the Peclet number, D_p (m²·s⁻¹) is the particle diffusion coefficient, the Stefan-Boltzmann constant is K_B = 1.38 × 10⁻²³ J·K⁻¹, T_filter (K) is the temperature of filter. Diffusion-inertial impaction efficiency (E_ID) can be calculated by: [52] EID=1.23G23(KuPe)12 {E_{ID}} = {{1.23{G^{{2 \over 3}}}} \over {{{\left( {{K_u}Pe} \right)}^{{1 \over 2}}}}} The combined single fiber efficiency can be calculated as equation: [53] εs=EIN+EIM+ED+EID {\varepsilon _s} = {E_{IN}} + {E_{IM}} + {E_D} + {E_{ID}} The overall filtration efficiency of a filter can be calculated from the single fiber efficiency by equation [23]: [54] ε=1−exp(−4αfεsLfπ(1−αf)df) \varepsilon = 1 - \exp \left( { - {{4{\alpha _f}{\varepsilon _s}{L_f}} \over {\pi \left( {1 - {\alpha _f}} \right){d_f}}}} \right) where L_f is the filter length, 0.027 m. The overall filtration efficiency at different positions of the filter can be calculated by: [55] dεdLx=4αfεsπ(1−αf)dfexp(−4αfεsLxπ(1−αf)df) {{d\varepsilon } \over {d{L_x}}} = {{4{\alpha _f}{\varepsilon _s}} \over {\pi \left( {1 - {\alpha _f}} \right){d_f}}}\exp \left( { - {{4{\alpha _f}{\varepsilon _s}{L_x}} \over {\pi \left( {1 - {\alpha _f}} \right){d_f}}}} \right) Therefore, the source term of “tar” in the filter, S“tar”′ S_{"tar"}^\prime (kg·m⁻³·min⁻¹) can be written as: [56] S“tar”'=f“tar”ρgvdεdLx=f“tar”ρgv4αfεsπ(1−αf)dfexp(−4αfεsLxπ(1−αf)df) \matrix{S_{"tar"}^\prime &= {f_{"tar"}}{\rho _g}v{{d\varepsilon } \over {d{L_x}}} \hfill\cr &= {f_{"tar"}}{\rho _g}v{{4{\alpha _f}{\varepsilon _s}} \over {\pi \left( {1 - {\alpha _f}} \right){d_f}}}\exp \left( { - {{4{\alpha _f}{\varepsilon _s}{L_x}} \over {\pi \left( {1 - {\alpha _f}} \right){d_f}}}} \right) \hfill\cr} The source term of “tar” in the combustion of cigarette has been established in Section 3.5.1 as: S“tar”=−10−3Y1∂ρs∂tsee [27] {S_{"tar"}} = - {10^{ - 3}}{Y_1}\matrix{{{{\partial {\rho _s}} \over {\partial t}}} \hfill & {{\rm{see\;[27]}}} \hfill \cr } According to the source term equations, the flow of “tar” and air over the entire cross section of the cigarette at the inlet of the filter can be obtained in FLUENT. The ratio of the two is taken as the mass fraction of the “tar” at the inlet of the filter, f_“tar” (%).

The filtration efficiency of filter can be calculated by using the following equation: [57] η=min−moutmin \eta = {{{m_{in}} - {m_{out}}} \over {{m_{in}}}} where m_in (mg) is the “tar” released at the inlet of the filter before puffing, m_out (mg) is the “tar” released at the outlet of the filter after puffing.

The four domains were discretized into structured control volumes over which the conservation equations for mass, momentum, energy and chemical species were discretized. A non-uniform structured mesh of 15854 control volumes was used in the domains, and the meshes near the surfaces were refined. All the exchange of heat, mass and momentum among the four domains were based on the laminar flow model with the energy equation and species transport models with the effect of diffusion. Standard second-order spatial discretization schemes were used for convective operators, with a second-order upwind discretization of the diffusion terms. Based on that, the gaseous emission towards the environment was rapid compared with the emission within the cigarette. The cigarette burning process was assumed to be free of gravity effects which allowed us to use symmetrical conditions, hence considerably shortening the computation time.

The length of the time interval was controlled by the time scale of the reaction. The number of time steps was set to 60,000, with a relative time step size of 0.01 s. Typically 10 iterations per time step were required. The discretization of the unsteady terms was carried out using a second-order implicit scheme.

In FLUENT, the Navier-Stokes equation was solved for conservation of mass, momentum and energy and other scalars such as turbulence and the chemical species using a pressure-based solver. A pressure-based solver requires pressure-velocity coupling, thus the Semi-Implicit Method for Pressure Linked equations (SIMPLE) (24) algorithm was recommended in order to improve the rate of convergence and to reduce simulation time. Basically, the governing equations were coupled to each other in a manner that the solution process required iterations wherein the entire set of governing equations were solved repeatedly until the solution converged.

All equations were solved sequentially and iteratively in keeping with the FLUENT algorithms. The equations described above were incorporated through the user-subroutines available in FLUENT, though some manipulations were not possible through these subroutines and had to be carried out by making changes to the source.

The permeability of cigarette paper plays an important and critical role during the cigarette burning process, especially the permeability of the cigarette paper behind the char line, which determines the resistance to the gas flow entering from the cigarette paper during puffing and consequently the char combustion rate.

The permeability of the cigarette paper increased almost exponentially with the temperature (25). We have reported that the pyrolysis temperature range of cigarette paper was 473 K – 623 K, and the temperature of the cigarette paper behind the char line could reach 623 K (26). Normally, the permeability of cigarette paper is measured at room temperature, but the actual permeability of the cigarette paper behind the char line was difficult to obtain. Hence, in this model, when the temperature of the cigarette paper was < 473 K, the permeability of the unburned cigarette paper was K_up = 5 × 10⁻¹⁵ m², (≈ 30 CU, the thickness of the cigarette paper was about 0.06 mm), when the temperature of the cigarette paper was > 623 K, the permeability of burned cigarette paper was K_bp = 10⁵ m² (18). The cigarette paper behind the char line was set locating at the position where the temperature ranged from 473 K to 623 K, and its permeability was K_p.

In order to obtain the proper K_p in this model, four cases were considered, K_p = 5 × 10⁻¹⁵ m², 10⁻⁹ m², 1.5 × 10⁻⁹ m² and 10⁻³ m², namely cases 1, 2, 3 and 4. Figure 9 shows the char density fields in the four cases during puffing (240 s, 241 s, and 242 s). The development of the coal shape could be checked visually to verify the adequacy of the simulations. When K_p = 5 × 10⁻¹⁵ m², most char was consumed and could not form the combustion coal, resulting from most of the air flowing into the combustion center, leading to the acceleration of the char combustion reaction. In contrast, when K_p = 10⁻³ m², the combustion coal was more obtuse and the shape nearly remained the same, which was due to most of the air entering from the cigarette paper behind the char line, resulting in the deceleration of the char combustion reaction.

Figure 9

These results show the two cases are unrealistic. When K_p = 10⁻⁹ m² and 1.5 × 10⁻⁹ m², the shapes of combustion coals were conical, which was similar to the realistic situation.

The proper value of K_p was determined by matching the numerical and experimental results. Table 6 shows the numerical and experimental yields of “tar” and CO in the mainstream. It can be seen that case 3 has produced a good match with the experimental data. The predicted yields of “tar” and CO were 12.2 mg/cig and 14.6 mg/cig, and the experimental values were 11.2 mg/cig and 13.2 mg/cig, respectively, with relative deviations of 8.9% and 10.6%, respectively. Therefore, K_p = 1.5 × 10⁻⁹ m² was used in the numerical simulations.

In practice, a full comparison of the experimental results with all predictions of the model is difficult. In order to check the validity of the proposed model, five criteria were chosen. These were: puff number, temperatures, flow velocity, filtration efficiency of the filter and the yields of “tar” and CO under different puff intensities.

Figure 10 shows the development of char density with time in case 3. We see that the combustion cone has formed at 60 s. As the combustion progresses, the shape of the combustion center becomes gradually more conical and moves in the smoldering direction, indicating that the cigarette can keep burning.

Figure 10

Figure 11 shows the permeability changes of the cigarette paper at different times in case 3. It can be noticed that the changes in permeability of the cigarette paper were closely synchronized with the changes in char density, which indicates setting K_p = 1.5 × 10⁻⁹ m² when the surface temperature in the range of 473 K – 623 K was reasonable.

Figure 11

From the ignition point to the point located 3 mm from the tipping paper (the length of the tipping paper of the actual cigarette was 31 mm), the experimental puff number was 6.8. This model terminated smoking at 420.5 s when the char line was located at 50.0 mm, the predicted puff number was 7.3. The agreement between the predicted and experimental puff numbers suggests that the burning speed predicted by the cigarette model is basically consistent with the actual cigarette burning speed.

Figure 12 shows the gas temperature fields during puffing (240 s, 241 s, and 242 s) in case 3. It can be seen that the gas temperature increased rapidly, and the maximum gas temperature of the center of the cigarette combustion cone could rise up to 1292 K at 242 s. In order to verify the accuracy of the predicted cigarette combustion temperature, the temperatures at specific positions were measured. Eight thermocouples were inserted into the center of the cigarette, at 22 mm, 24 mm, 26 mm, 28 mm, 30 mm, 32 mm, 34 mm and 36 mm from the lighting end of the cigarette, respectively. After igniting, the cigarette kept smoldering. When the char line migrated at 26 mm, puffing started. The initial gas temperature at the 26-mm position before puffing was 933 K as shown in Figure 13. In the cigarette combustion model, the cigarette also kept smoldering after ignition. When the predicted temperature at 26 mm reached 933 K, a puff was simulated using the model. The real-time temperatures of the center of the cigarette at the corresponding positions were monitored.

Figure 12

Figure 13 compares the predicted gas temperature at each position from 0 s to 350 s with the experimental data. Table 7 shows the standard root mean square error (NRMSE) of the predicted gas temperatures and experimental gas temperatures for eight positions, and the overall deviation was less than 18%. The results show that this model has reproduced the basic features of the temperature distribution, and is also in good agreement with the experimental data quantitatively.

Figure 13

Figure 14 shows the flow velocity fields during puffing (240 s, 241 s, and 242 s) in case 3. It can be seen that the flow velocity increased at the first second, reached the maximum flow velocity at 241 s, and then dropped significantly. From 241 s to 242 s, the flow velocity near the cigarette paper area was higher than in the central region of the combustion cone. A similar tendency was experimentally observed by Li et al. (27). Because of the increase in the permeability of cigarette paper behind the char line, a large amount of flow avoids the center of the combustion cone and enters from the cigarette paper area.

Figure 14

The density fields of “tar” and CO during puffing (240 s, 241 s, and 242 s) in case 3 are shown in Figures 15 and 16. It is clear from the figures that both “tar”- and CO-release amounts increased significantly during puffing. This situation could be explained from the above simulations based on the char density field, temperature field and flow velocity field. The main difference between smoldering and puffing is the increase in the flow velocity (Figure 14), which resulted in an increase of oxygen supply. Therefore the char combustion reaction rate increased while more heat was released and the temperature increased (Figure 12). After a temperature increase, the tobacco pyrolysis would be accelerated too and more tobacco would be involved in the reactions and produce more char for combustion. Thus, the rapid increase in “tar” and CO release can be attributed to the acceleration of the reactions due to the increase in flow velocity and temperature.

Figure 15

Figure 16

Figure 17 shows the release amounts of “tar” at the inlet and outlet of the filter rod and the filtration efficiency during puffing (240 s – 242 s) in case 3. The overall filtration efficiency varied on a degree of about 44.9%. There were two kinds of effects during this period. On the one hand, the flow velocity was high, which lead to a decrease in the filtration efficiency. On the other hand, the concentration of aerosol particles grew higher, which resulted in an increase in the filtration efficiency. Accordingly, these two effects tended to compensate each other and the filtration efficiency remained stable.

Figure 17

The total “tar” amounts released at the inlet and outlet of the filter rod were calculated to be 22.63 mg/cig and 12.19 mg/cig, respectively, and the predicted filtration efficiency was 46.1%. In this study, nicotine retention efficiency was chosen to represent the experimental filtration efficiency of the filter. The experimentally determined filtration efficiency for nicotine was 44.5%.

In order to reinforce the validity and robustness of the model, case 5 was numerically simulated in which just 25 mL air was drawn for 2 s, all other conditions were kept constant. In this case, the model terminated smoking at 421 s when the char line was at 50.0 mm. The predicted puff number was 7.5, and the actual puff number was 7.

Figures 18 and 19 compare the released puff-by-puff amounts of “tar” and CO in the mainstream smoke under different puff intensities. As the puff intensity decreased, “tar” and released CO amounts also decreased. The predicted yields of “tar” and CO were 9.9 mg/cig and 11.6 mg/cig in case 5, and the experimental values were 9.3 mg/cig and 11.0 mg/cig, with the relative deviations of 6.5% and 5.5%, respectively. The numerical results were in good agreement with the experimental ones.

Figure 18

Figure 19