Методы ввода технических параметров и прогнозирования времени горения в компьютерном моделировании процесса горения сигары

В исследованиях численного моделирования в табачной науке моделирование процесса горения сигары всегда считалось чрезвычайно сложной сопряжённой мультифизической задачей. В отличие от сигарет — относительно однородных бумажных обёрток простой структуры, сигара представляет собой сложную систему пористой среды, состоящую из набивочного табака (Filler), связующего табака (Binder) и покровного табака (Wrapper). Эта структурная неоднородность приводит к высокой степени сложности путей диффузии кислорода, эффективности теплопередачи и морфологии фронта горения. Для исследователей то, как преобразовать эти сложные физические явления в точные математические модели и обеспечить высокое соответствие рассчитанного времени горения лабораторным измерениям, является наиболее центральным критерием оценки точности моделирования.

Основы физического моделирования в рамках пористой среды

В этом разделе излагаются основы моделирования сигары как континуума пористой среды, охватывающие уравнение импульса, послойное задание проницаемости и коррекцию эффективной теплопроводности.

Чтобы точно описать горение сигары, прежде всего необходимо отказаться от упрощённого представления о сигаре как о «твёрдом блоке» и принять континуальную модель пористой среды (Porous Media Continuum Model). В рамках этой модели сигара рассматривается как композитная система, состоящая из твёрдых частиц (волокон табака, золы) и пор (заполненных воздухом или газами пиролиза).

При обработке уравнения импульса мы не можем просто использовать уравнение Навье–Стокса (Navier-Stokes); необходимо ввести член Дарси–Форхгеймера (Darcy-Forchheimer) для описания перепада давления при прохождении жидкости через слой табака. Выражение имеет вид:

$$\nabla p = -\frac{\mu}{K} \mathbf{u} - \frac{1}{2} \rho C_F |\mathbf{u}| \mathbf{u}$$ где $K$ — проницаемость (Permeability), а $C_F$ — коэффициент инерционного сопротивления. В практическом моделировании задание проницаемости является ключевым фактором, определяющим морфологию горения (коническое горение или туннельное горение). Учитывая особенности различных частей сигары, нам необходимо задавать $K$ послойно: проницаемость набивочного слоя обычно выше (примерно $10^{-10}$ до $10^{-9} \text{ м}^2$), тогда как проницаемость покровного слоя из-за плотности листовых волокон часто на порядок ниже. Процесс теплопроводности также требует коррекции с помощью эффективной теплопроводности (Effective Thermal Conductivity, $k_{eff}$). Поскольку внутри слоя табака существует значительная конвекция и излучение в газовой фазе, одна лишь теплопроводность твёрдой фазы серьёзно занизила бы скорость переноса энергии. Обычно мы оцениваем её методом объёмного взвешивания: $$k_{eff} = \varepsilon k_g + (1-\varepsilon) k_s$$

где $\varepsilon$ — пористость, а $k_g$ и $k_s$ — теплопроводности газовой и твёрдой фаз соответственно. Вблизи фронта горения из-за резких колебаний температуры эта эквивалентная обработка должна быть сопряжена в реальном времени с локальной эволюцией пористости.

Схема сопряжённого мультифизического моделирования процесса горения сигары

10⁻¹⁰–10⁻⁹ m²

Типичный диапазон проницаемости набивочного слоя

12%

Единое задание влагосодержания всей сигары в первоначальной модели

18.5% → 3.2%

Изменение ошибки прогноза до и после введения поправки на испарение воды

2.26 MJ/kg

Справочное значение скрытой теплоты парообразования воды

Неопределённость ввода параметров: ловушка сопряжения влажности и плотности

В этом разделе показано, как пространственная неоднородность влагосодержания и эффект тепловой ловушки становятся ключевой причиной искажения модели.

За годы практического моделирования я обнаружил, что к отказу модели чаще всего приводит не сложность кинетических уравнений, а «чрезмерное упрощение» ввода базовых физических параметров.

В исследовании 2023 года, посвящённом характеристикам горения одной высококлассной кубинской сигары, мы столкнулись с трудной проблемой: рассчитанное время горения оказалось примерно на 18.5% меньше экспериментального значения, а распределение температуры на фронте горения имело неестественную «пиковую» форму, полностью не соответствующую наблюдавшемуся в реальности плавному волновому фронту горения.

После углублённого ретроспективного анализа мы выяснили, что корень проблемы кроется в игнорировании пространственной неоднородности влагосодержания табака (Moisture Content, MC) и вызванного им эффекта тепловой ловушки. В первоначальной модели для упрощения расчёта мы задали влагосодержание всей сигары постоянным — 12%. Однако в реальном производстве покровный слой сигары, испытывающий более прямое влияние влажности окружающей среды, часто имеет градиентное распределение влагосодержания снаружи внутрь. Более того, испарение воды — это огромный эндотермический процесс (скрытая теплота парообразования воды составляет примерно $2.26 \text{ МДж/кг}$).

При высоком влагосодержании энергия, выделяемая реакциями пиролиза, сначала расходуется на испарение воды, что формирует на фронте горения отчётливую «температурную буферную зону». В усовершенствованной модели мы ввели модель переходного испарения воды на основе уравнения скорости Аррениуса (Arrhenius) и применили метод стохастического поля (Stochastic Field Method) для моделирования пространственной неоднородности влагосодержания. Путём разделения испарения воды и реакций пиролиза и введения локального члена теплоты испарения ошибка прогноза быстро снизилась с 18.5% до 3.2%. Этот урок научил нас, что при работе с биомассовыми материалами, такими как табак, влага — это отнюдь не просто переменная массовой доли, а динамический термодинамический регулятор.

