1 前言

干冰膨胀烟丝(DIET)技术是烟草工业中普遍使用的工艺技术。这项技术的研究始于20世纪70年代^[1-2]，其原理是用液态二氧化碳浸泡干烟丝，然后减压，形成干冰烟丝，再与高温工艺气体接触，干冰升华导致烟丝体积急速膨胀，从而增大烟丝的填充性能，由此可降低单箱卷烟的烟丝消耗量、减少焦油等危害性成分释放量、提升烟气透发性^[3-4]。

国内外关于DIET技术研究已有不少报道，主要涉及工艺改进、设备改造和膨胀前后化学组分变化等。Utsch等^[5]，de la Burde等^[6]在其专利中提出了以下工艺：用液态二氧化碳浸润烟丝，然后快速冷却使二氧化碳固化，再通过加热使二氧化碳升华，促使烟丝细胞膨胀，并对设备进行了详细介绍。在此基础上，Cho等^[7]提出了DIET生产的具体流程和完整设备的设计方案，使其适用于高堆积密度的烟丝处理过程。Enzell等^[8]对膨胀后烟丝组分的研究表明，膨胀后烟丝中焦油等有害有机物和卷烟的抽吸次数均有所降低。Halter等^[9]比较了二氧化碳、氨气和氟利昂等3种烟丝膨胀方法，结果表明，3种方法得到的烟丝膨胀率依次降低。Theophilus等^[10]对DIET进行毒理学研究，发现DIET与传统的三氯氟甲烷膨胀后的烟丝具有相似的生物活性。王志华等^[11]对二氧化碳浸渍缸进行了优化设计，使浸渍后的烟草较松散、不易结块，从而减少碎丝率，保证膨胀后烟丝的品质。姜建平^[12]对二氧化碳膨胀升华管的设计进行了改良，使整丝率提高了8.4%、碎丝率降低了1.54%。崇黎明等^[13]以烟丝的感官评价为评价指标，研究了运行参数对评价指标的影响。资文华等^[14]指出膨胀后烟丝含水率w₂(质量分数，下同)对膨胀烟丝感官质量有显著影响，当w₂为8.0% 时，感官评吸结果最好。韩金民^[15]通过Fluent软件对“Z”型膨胀塔进行了数值模拟，得到干冰烟丝在膨胀塔内运动的物流模型。然而，目前大部分研究以实验为主，缺少对DIET过程的相关数学模型研究，导致实际生产中难以对w₂进行精准控制，例如出现w₂逐步降低的异常现象。

本研究通过建立DIET过程的数学模型，着重分析工艺条件对w₂这一重要技术指标的影响，以期为改进生产工艺提供相关理论依据。

2 数学模型

如图 1所示，DIET过程膨胀段的主要流程为：饱和蒸气与回风(循环废气)混合后加热到300 ℃左右，形成工艺气体(其中水蒸气质量分数y_w约为76.4%，空气质量分数y_a约为23.6%)，工艺气体进入膨胀管；干冰烟丝(其中含水率w₁为18.9%，湿基干冰含量w_CO2为5.6%，均为质量分数，下同)与工艺气体在膨胀管中进行气流干燥；之后经旋风分离器进行气、固分离，出口工艺气体废气(质量流量为图 1中的L_o)部分循环作为回风，膨胀后的烟丝则进入下一工序，此时w₂降至8% 左右，烟丝体积膨胀率为70%~100%。为确保烟丝能够顺利进入膨胀管，管内需要维持一定的负压(约为−1kPa)，这会导致烟丝进口处吸进一部分空气(质量流量为图 1中的L_leak)。在实际操作中，每天需先对膨胀管系统进行一段时间的吹扫，然后再进料，1个操作周期约为5 h。一般情况下1 d的总操作时长大约10 h，即2个周期。

DIET过程本质上是一个过热蒸气气流干燥过程。烟丝在膨胀管内的停留时间很短^[16]，分别经历了干冰升华和烟丝干燥除去水分等2个阶段。第1阶段，烟丝中的干冰升华，烟丝体积迅速增大、温度恒定为−78.5 ℃(干冰升华温度)，干冰烟丝湿基干冰含量由5.6% 降至0；第2阶段，烟丝升温的同时进行干燥以除去水分，烟丝含水率由18.9% 降至8% 左右。由烟丝的干燥特性^[17]可知，烟丝干燥过程主要为恒速干燥。若将膨胀后的烟丝(w₂=8%)视为一个物质，且忽略设备热损失、膨胀过程中的烟丝比热容变化，2个阶段中烟丝吸收的热量Q_t可表示为：

