Methods for Technical Parameter Input and Simulation-Based Prediction of Cigar Combustion Time

In numerical simulation research in tobacco science, the simulation of cigar combustion has long been regarded as a highly challenging multiphysics coupling problem. Unlike machine-made cigarettes—relatively uniform, structurally simple paper-wrapped cylinders—a cigar is a complex porous medium system composed of filler tobacco (Filler), binder tobacco (Binder), and wrapper tobacco (Wrapper). This structural heterogeneity gives rise to a high degree of complexity in oxygen diffusion pathways, heat transfer efficiency, and combustion front morphology. For researchers, how to translate these complex physical phenomena into precise mathematical models and ensure that the calculated burn time closely matches laboratory measurements is the most central criterion by which simulation accuracy is judged.

Physical Modeling Foundations in the Porous Media Framework

This section explains the modeling basis for treating a cigar as a porous medium continuum, covering the momentum equation, layered permeability settings, and effective thermal conductivity correction.

To accurately describe cigar combustion, one must first abandon the oversimplified view of the cigar as a "solid block" and instead adopt the Porous Media Continuum Model. Under this framework, the cigar is treated as a composite system of solid particles (tobacco fibers, ash) and pores (filled with air or pyrolysis gases).

In handling the momentum equation, we cannot simply use the Navier–Stokes equation; instead, we must introduce a Darcy–Forchheimer term to describe the pressure drop when fluid passes through the tobacco layer. The expression is:

$$\nabla p = -\frac{\mu}{K} \mathbf{u} - \frac{1}{2} \rho C_F |\mathbf{u}| \mathbf{u}$$ where $K$ is the permeability (Permeability) and $C_F$ is the inertial resistance coefficient. In practical modeling, the setting of permeability is the key factor determining the combustion morphology (whether combustion is conical or tunnel-type). Given the differing characteristics of each part of the cigar, we need to set $K$ in layers: the permeability of the filler layer is usually higher (approximately $10^{-10}$ to $10^{-9} \text{ m}^2$), while the wrapper layer, owing to the density of its leaf fibers, often has a permeability one order of magnitude lower. The heat conduction process also requires correction using the effective thermal conductivity (Effective Thermal Conductivity, $k_{eff}$). Because there is substantial gas-phase convection and radiation inside the tobacco layer, a single solid thermal conductivity would severely underestimate the rate of energy transfer. We usually estimate it using a volume-weighted method: $$k_{eff} = \varepsilon k_g + (1-\varepsilon) k_s$$

where $\varepsilon$ is the porosity, and $k_g$ and $k_s$ are the thermal conductivities of the gas and solid phases, respectively. Near the combustion front, owing to violent temperature fluctuations, this equivalent treatment must be coupled in real time with local porosity evolution.

Schematic of multiphysics-coupled simulation of the cigar combustion process

10⁻¹⁰–10⁻⁹ m²

Typical permeability range of the filler layer

12%

Uniform moisture content setting of the entire cigar in the initial model

18.5% → 3.2%

Prediction error change before and after introducing the moisture evaporation correction

2.26 MJ/kg

Reference value of the latent heat of vaporization of water

Parameter Input Uncertainty: The Coupling Trap of Moisture and Density

This section reveals how the spatial non-uniformity of moisture content and the heat sink effect become a key pitfall of model distortion.

In years of hands-on modeling, I have found that what most often causes model failure is not the complex kinetic equations but the "oversimplification" of basic physical parameter inputs.

In a 2023 study of the combustion characteristics of a certain high-grade Cuban cigar, we encountered a thorny problem: the simulated burn time was about 18.5% faster than the experimental measurement, and the temperature distribution at the combustion front exhibited an unnatural "spike" shape, completely inconsistent with the smooth combustion wavefront actually observed.

After deep retrospective analysis, we found that the root cause lay in the neglect of the spatial non-uniformity of tobacco moisture content (Moisture Content, MC) and the heat sink effect it produces. In the initial model, to simplify the calculation, we set the moisture content of the entire cigar to a constant 12%. However, in actual production, the wrapper of a cigar, being more directly affected by ambient humidity, often exhibits a moisture gradient from the outside inward. More importantly, water evaporation is a huge heat-absorbing process (the latent heat of vaporization of water is approximately $2.26 \text{ MJ/kg}$).

When the moisture content is high, the energy released by pyrolysis reactions is first used for water evaporation, forming an obvious "temperature buffer zone" at the combustion front. In our improved model, we introduced a transient water evaporation model based on the Arrhenius rate equation and adopted the Stochastic Field Method to simulate the spatial non-uniformity of moisture content. By decoupling water evaporation from pyrolysis reactions and introducing a local heat-of-vaporization term, the prediction error quickly dropped from 18.5% to within 3.2%. This lesson taught us that, when dealing with biomass materials such as tobacco, moisture is by no means merely a mass-fraction variable; it is a dynamic thermodynamic regulator.