Математическая обработка кинетики горения и подвижной границы

В этом разделе рассматриваются кинетика горения с чередованием пиролиза и окисления, а также математическая обработка задачи о подвижной границе и техники адаптивной сетки.

Ядро горения сигары заключается в чередовании пиролиза (Pyrolysis) и окисления (Oxidation). Полный цикл горения обычно следует такой логической цепочке: твёрдый табачный лист нагревается $\rightarrow$ вода испаряется $\rightarrow$ органическое вещество подвергается пиролизу с образованием летучих веществ $\rightarrow$ летучие вещества поступают в газовую фазу и вступают в реакции окисления $\rightarrow$ остаточный уголь (Char) окисляется на поверхности под контролем диффузии кислорода.

При вводе кинетических параметров мы должны применять многоступенчатую модель реакций. Одна энергия активации не может описать сложный химический состав табака. Обычно мы используем многокомпонентную модель пиролиза, разделяя табак на целлюлозоподобные, лигниноподобные компоненты и остаточную золу и присваивая каждому компоненту различные параметры Аррениуса ($A, E_a$).

Наиболее центральный алгоритм прогнозирования времени горения по сути является решением «задачи о подвижной границе» (Moving Boundary Problem), математически сходной с задачей Стефана (Stefan Problem). Положение фронта горения $\mathbf{s}(t)$ перемещается с течением времени, и его скорость перемещения $\dot{r}$ зависит от локального массопереноса и энергетического баланса:

$$\rho_s (1-\varepsilon) \frac{d\mathbf{s}}{dt} = \dot{m}_{pyro} + \dot{m}_{ox}$$ где $\dot{m}_{pyro}$ — скорость пиролиза, а $\dot{m}_{ox}$ — скорость окисления угля. В реальном процессе решения CFD для захвата этого чрезвычайно тонкого и интенсивного реакционного слоя мы должны применять технику адаптивного измельчения сетки (Adaptive Mesh Refinement, AMR). Если разрешение сетки недостаточно для покрытия градиента фронта горения, численная диффузия искусственно «сглаживает» скорости реакций, что приводит к завышению прогнозируемого времени горения. При решении таких сильно нелинейных задач выбор временного шага также является искусством — он должен удовлетворять ограничению числа Куранта (Courant Number), а также учитывать временной масштаб тепловой диффузии, иначе вычисление легко расходится в момент вспышки горения. ## Интеграция и инженерная реализация алгоритма прогнозирования Чтобы в конечном счёте получить точное время горения, мы построили интегральную модель от микроскопических реакций до макроскопической временной истории. Сначала с помощью переходного CFD-моделирования мы получаем поле температуры $T(x,t)$, поле концентрации кислорода $C_{O_2}(x,t)$ и локальное поле влагосодержания $M(x,t)$ в каждом положении микроэлемента в разные моменты времени. Затем, на основе этих полевых переменных, используя локальную скорость горения $\dot{r}(x,t)$, мы выполняем пространственно-временное интегрирование по геометрической длине $L$ сигары. Формула для прогнозирования общего времени горения $t_b$ может быть упрощённо выражена как: $$t_b = \int_0^L \frac{1}{\bar{\dot{r}}(x)} dx$$

где $\bar{\dot{r}}(x)$ — средняя линейная скорость вдоль осевого направления. Однако в инженерной практике, поскольку процесс горения сигары обладает импульсной характеристикой «зажигание-затухание-повторное зажигание» (особенно при моделировании процесса затяжки), этот интегральный процесс должен накапливаться во временном измерении. Нам необходимо точно сшивать во временной последовательности прогресс горения в течение каждого цикла затяжки (Puff) и стадию тления (Smoldering) между двумя затяжками.

При решении прогноза «стадии тления» другим ключевым моментом является задание граничных условий. После прекращения затяжки сигара находится в состоянии медленного окисления в среде с низким содержанием кислорода. Если граничное условие окружающей среды просто задать как воздух при нормальной температуре, смоделированная температура тления будет слишком высокой, что приведёт к занижению прогнозируемого времени горения. Нам необходимо ввести динамическую границу, учитывающую коэффициент конвективного теплообмена с окружающей средой $h$, а также радиационное рассеивание тепла с поверхности сигары.

Резюме и перспективы

В этом разделе представлен путь инженерной реализации от переходных полевых переменных CFD до интегрального прогнозирования времени горения.

Благодаря сложному вводу параметров и построению алгоритмов, описанным выше, мы можем не только точно прогнозировать время горения сигары, но, что более важно, наблюдать с помощью моделирования физические детали процесса горения — например, явление «неравномерного горения», вызванное слишком низкой проницаемостью покровных листьев, или явление «остановки горения», вызванное градиентами влажности.

Современные технологии моделирования развиваются в сторону более высоких измерений. С распространением технологии микротомографии (Micro-CT) мы можем получать чрезвычайно точные данные о микроскопической структуре пор табачных листьев и напрямую импортировать их в модели CFD, совершая переход от «приближения сплошной среды» к «реальной геометрической структуре». В то же время сопряжение алгоритмов машинного обучения с традиционной CFD — использование известных экспериментальных данных для обучения нейронных сетей с целью прогнозирования сложных параметров скорости реакций — возможно, станет следующим технологическим прорывом в решении проблемы сложности получения кинетических параметров горения табака и их высокой неопределённости.