$ Q_{\mathrm{t}}=M_{\mathrm{CO}_2} \gamma_{\mathrm{CO}_2}+W \gamma_{\mathrm{w}}+G_2 c_{\mathrm{m}}\left(t_{\mathrm{m} 2}-t_{\mathrm{m} 1}\right) $

(1)

式中：烟丝干冰升华量M_CO2=G₁w_CO2，烟丝干燥去除的总水量W=G₁w₁− G₂w₂。

再由烟丝进、出系统的总物料衡算可得：

$ G_1=G_2+W+M_{\mathrm{CO}_2}=G_2+\left(G_1 w_1-G_2 w_2\right)+G_1 w_{\mathrm{CO}_2} $

(2)

整理式(2)可得绝干烟丝质量流量G_c表达式：

$ G_{\mathrm{c}}=G_1\left(1-w_1-w_{\mathrm{CO}_2}\right)=G_2\left(1-w_2\right) $

(3)

烟丝由计量带送入膨胀塔，故G₁由计量带频率f和计量带上烟丝高度H控制，为烟丝密度ρ_p、烟丝流股截面宽度、H和带速等4者之积，其中前2者为定值，带速与f呈正比，因此，G₁可用常数ϕ(计量带流量系数)与H、f之积来表示：

$ {G_1} = \phi Hf $

(4)

热风干燥下工艺气体与烟丝间传热方式主要为对流传热，故Q_t亦可表示为：

$ Q_{\mathrm{t}}=\alpha S \Delta t $

(5)

式中：工程实际中，气体与烟丝间的平均传热温差∆t可采用算术平均，即

$ \Delta t=\frac{\left(t_{\mathrm{i}}-t_{\mathrm{m} 1}\right)+\left(t_{\mathrm{o}}-t_{\mathrm{m} 2}\right)}{2} $

(6)

联立式(1) ~(6)可得：

$ w_2=1-\frac{\phi H f \times\left(1-w_1-w_{\mathrm{CO}_2}\right)\left[\gamma_{\mathrm{w}}-c_{\mathrm{m}}\left(t_{\mathrm{m} 2}-t_{\mathrm{m} 1}\right)\right]}{\phi H f \times\left[w_{\mathrm{CO}_2} \gamma_{\mathrm{CO}_2}+\left(1-w_{\mathrm{CO}_2}\right) \gamma_{\mathrm{w}}\right]-\frac{\alpha S\left(t_{\mathrm{i}}+t_{\mathrm{o}}-t_{\mathrm{m} 1}-t_{\mathrm{m} 2}\right)}{2}} $

(7)

由式(7)可知，需先确定气体与烟丝间的对流传热系数α、传热面积S、膨胀后烟丝出口温度t_m2、工艺气体出口温度t_o等物理量数值，才能得到w₂的数学模型。

2.1 对流传热系数α及传热面积S 2.1.1 对流传热系数α

设每天开始进料的时刻为0时刻，此时的对流传热系数记作初始对流传热系数α₀，可用白丽青^[18]和Lee等^[19]总结的公式进行计算：

$ \alpha_0=\frac{N u \cdot \lambda_{\mathrm{g}}}{d_{\mathrm{e}}}=\frac{\lambda_{\mathrm{g}}}{d_{\mathrm{e}}}\left(2+0.74 {Re}_{\mathrm{r}}^{1 / 2} {Pr}^{1 / 3}\right)$

(8)

式中：$R e_{\mathrm{r}}=\frac{d_{\mathrm{e}}\left(u_{\mathrm{g}}-u_{\mathrm{p}}\right) \rho_{\mathrm{g}}}{\mu_{\mathrm{g}}} $；$ {Pr}=\frac{c_{\mathrm{g}} \mu_{\mathrm{g}}}{\lambda_{\mathrm{g}}} $。对于所涉及的工艺气体相关物性参数ρ_g、μ_g、λ_g、c_g，由于系统进、出口处的工艺气体浓度差别很小，均可按入口浓度估算，0时刻入口处工艺气体仅含有水蒸气和空气。具体估算式^[20-23]如下：

$ \rho_{\mathrm{g}}=y_{\mathrm{w}} \rho_{\mathrm{w}}+y_{\mathrm{a}} \rho_{\mathrm{a}} $

(9)