Mathematical Treatment of Combustion Kinetics and the Moving Boundary

This section explains the combustion kinetics of alternating pyrolysis and oxidation, along with the mathematical treatment of the moving boundary problem and adaptive mesh techniques.

The core of cigar combustion lies in the alternation of pyrolysis (Pyrolysis) and oxidation (Oxidation). A complete combustion cycle usually follows this logical chain: solid tobacco leaves are heated $\rightarrow$ water evaporates $\rightarrow$ organic matter pyrolyzes to produce volatiles $\rightarrow$ volatiles enter the gas phase and undergo oxidation reactions $\rightarrow$ residual char (Char) undergoes surface oxidation under oxygen diffusion control.

When inputting kinetic parameters, we must adopt a multi-step reaction model. A single activation energy cannot describe the complex chemical composition of tobacco. We usually adopt a multicomponent pyrolysis model, dividing tobacco into cellulose-like, lignin-like, and residual ash components, and assign different Arrhenius parameters ($A, E_a$) to each.

The core algorithm for predicting burn time is essentially a solution to the "Moving Boundary Problem," mathematically similar to the Stefan Problem. The position of the combustion front $\mathbf{s}(t)$ moves over time, and its moving velocity $\dot{r}$ depends on local mass transfer and energy balance:

$$\rho_s (1-\varepsilon) \frac{d\mathbf{s}}{dt} = \dot{m}_{pyro} + \dot{m}_{ox}$$ where $\dot{m}_{pyro}$ is the pyrolysis rate and $\dot{m}_{ox}$ is the char oxidation rate. In actual CFD solution processes, to capture this extremely thin and intense reaction layer, we must employ Adaptive Mesh Refinement (AMR). If the mesh resolution is insufficient to cover the gradient at the combustion front, numerical diffusion artificially "flattens" the reaction rates, causing the predicted burn time to be too long. When handling such highly nonlinear problems, the selection of the time step is also an art—it must satisfy the Courant Number constraint while also accounting for the time scale of thermal diffusion; otherwise, the computation easily diverges at the moment of combustion burst. ## Integration and Engineering Implementation of the Prediction Algorithm To finally obtain an accurate burn time, we established an integral model spanning from micro-scale reactions to macro-scale time history. First, through transient CFD simulation, we obtain the temperature field $T(x,t)$, oxygen concentration field $C_{O_2}(x,t)$, and local moisture content field $M(x,t)$ at each micro-element location at different times. Then, based on these field variables, using the local combustion rate $\dot{r}(x,t)$, we perform spatiotemporal integration over the geometric length $L$ of the cigar. The formula for predicting the total burn time $t_b$ can be simplified as: $$t_b = \int_0^L \frac{1}{\bar{\dot{r}}(x)} dx$$

where $\bar{\dot{r}}(x)$ is the average linear velocity along the axial direction. However, in engineering practice, because cigar combustion exhibits a pulsed "ignition–extinction–reignition" characteristic (especially when simulating the puffing process), this integral process must be accumulated in the time dimension. We need to precisely splice, in time sequence, the combustion progress within each puff cycle (Puff) and the smoldering (Smoldering) stage between two puffs.

When solving the prediction of the "smoldering stage," another key point is the setting of boundary conditions. After puffing stops, the cigar is in a state of slow oxidation in a low-oxygen environment. If the ambient boundary condition is simply set to normal-temperature air, the simulated smoldering temperature will be too high, causing the predicted burn time to be too short. We need to introduce a dynamic boundary that accounts for the ambient convective heat transfer coefficient $h$, as well as the radiative heat loss from the cigar surface.

Summary and Outlook

This section describes the engineering implementation path from transient CFD field variables to integral prediction of burn time.

Through the complex parameter inputs and algorithm construction described above, we can not only achieve accurate prediction of cigar burn time but, more importantly, observe physical details during combustion through simulation—such as the "uneven combustion" phenomenon caused by excessively low wrapper leaf permeability, or the "combustion stall" phenomenon caused by moisture gradients.

Current simulation technology is moving toward higher dimensions. With the popularization of micro-CT (Micro-CT) scanning technology, we can obtain extremely precise microscopic pore structure data of tobacco leaves and import it directly into CFD models, achieving a leap from "continuum approximation" to "real geometric structure." At the same time, coupling machine learning algorithms with traditional CFD—using known experimental data to train neural networks to predict complex reaction rate parameters—may well be the next technological breakthrough in solving the difficulties of obtaining tobacco combustion kinetic parameters and their high uncertainty.