$\mu_{\mathrm{g}}=\frac{y_{\mathrm{w}} \mu_{\mathrm{w}} M_{\mathrm{w}}^{1 / 2}+y_{\mathrm{a}} \mu_{\mathrm{a}} M_{\mathrm{a}}^{1 / 2}}{y_{\mathrm{w}} M_{\mathrm{w}}^{1 / 2}+y_{\mathrm{a}} M_{\mathrm{a}}^{1 / 2}} $

(10)

$ \lambda_{\mathrm{g}}=\frac{y_{\mathrm{w}} \lambda_{\mathrm{w}} M_{\mathrm{w}}^{1 / 3}+y_{\mathrm{a}} \lambda_{\mathrm{a}} M_{\mathrm{a}}^{1 / 3}}{y_{\mathrm{w}} M_{\mathrm{w}}^{1 / 3}+y_{\mathrm{a}} M_{\mathrm{a}}^{1 / 3}} $

(11)

$ c_{\mathrm{g}}=y_{\mathrm{w}} c_{\mathrm{w}}+y_{\mathrm{a}} c_{\mathrm{a}} $

(12)

Re_r计算式中的烟丝颗粒的运动速度u_p可由烟丝单颗粒的自由沉降速度u_t(式(13))按式(14)~(16)修正得到。

对于气流干燥，其雷诺数一般为1~500，属过渡流状态，故u_t计算公式为^[24]：

$ u_{\mathrm{t}}=1.195 d_{\mathrm{p}}\left[\frac{\left(\rho_{\mathrm{p}}-\rho_{\mathrm{g}}\right)^2}{\rho_{\mathrm{g}} \mu_{\mathrm{g}}}\right]^{\frac{1}{3}} $

(13)

干冰烟丝膨胀过程属于烟丝颗粒群在工艺管道内的气力输送，直接使用式(13)估算u_p误差较大，因此需对其进行以下几方面修正。

(1) 竖直管内过渡流状态下，考虑到颗粒群对烟丝运动的影响，可将式(13)中的u_t修正为式(14)中的u'_p ^[24]：

$ u_{\mathrm{p}}^{\prime}=\frac{1}{1.5}\left[1-\left(\frac{u_{\mathrm{t}}}{u_{\mathrm{g}}}\right)^{1.5}\right] \cdot u_{\mathrm{g}} $

(14)

(2) 考虑到管壁对烟丝的运动速度也会产生影响，可将式(14)中的u_p′修正为式(15)中的$ u_{\mathrm{p}}^{\prime \prime} $ ^[25]：

$ u_{\mathrm{p}}^{\prime \prime}=\left[1-\left(\frac{d_{\mathrm{e}}}{D}\right)^{1.5}\right] \cdot u_{\mathrm{p}}^{\prime} $

(15)

(3) 其他修正：烟丝为丝状柔性细长颗粒，在运动中会发生“结绳现象”^[26]，另外，被二氧化碳浸渍形成的干冰烟丝会发生部分结块，以上2个原因均会导致烟丝平均粒径增大。此外，烟丝颗粒在经过弯管时，会发生聚集并贴壁流动，壁面会对颗粒产生摩擦阻力，同时烟丝在弯管处撞击壁面导致其运动轨迹不再是直线，十分复杂^[16]。以上各影响因素均会导致烟丝的实际运动速度降低，本研究将这些影响用1个小于1的修正系数θ表示，于是式(15)中的u_p″修正如下：

$u_{\mathrm{p}}=\theta u_{\mathrm{p}}^{\prime \prime} $

(16)

由于干冰烟丝在膨胀过程中会产生碎丝^[12]，且旋风分离器无法完全分离粒径较小的碎丝颗粒，部分碎丝颗粒会随废气进入回风，所以在1天两个操作周期内，设备中的工艺气体里的碎丝颗粒含量会逐渐变大，而碎丝颗粒含量的增大又会增大颗粒群间的撞击频率，引起α增大^[27]。所以本研究认为，在1天两个操作周期内，α是逐步增大的，可表示为：

$ \alpha = {\alpha _0} + \sigma \tau $

(17)

式(17)中，烟丝颗粒停留时间τ=0时，α =α₀，表示工艺气体中只有气相不含碎丝。联立式(8)~(17)可得：

$ \alpha=\frac{\lambda_{\mathrm{g}}}{d_{\mathrm{e}}}\left\{2+0.74\left[\frac{d_{\mathrm{e}}\left(u_{\mathrm{g}}-u_{\mathrm{p}}\right) \rho_{\mathrm{g}}}{\mu_{\mathrm{g}}}\right]^{1 / 2} \cdot\left(\frac{c_{\mathrm{g}} \mu_{\mathrm{g}}}{\lambda_{\mathrm{g}}}\right)^{1 / 3}\right\}+\sigma \tau $

(18)

式中：u_p由式(13)~(16)整理得到，$ u_{\mathrm{p}}=\frac{\theta}{1.5}\left[1-\left(\frac{d_{\mathrm{e}}}{D}\right)^{1.5}\right]\left[1-\frac{1.31 d_{\mathrm{p}}^{1.5}\left(\rho_{\mathrm{p}}-\rho_{\mathrm{g}}\right)}{\rho_{\mathrm{g}}^{0.5} \mu_{\mathrm{g}}^{0.5} u_{\mathrm{g}}^{1.5}}\right] \cdot u_{\mathrm{g}} $计算。

2.1.2 传热面积S

S为塔内烟丝的总外表面积，等于烟丝颗粒数N乘以单个烟丝颗粒的外表面积S_p，计算式为：

$ S=N \cdot S_{\mathrm{p}} $

(19)

式中：N可由烟丝的总质量(G₁τ)除以入口处单个烟丝颗粒质量m₀求取(见式(20))，S_p则按烟丝颗粒为片状来估算(见式(21))。

$ N=\frac{G_1}{m_0} \cdot \tau=\frac{G_1}{\frac{1}{6} \mathtt{π} d_{\mathrm{e}}^3 \rho_{\mathrm{p}}} \cdot \tau $

(20)

$S_{\mathrm{p}}=2\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right) $

(21)

烟丝颗粒在膨胀管中的流动可分为加速段和匀速段，由于加速段比匀速段要短得多，所以近似认为颗粒在膨胀塔内只做匀速运动，故式(18)、(20)中τ的表达式为：

$ \tau=l / u_{\mathrm{p}} $

(22)

将式(4)、式(20)~(22)代入式(19)可得：

$ S=\frac{12\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right) \cdot \phi H f}{\pi d_{\mathrm{e}}^3 \rho_{\mathrm{p}}} \cdot \frac{l}{u_{\mathrm{p}}} $

(23)

2.2 烟丝出口温度t_m2

在实际生产中，t_m2的监控十分困难，所以需要给出其计算式，以便估计其数值。以下计算假设：单个烟丝颗粒的几何尺寸很小，可近似认为颗粒内部温度分布均匀；烟丝中的干冰升华完全后才进行干燥；为简化计算，α采用不含碎丝时的α₀。

(1) 对于干冰升华段：

$ -\gamma_{\mathrm{CO}_2} \frac{\mathrm{d} m_{\mathrm{p}}}{\mathrm{d} \tau}=\alpha_0 S_{\mathrm{p}}\left(t_{\mathrm{g}}-t_{\mathrm{m} 1}\right) $

(24)

式(24)的初始条件为：τ=0时，$ m_{\mathrm{p}}=m_0=\frac{1}{6} \mathtt{π} d_{\mathrm{e}}^3 \rho_{\mathrm{p}} $。烟丝中的干冰升华完全时，干冰烟丝湿基干冰含量为0，m_p=m₀(1−w_CO2)，此时所需时间$ \tau_1=\frac{m_0 w_{\mathrm{CO}_2} \gamma_{\mathrm{CO}_2}}{\alpha_0 S_{\mathrm{p}}\left(t_{\mathrm{g}}-t_{\mathrm{m} 1}\right)} $。

(2) 对于烟丝干燥过程：

$\alpha_0 S_{\mathrm{p}}\left(t_{\mathrm{g}}-t\right)=-\gamma_{\mathrm{w}} \frac{\mathrm{d} m_{\mathrm{p}}}{\mathrm{d} \tau}+m_{\mathrm{p}} c_{\mathrm{m}} \frac{\mathrm{d} t}{\mathrm{d} \tau} $

(25)

由于烟丝干燥过程为恒速干燥，故dm_p/dτ=k为定值，代入式(25)得：

$ \alpha_0 S_{\mathrm{p}}\left(t_{\mathrm{g}}-t\right)=-k \gamma_{\mathrm{w}}+m_{\mathrm{p}} c_{\mathrm{m}} \frac{\mathrm{d} t}{\mathrm{~d} \tau} $

(26)

式(26)的定解条件为：τ=τ₁时，m_p=m₀(1−w_CO2)，t=t_m1。

从工程角度看，t_g可取为工艺气体的出、入口温度的算术平均值：

$ t_{\mathrm{g}}=\frac{\overline{t_{\mathrm{i}}}+\overline{t_{\mathrm{o}}}}{2} $

(27)

由以上微分方程组可得式(28)，通过式(28)可求取烟丝出口温度t_m2。

$ \left\{\begin{array}{l} t=t_{\mathrm{m} 1} \quad\left(\tau \leqslant \tau_1\right) \\ t=t_{\mathrm{g}}-\frac{k \gamma_{\mathrm{w}}}{2 \alpha_0\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right)}+a\left[\frac{1}{6 k} \mathtt{π} d_{\mathrm{e}}^3 \rho_{\mathrm{p}}\left(1-w_{\mathrm{CO}_2}\right)-\tau\right]^{\frac{2 \alpha_0\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right)}{k c_{\mathrm{m}}}} \quad\left(\tau>\tau_1\right) \end{array}\right. $

(28)

式中：$a=\left[\frac{k \gamma_{\mathrm{w}}}{2 \alpha_0\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right)}+t_{\mathrm{m} 1}-\frac{\bar{t}_{\mathrm{i}}+\bar{t}_{\mathrm{o}}}{2}\right] \cdot\left[\frac{1}{6 k} \mathtt{π} d_{\mathrm{e}}^3 \rho_{\mathrm{p}}\left(1-w_{\mathrm{CO}_2}\right)-\frac{m_0 w_{\mathrm{CO}_2} \gamma_{\mathrm{CO}_2}}{\alpha_0 S_{\mathrm{p}}\left(t_{\mathrm{g}}-t_{\mathrm{m} 1}\right)}\right]^{-\frac{2 \alpha_0\left(l_{\mathrm{p}} b+l_{\mathrm{p}} \delta+b \delta\right)}{k c_{\mathrm{m}}}} $。

2.3 回风温度

回风温度可近似视为与t_o相同。根据整个系统的热量衡算可得：

$ L_{\mathrm{i}} c_{\mathrm{g}} t_{\mathrm{i}}+L_{\text {leak }} c_{\mathrm{a}} t_{\infty}+\left(G_2 c_{\mathrm{m}}+M_{\mathrm{CO}_2} c_{\mathrm{d}}+W c_{\text {ice }}\right) t_{\mathrm{m} 1}=L_{\mathrm{o}} c_{\mathrm{g}} t_{\mathrm{o}}+G_2 c_{\mathrm{m}} t_{\mathrm{m} 2}+M \gamma_{\mathrm{CO}_2}+W \gamma_{\mathrm{w}}+Q_1 $

(29)

再由干燥过程总物料衡算式$ L_{\mathrm{i}}+L_{\text {leak }}+G_1=L_{\mathrm{o}}+G_2 $及式(2)可得：

$ L_{\mathrm{o}}=L_{\mathrm{i}}+M_{\mathrm{CO}_2}+W+L_{\text {leak }} $

(30)

式中：L_i=ρ_gu_gπD²/4，L_leak可由烟丝进口管道尺寸与壁内风速估算出来，正常工况下约为219 kg⋅h⁻¹或0.0608 kg⋅s⁻¹。为简化运算，式(29)的等号左侧的W和式(30)中的W按式(31)计算：

$ W = {G_1}{w_1} - {G_2}{w_2} $

(31)

由热量衡算可知，式(29)中的膨胀塔管壁与环境之间的热损失Q_l，与工艺气体与污垢层间的对流传热量、污垢层内的导热量、金属管壁内的导热量、保温层内的导热量、保温层外壁与环境之间的对流传热量均相等，于是得到总传热方程为：

$ Q_1=K_{\mathrm{o}} \cdot \mathtt{π} D_2 l \cdot\left(t_{\mathrm{g}}-t_{\infty}\right) $

(32)

式中：K_o的表达式为

$ K_{\mathrm{o}}=\frac{1}{\frac{1}{\alpha_{\text {in }}} \cdot \frac{D_2}{D}+R_{\text {si }} \cdot \frac{D_2}{D}+\frac{D_1-D}{\lambda_{\text {wall }}} \cdot \frac{D}{D_{\mathrm{m}}}+\frac{D_2-D_1}{\lambda_{\text {ins }}} \cdot \frac{D}{D_{\text {m }}}+\frac{1}{\alpha_{\text {out }}}} $

(33)

式中：管内工艺气体与污垢层间的对流传热系数α_in可用齐德-泰特公式^[28](式(34))计算、保温层外壁与环境之间的对流传热系数α_out可用经验公式(式(35))计算^[28]：

$ \alpha_{\mathrm{in}}=\frac{\lambda_{\mathrm{g}}}{D} N u=\frac{\lambda_{\mathrm{g}}}{D} \cdot\left[0.027 {Re}_{\mathrm{g}}^{0.8} {Pr}^{1 / 3}\left(\frac{\mu_{\mathrm{g}}}{\mu_{\mathrm{w}}}\right)^{0.14}\right]$

(34)

$ \alpha_{\text {out }}=9.8+0.052\left(t_{\text {out }}-t_{\infty}\right) $

(35)

式中：Re_g=Du_gρ_g/μ_g，(μ/μ_w)^0.14=0.95。

联立式(24)~(35)可得到t_o表达式：

$ t_{\mathrm{o}}=\frac{1}{B}\left[A t_{\mathrm{i}}+C-\frac{C+(A+B) t_{\mathrm{i}}-2 B t_{\infty}}{\frac{2 B}{\pi D_2 l}\left(\frac{1}{\alpha_{\text {in }}} \cdot \frac{D_2}{D}+R_{\mathrm{si}} \cdot \frac{D_2}{D}+\frac{D_1-D}{\lambda_{\mathrm{w}}} \cdot \frac{D}{D_{\mathrm{m}}}+\frac{D_2-D_1}{\lambda_{\text {ins }}} \cdot \frac{D}{D_{\mathrm{m}}}+\frac{1}{\alpha_{\text {out }}}\right)+1}\right]$

(36)

式中：$ A=\frac{\pi D^2}{4} \cdot \rho_{\mathrm{g}} u_{\mathrm{g}} c_{\mathrm{g}}-0.5 \alpha_0 S $，$B=\left[\frac{\pi D^2}{4} \cdot \rho_{\mathrm{g}} u_{\mathrm{g}}+L_{\text {leak }}+\phi H f \cdot\left(w_{\mathrm{CO}_2}+\frac{w_1+w_{\mathrm{CO}_2}-w_2}{1-w_2}\right)\right] c_{\mathrm{g}}+0.5 \alpha_0 S $，$ C=L_{\text {leak }} c_{\mathrm{a}} t_{\infty}+\phi H f \cdot\left(w_{\mathrm{CO}_2} c_{\mathrm{d}}+\frac{w_1+w_{\mathrm{CO}_2}-w_2}{1-w_2} c_{\mathrm{i}}\right) t_{\mathrm{m} 1}+0.5 \alpha_0 S\left(t_{\mathrm{m} 1}+t_{\mathrm{m} 2}\right) $。

将式(18)、(23)、(28)、(36)代入式(7)，可得到w₂的最终数学模型，其中，θ和σ作为拟合参数需通过实验数据拟合得到。

3 模型参数的标定 3.1 DIET实验

数学模型需要的设备结构参数、烟丝和工艺气体等物性参数以及膨胀烟丝实验的工艺条件列于表 1。

由表 1可见，工艺条件t_∞、t_out、L_leak、w₁、w_CO2、t_m1均为定值，而DIET实验中属于操作控制的工艺条件f、H、u_g、t_i在1个操作周期(5 h即300 min)内的变化情况如图 2、3所示。由图 2可见，f为28~30 s⁻¹、H为1.60~1.76 m；由图 3可见，u_g为35.9~36.2 m⋅s⁻¹、t_i为300~305 ℃。生产过程中，每10 s取1组数据点；处理数据时，每6个点即1 min取1个均值，共监控300 min即300个数据点。

3.2 实验结果及模型标定

通过实验得到的t_o和w₂的数值见图 4、5。

结合图 2、3及表 1中的参数与工艺条件，使用Matlab软件、调节拟合参数θ和σ，对实验数据进行拟合，求出过程目标函数F的最小值，即可得到最佳θ和σ。F表达式为：

$ F=\frac{1}{N_{\mathrm{p}}} \sum\limits_{i=1}^{N_{\mathrm{p}}}\left\{\frac{\left[\left(\hat{w}_2\right)_i-\left(w_2\right)_i\right]^2}{\left(w_2\right)_i^2}+\frac{\left[\left(\hat{t}_{\mathrm{o}}\right)_i-\left(t_{\mathrm{o}}\right)_i\right]^2}{\left(t_{\mathrm{o}}\right)_i^2}\right\}$

(37)

最终拟合结果为：当θ为0.992 2、σ为1.833×10⁻⁷ kW⋅(m²⋅℃⋅s)⁻¹时模型效果最好。

对实验中点(第150 min)的数据应用第2节得到的数学模型，计算得到$ t_{\mathrm{o}} $和$ {w_2} $的关联值及若干中间参数(见表 2)，表 2同时还列出了实验值。

使用平均绝对误差(MAPE)来评估通过数学模型计算出的关联值$ {\hat{t}_{\mathrm{o}}} $、$ {\hat{w}_2} $，MAPE定义式如下：

$ \mathrm{MAPE}_{t_{\mathrm{o}}}=\frac{100 \%}{N_{\mathrm{p}}} \sum\limits_{i=1}^{N_{\mathrm{p}}} \frac{\left|\left(\hat{t}_{\mathrm{o}}\right)_i-\left(t_{\mathrm{o}}\right)_i\right|}{\left(t_{\mathrm{o}}\right)_i} $

(38)

$ \mathrm{MAPE}_{w_2}=\frac{100 \%}{N_{\mathrm{p}}} \sum\limits_{i=1}^{N_{\mathrm{p}}} \frac{\left|\left(\hat{w}_2\right)_i-\left(w_2\right)_i\right|}{\left(w_2\right)_i} $

(39)

评估结果为：对于t_o，MAPE为0.24%，最大误差为0.83%；对于w₂，MAPE为1.99%，最大误差为6.33%。可见，t_o和w₂的拟合精度都较好，同时$ {\hat w_2} $也有逐步降低趋势，与实际情况相符。

4 结论

本研究针对卷烟厂制丝车间DIET工艺，基于物料衡算和热量衡算，推导出了τ、f、u_g、t_i等操作条件与t_m2、t_o、w₂之间的数学模型。t_o和w₂的MAPE分别为0.24% 和1.99%，表明模型关联值与实验值吻合良好，能较准确地估算干冰烟丝膨胀前、后的含水率变化值，对DIET工艺过程w₂的控制具有较好的指导意义。

符号说明：

A   — 式(36)中的指代量    t_o   — 工艺气体出口或回风温度，℃

a   — 式(28)中的指代量    $ \bar t_{\mathrm{o}} $   — 工艺气体出口或回风平均温度，℃

B   — 式(36)中的指代量    $ \hat t_{\mathrm{o}} $   — t_o的关联值，℃

b   — 烟丝颗粒宽度，m     t_out   — 保温层外壁温度，℃

C   — 式(36)中的指代量    t_s   — 饱和蒸气温度，℃

c_a   — 空气的比热容，kJ⋅kg⁻¹⋅℃⁻¹     t_∞   — 环境空气温度，℃

c_d   — 干冰的比热容，kJ⋅kg⁻¹⋅℃⁻¹     ∆t   — 气体与烟丝间的平均传热温差，℃

c_g   — 工艺气体的比热容，kJ⋅kg⁻¹⋅℃⁻¹     u_g   — 工艺气体风速，m⋅s⁻¹

c_ice   — 冰的比热容，kJ⋅kg⁻¹⋅℃⁻¹     u_p   — 烟丝颗粒运动速度，m⋅s⁻¹

c_m   — 不含干冰的烟丝比热容，kJ⋅kg⁻¹⋅℃⁻¹     u_p′   — 第一次修正后竖直管颗粒群的运动速度，m⋅s⁻¹

c_p   — 干冰烟丝比热容，kJ⋅kg⁻¹⋅℃⁻¹     u_p″   — 第二次修正后竖直管颗粒群的运动速度，m⋅s⁻¹

c_w   — 水蒸气比热容，kJ⋅kg⁻¹⋅℃⁻¹     u_t   — 烟丝颗粒自由沉降速度，m⋅s⁻¹

D   — 膨胀塔内径，m     W   — 烟丝干燥总水量，kg⋅s⁻¹

D₁   — 膨胀塔外径，m     w₁   — 入口处干冰烟丝含水率

D₂   — 保温层外径，m     $ {\bar w_1} $   — 入口处干冰烟丝平均含水率

D_m   — 膨胀塔平均直径，m     w₂   — 出口处膨胀后烟丝含水率

d_e   — 烟丝颗粒的当量直径，m     $ {\bar w_2} $   — 出口处膨胀后烟丝平均含水率

d_p   — 颗粒直径，m     $ {\hat w_2} $ ⎯ w₂的关联值

F   — 目标函数    w_CO2   — 入口处干冰烟丝湿基干冰含量，无因次

f   — 计量带频率，s⁻¹     y_a   — 工艺气体含空气率，无因次

G₁   — 入口处干冰烟丝质量流量，kg⋅s⁻¹     y_w   — 工艺气体含水率，无因次

G₂   — 出口处膨胀后烟丝质量流量，kg⋅s⁻¹     α   — 气体与烟丝间的对流传热系数，kW⋅m⁻²⋅℃⁻¹

G_c   — 绝干烟丝质量流量，kg⋅s⁻¹     α₀   — 气、固两相间的初始对流传热系数，kW⋅m⁻²⋅℃⁻¹

H   — 计量带上的烟丝高度，m     α_in   — 气体与内壁间的对流传热系数，kW⋅m⁻²⋅℃⁻¹

K_o   — 膨胀塔总传热系数，kW⋅m⁻²⋅℃⁻¹     α_out   — 保温层与环境间的对流传热系数，kW⋅m⁻²⋅℃⁻¹

k   — 烟丝的质量变化率，kg⋅s⁻¹     δ   — 烟丝颗粒厚度，m

L_i   — 进风量，kg⋅s⁻¹     ϕ   — 计量带流量系数，kg⋅m⁻¹

L_leak   — 设备漏风量，kg⋅s⁻¹     γ_CO2   — 干冰升华热，kJ⋅kg⁻¹

L_o   — 回风量，kg⋅s⁻¹     γ_w   — 水的相变热(包括融化与汽化)，kJ⋅kg⁻¹

L_s   — 入口饱和蒸气流量，kg⋅s⁻¹     λ_a   — 空气的导热系数，kW⋅m⁻¹⋅℃⁻¹

l   — 膨胀塔长度，m     λ_g   — 工艺气体的导热系数，kW⋅m⁻¹⋅℃⁻¹

l_p   — 烟丝颗粒长度，m     λ_ins   — 保温层的导热系数，kW⋅m⁻¹⋅℃⁻¹

MAPE   — 绝对百分比误差，无因次    λ_wall   — 金属管壁的导热系数，kW⋅m⁻¹⋅℃⁻¹

M_a   — 空气的摩尔质量，kg⋅mol⁻¹     λ_w   — 水蒸气的导热系数，kW⋅m⁻¹⋅℃⁻¹

M_CO2   — 烟丝干冰升华量，kg⋅s⁻¹     μ_a   — 空气的黏度，Pa⋅s

M_w   — 水的摩尔质量，kg⋅mol⁻¹     μ_g   — 工艺气体的黏度，Pa⋅s

m_p   — 单个烟丝颗粒质量，kg     μ_w   — 水蒸气的黏度，Pa⋅s

m₀   — 入口处单个干冰烟丝颗粒质量，kg     θ   — 烟丝颗粒运动速度修正系数，无因次

N   — 膨胀塔中烟丝的颗粒数，无因次    ρ_a   — 空气密度，kg⋅m⁻³

N_P   — 数据总量，无因次    ρ_g   — 工艺气体密度，kg⋅m⁻³

Nu   — 努塞尔数，无因次    ρ_p   — 烟丝密度，kg⋅m⁻³

n   — 每秒进入膨胀塔的烟丝颗粒数，s⁻¹     ρ_w   — 水蒸气密度，kg⋅m⁻³

Pr   — 普朗特数，无因次    σ   — 对流传热系数变化速率，kW⋅m⁻²⋅℃⁻¹⋅s⁻¹

p_s   — 入口饱和蒸气压力，kPa     τ   — 烟丝颗粒停留时间，s

Q_l   — 过程热损失，kW     τ₁   — 干冰升华时间，s

Q_t   — 烟丝吸收的热量，kW     下标

Re_g   — 工艺气体的雷诺数，无因次    a   — 空气

Re_r   — 烟丝颗粒的真实雷诺数，无因次    g   — 工艺气体

R_si   — 内壁污垢热阻，m²⋅℃⋅kW⁻¹     in   — 塔壁内侧

S   — 塔内烟丝的总外表面积即传热面积，m²     out   — 保温层外侧

S_p   — 单个烟丝颗粒的外表面积，m²     p   — 烟丝颗粒

t   — 烟丝温度，℃     s   — 饱和蒸气

t_g   — 工艺气体平均温度，℃     w   — 水蒸气

t_i   — 工艺气体入口温度，℃     d   — 干冰

$ {\bar t_i} $   — 工艺气体入口平均温度，℃     上标

t_m1   — 干冰烟丝进口温度，℃     ^   — 关联值

t_m2   — 膨胀后烟丝出口温度，℃     —   — 平均值