多孔人造板VOC释放的时间分数阶传质模型

A Time-Fractional Mass Transfer Model for Predicting VOC Emissions from Porous Wood-Based Panel

摘要: 具有分形结构的多孔人造板散发的挥发性有机化合物(VOC)严重污染室内环境。本文首次提出时间分数阶传质模型分析人造板内VOC散发的反常扩散。并首次建立时间分数阶对流传质边界条件和质量平衡方程揭示VOC在人造板与空气界面处的散发特性。利用有限差分格式和改进的Nelder-Mead单纯形搜索与粒子群优化方法对模型进行数值优化。与Deng和Kim的模型作对比可知，本模型的数值模拟结果与实验数据吻合更好，且其相对误差Re (0.0034%)远小于先前模型(0.0257%)。这表明现有模型能更准确地揭示人造板中VOC释放的拖尾现象。此外，还分析了关键释放参数对VOC释放的影响，结果表明，较高的α、ε、D、N和较低的Kma均促进VOC的散发，这可为改善室内空气质量提供理论指导。

Abstract: Volatile organic compounds (VOC) emissions from porous wood-based panel with fractal structure seriously pollute indoor environment. In this paper, a time-fractional mass transfer model is proposed to analyze the anomalous diffusion of VOC emissions from wood-based panel. The time-fractional convective mass transfer boundary condition and mass balance equation are both developed for the first time. The finite difference scheme and the improved Nelder-Mead simplex search and particle swarm optimization are used to optimize the model numerically. Compared with Deng and Kim’s model, the numerical simulation results of present model are in better agreement with the experimental data, and the relative error Re (0.0034%) is much smaller than that of the previous model (0.0257%). This shows that the present model can reveal better the heavy-tailed phenomenon of VOC releases from wood-based panel more accurately. In addition, the influences of the key release parameters on VOC emission are analyzed. The results indicate that higher α, ε, D, N and lower Kma all can promote VOC emission, which can provide theoretical guidance for improving indoor air quality.

1. 引言

室内建筑材料人造板，如中密度纤维板(MDF)和刨花板(PB)，通常用于房屋内装饰、家具制造、地板、墙板等 [1]。人造板是由木材纤维中加入粘合剂再经热压而制成。VOCs主要从生产过程使用的粘合剂中释放，少量来源于木材本身 [2] [3]，易对人类健康产生严重的负面影响。在高浓度VOCs的室内环境中长期生活或工作会对人体造成不良建筑综合征，如呼吸系统疾病 [4] [5] [6]。此外，有研究结果表明癌症风险随着挥发性有机化合物浓度的升高而增加 [7] [8]。具有复杂结构的多孔人造板内VOC的释放过程复杂且周期较长。为了改善室内空气质量，阐明VOCs的释放机理和特性具有重要意义。

基于传质理论的数学模型已被广泛用于描述挥发性有机物在室内建筑材料中的扩散过程。例如，Little等人 [9] 提出了传质模型，以揭示VOCs从地毯中释放至环境舱扩散的行为。基于此，Zhang和Xu [10] 利用界面处的对流传质进一步改进了传质模型。随后，考虑到材料-空气界面处VOCs浓度的平衡关系，Deng和Kim [11] 不仅提出了一个具有充分解析优点的传质模型，而且还获得数值解。基于上述模型，Zhang和Xu [12] 讨论了建筑材料中初始浓度分布不均匀的影响。Wang等人 [13] 和Yang等人 [14] 研究了通风舱内建筑材料中VOCs的双面散发。He等人 [15] 提出了用于分析单层和复合木制品中VOCs释放特征的传质模型。虽然上述研究有效地模拟了VOCs的释放过程，但精度仍有待提高。

人造板是一种多孔建筑材料，其内部孔道排列无序，孔径大小不一。挥发性有机化合物的慢扩散与这种复杂的几何分形结构有关 [16]。Blondau等人 [17] 通过压汞法测试探索了多孔固体建筑材料中的孔径分布。Xiong等人 [18] 开发了一个双尺度模型，用于计算具有分形特征的室内多孔建筑材料的扩散系数。通过考虑干燥建筑材料的孔隙率影响，Deng等人 [19] 创建了一个多相传质模型，以揭示多孔材料VOCs释放。Cai等人 [20] 提出了一个解析模型，以刻画润湿液体进入多孔材料的自发吸入过程。Liu等人 [21] 根据分形理论详细揭示了中密度纤维板(MDF)和刨花板(PB)的分形结构。此外，Zhou等人 [22] 提出了一个用混合吸附剂修改多孔建筑材料中VOCs传质特性的预测系统。Zhou等人 [23] 提出的传质模型描述了多孔材料上漆膜从湿到干过程中的VOCs扩散。Su等人 [24] 从微观角度分析了不同类型多孔壁中VOC的释放过程。上述研究主要集中于人造板多孔分形结构对扩散系数和吸附的影响，而没有讨论复杂分形通道对VOCs扩散的时滞效应。

事实上，Fick定律不能准确描述自然界中许多反常扩散现象，比如超扩散和慢扩散。而分数阶微积分更适于描述反常扩散的空间非局部性和时间依赖性。例如，Yao等人 [25] 提出了一个时间分数平流扩散模型，以揭示腐植酸铁胶体负载砷在石英砂柱中的慢扩散行为。Chang等人 [26] 提出了一个时间分数阶对流扩散模型来描述多孔油气藏中气体扩散的后期拖尾现象。Kang等人 [27] 通过时空分数阶模型探索了煤基质中气体扩散的反常特征。基于分形理论，Zhokh等人 [28] 讨论了异丙醇在二氧化硅中的反常扩散行为。Zhang等人 [29] 提出了时空分数传质模型。Bai等人 [30] 将Fick定律扩展到多孔建筑材料中VOCs扩散的空间分数阶模型。到目前为止，出现了大量用于模拟反常扩散过程的分数阶导数模型 [31] [32] [33]。然而，人造板表面扩散通道对VOCs扩散的影响研究很少。

准确预测建筑材料释放的VOC是改善室内空气质量的基础。而基于Fick扩散的传质模型在拟合环境舱实验数据时存在较大误差，为了更准确地描述VOC的释放特性，本研究建立了时间分数阶传质模型，以探讨分形结构对VOC在人造板中反常扩散的影响。建立时间分数阶对流传质方程和质量平衡方程，利用有限差分和参数估计方法可以精确地模拟拖尾现象。此外，还以图形形式说明了关键释放参数对VOC释放浓度的影响。

2. 模型建立

将长度和厚度分别为l m和δ m的人造板作为测试材料放置在通风环境舱中进行试验。人造板的上表面暴露在空气中，在下表面不存在质量传递。由于其他四个侧面是用铝箔纸密封的，VOC在侧面不发生扩散。因而人造板中VOC散发的过程可看作一维瞬态传质。图1是测试材料在环境舱中释放VOC的物理示意图。x轴垂直于板的上表面，坐标原点位于板的底部。

Figure 1. The physical schematic diagram of VOC diffuses from wood-based panel in the environmental chamber

图1. 环境舱内人造板中VOC扩散的物理示意图

VOC分子与人造板孔壁间及分子与分子间的剧烈碰撞降低了分子的活性和自由选择扩散通道的机会，阻碍了VOC分子的扩散，进而导致多孔建筑材料中VOCs反常慢扩散。因此，为了描述反常扩散，将菲克定律推广为以下时间分数阶Fick定律 [28]：

$J=-{\lambda }^{1-\alpha }\frac{{\partial }^{1-\alpha }}{\partial t^{1-\alpha }}\left(D\frac{\partial C}{\partial x}\right)$, $0<\alpha <1$ (1)

其中J为VOC扩散质量通量，λ为延迟时间参数，其单位为秒(s)，α为时间分数阶导数参数，D为扩散系数，C为VOC浓度。这里定义Caputo分数阶导数算子 [34]：

${}_0{}^{C}D{}_t^{\alpha }C\left(x,t\right)=\frac{{\partial }^{\alpha }C\left(x,t\right)}{\partial t^{\alpha }}=\{\begin{array}{ll}\frac{1}{\Gamma \left(m-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\zeta \right)}^{m-\alpha -1}{C}^{\left(m\right)}\left(x,\zeta \right)\text{d}\zeta ,}\hfill & m-1<\alpha <m,\hfill \\ \frac{{\text{d}}^{m}C\left(x,t\right)}{\text{d}{t}^m},\alpha =m\in N^{{}_{*}}.\hfill & \hfill \end{array}$ (2)

这里m是一个正整数，Г(·)是伽玛函数。

描述VOC在人造板中扩散过程的传质方程为 [29]：

$\frac{\partial C}{\partial t}+\nabla \cdot J=0$ (3)

将式(1)代入式(3)得到一维瞬态传质方程：

$\frac{{\partial }^{\alpha }C}{\partial t^{\alpha }}={\lambda }^{1-\alpha }D\frac{{\partial }^{2}C}{\partial x^2}$ (4)

人造板内VOC浓度的初始条件为

$C\left(x,0\right)=C_0$, $0\le x\le \delta$ (5)

人造板材的下表面不存在传质通量，故边界条件为：

${\frac{\partial C}{\partial x}|}_{x=0}=0$ (6)

此外，VOC是通过人造板的表面孔隙释放的，因此面孔隙率ε对VOC释放的影响不可忽视。又研究表明多孔介质和许多合成多孔材料在多个尺度上表现出统计上的自相似分形尺度规律，多孔介质横截面积上的孔径分布服从分形尺度定律 [35]，那么人造板材表面的孔隙结构具有分形特征。根据以上分析，可以建立界面处的时间分数阶对流传质边界条件：

$-\epsilon D{\lambda }^{1-\alpha }{\frac{{\partial }^{1-\alpha }}{\partial t^{1-\alpha }}\left(\frac{\partial C}{\partial x}\right)|}_{x=\delta }=h\left(\frac{C\left(\delta ,t\right)}{K_{ma}}-C_a\right)$ (7)

其中h为对流传质系数，K_ma为板材上表面与空气界面处的分配系数，C_a为环境舱内空气中的VOC浓度。

为了更好地简化模型，保证VOC在实验舱内的质量守恒，假设：1) VOC在环境实验舱内混合均匀；2) 入口处通入干净空气；3) VOC不吸附在环境舱内壁。从而建立舱内质量平衡方程为：

$\frac{\text{d}{C}_a}{\text{d}t}=-L\epsilon D{\lambda }^{1-\alpha }{\frac{{\partial }^{1-\alpha }}{\partial t^{1-\alpha }}\left(\frac{\partial C}{\partial x}\right)|}_{x=\delta }-N{C}_a$ (8)

其中L为承载率，N为环境舱的换气率。

假设环境舱中VOC的初始浓度为

$C_a\left(0\right)=0$ (9)

由式(4)~(9)组成的系统封闭，其可描述人造板向环境舱释放VOC的全过程。

3. 数值方法

令 $x_i=i{h}_x$， $i=0,1,2,\cdots ,M$ ； $t_k=k\tau$， $k=0,1,2,\cdots ,Z$，其中空间步长为 $h_x=\delta /M$，时间步长定义为 $\tau =t/Z$。设 $C_i^k$ 是式(4)在网格点(x_i, t_k)处的数值解。

定理1 Caputo时间分数阶导数算子 $\alpha \left(0<\alpha <1\right)$ 可用有限差分L1算法离散为 [36]：

$\frac{{\partial }^{\alpha }f\left(t_k\right)}{\partial t^{\alpha }}=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[a_0^{\left(\alpha \right)}f\left(t_k\right)-{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left(a_{k-n-1}^{\left(\alpha \right)}-a_{k-n}^{\left(\alpha \right)}\right)f\left(t_n\right)-a_{k-1}^{\left(\alpha \right)}f\left(t_0\right)}\right]+{\rm O}\left({\tau }^{2-\alpha }\right)$ (10)

式(4)的时间分数阶导数项在格点(x_i, t_k)处的离散形式为：

$\frac{{\partial }^{\alpha }C\left(x_i,t_k\right)}{\partial t^{\alpha }}=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[a_0^{\left(\alpha \right)}C\left(x_i,t_k\right)-{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left(a_{k-n-1}^{\left(\alpha \right)}-a_{k-n}^{\left(\alpha \right)}\right)C\left(x_i,t_n\right)}-a_{k-1}^{\left(\alpha \right)}C\left(x_i,t_0\right)\right]+\text{O}\left({\tau }^{2-\alpha }\right)$ (11)

其中 $a_n^{\left(\alpha \right)}={\left(n+1\right)}^{1-\alpha }-n^{1-\alpha }$， $n=0,1,2,\cdots ,Z$， $a_0^{\left(\alpha \right)}=1$。

传质方程(4)、时间分数阶对流传质边界条件和质量平衡方程(7)~(8)的数值格式分别为

$\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[C_i^k-{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left(a_{k-n-1}^{\left(\alpha \right)}-a_{k-n}^{\left(\alpha \right)}\right)C_i^n-a_{k-1}^{\left(\alpha \right)}{C}_i^0}\right]={\lambda }^{1-\alpha }D\left[\frac{C_{i+1}^k-2C_i^k+C_{i-1}^k}{h_x^2}\right]$ (12)

$-\epsilon D{\lambda }^{1-\alpha }\frac{{\tau }^{\alpha -1}}{h_x\Gamma \left(1+\alpha \right)}\left[C_M^k-C_{M-1}^k-{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left(a_{k-n-1}^{\left(1-\alpha \right)}-a_{k-n}^{\left(1-\alpha \right)}\right)\left(C_M^n-C_{M-1}^n\right)-a_{k-1}^{\left(1-\alpha \right)}\left(C_M^0-C_{M-1}^0\right)}\right]=h\left(\frac{C_M^k}{K_{ma}}-C_a^k\right)$ (13)

$\begin{array}{c}\frac{C_a^k-C_a^{k-1}}{\tau }=-L\epsilon D{\lambda }^{1-\alpha }\frac{{\tau }^{\alpha -1}}{h_x\Gamma \left(1+\alpha \right)}\left[C_M^k-C_{M-1}^k-{\displaystyle \underset{n=1}{\overset{k-1}{\sum }}\left(a_{k-n-1}^{\left(1-\alpha \right)}-a_{k-n}^{\left(1-\alpha \right)}\right)\left(C_M^n-C_{M-1}^n\right)-a_{k-1}^{\left(1-\alpha \right)}\left(C_M^0-C_{M-1}^0\right)}\right]\\ -N{C}_a^k\end{array}$ (14)

其中 $C_a^k$ 为舱内VOC浓度的数值解。

此外，将模型的上述数值格式与一种改进的Nelder-Mead单纯形搜索和粒子群优化方法(MH-NMSS-PSO) [37] [38] 相结合可估计以下未知参数：扩散系数D、分配系数K_ma、初始浓度C₀、时间分数阶导数参数α、延迟时间参数λ和面孔隙率ε。估计参数的步骤详见文献 [29]。

4. 结果与讨论

本节利用MATLAB 2014b软件拟合实验数据，并与Deng和Kim的模型 [11] 进行比较，验证了本文模型的可靠性。揭示了关键释放参数对VOCs散发特性的影响。

4.1. 模型验证

为了进一步验证本模型的有效性，我们数值模拟了Yang等人 [39] 获得的刨花板PB释放至舱内的α-蒎烯浓度。表1给出了Deng和Kim的模型和本文模型的参数值。由图2可以发现，与Deng和Kim的模型 [11] 相比，本模型与实验数据具有更好的一致性。Deng和Kim的模型在释放前期、后期分别严重低估、高估了α-蒎烯的释放。此外，本文模型的相对误差Re为0.0034%，比Deng和Kim的模型的精度提升了86.7704%，这也表明本文模型是可靠的。因此，所提出的模型能够比先前模型更准确地捕捉后期拖尾现象。

Table 1. The parameter values of models

表1. 模型的参数值

Figure 2. The concentration of α-Pinene released from PB into the chamber

图2. 从刨花板中释放至环境舱内的α-蒎烯浓度

4.2. 灵敏度分析

本节分析了VOC散发特性对关键释放参数的敏感性，如图3~7所示。图3~7分别揭示了时间分数阶导数参数α、面孔隙率ε、扩散系数D、分配系数K_ma和换气率N对从人造板释放至舱内VOC浓度的影响。

Figure 3. The influence of the time fractional derivative parameter α on the VOC concentration released from wood-based panel into the chamber

图3. 时间分数阶导数参数α对从人造板释放至舱内VOC浓度的影响

图3中VOC的舱内浓度在50小时之前随α的增大而增大，即在释放前期α影响VOC释放的早期突破现象。一个有趣的现象是在约50小时之后，舱内VOC浓度随α的增大而减小，且较高α的相邻两条曲线间的交点时间比低α的要大，也即是说明人造板材的分形程度越低，后期拖尾现象越不显著。因而α主要影响人造板VOC释放后期的拖尾特性。

从图4可以发现，随着ε的增加，舱内VOC浓度和峰值浓度均升高。因此，ε主要影响人造板VOC释放的早期突破现象。

Figure 4. The influence of the areal porosity ε on the VOC concentration released from wood-based panel into the chamber

图4. 面孔隙率ε对从刨花板中释放至舱内VOC浓度的影响

Figure 5. The influence of the diffusion coefficient D on the VOC concentration released from wood-based panel into the chamber

图5. 扩散系数D对从刨花板中释放至舱内VOC浓度的影响

图5中在约40小时之前，峰值浓度随着扩散系数D的增大而增大，尔后，高D影响下的VOC浓度的变化速率逐渐增大，并依次与次高、次低、低D的影响曲线相交，交点的时间依次增大。此现象表明高D促进人造板材内VOC的释放，且D主要影响后期拖尾的反常扩散现象，即D越小，人造板释放VOC的后期拖尾现象越显著。

Figure 6. The influence of the partition coefficient K_ma on the VOC concentration released from wood-based panel into the chamber

图6. 分配系数K_ma对从刨花板中释放至舱内VOC浓度的影响

Figure 7. The influence of the air exchange rate N on the VOC concentration released from wood-based panel into the chamber

图7. 换气率N对从刨花板中释放至舱内VOC浓度的影响

从图6可知，VOC峰值浓度在30小时之前的释放前期随着分配系数K_ma的增大而减小，但其对早期突破现象比D的影响更大。在释放后期，高K_ma的VOC释放浓度分别与低、次低、次高的K_ma的影响浓度相交的时间点依次增大，即高K_ma的后期拖尾特性更加明显。因而，K_ma不仅影响人造板内VOC散发的早期突破，而且对后期拖尾现象也有影响。

图7为换气率N对人造板内VOC释放浓度的影响。大的换气率N促进VOC的散发。因而，需要保持较高的通风速率降低室内空气污染程度，进而提高室内空气质量。

5. 结论

本文首次建立了时间分数阶传质模型模拟人造板中VOC的反常扩散现象。首次尝试建立了时间分数阶的对流传质边界条件和质量平衡方程。数值模拟结果可知，与Deng和Kim的模型对比，本文所建立的模型与实验数据吻合的更好，表明本模型优于先前模型，且能够更好地刻画人造板释放VOC的反常扩散。此外，时间分数阶导数参数α和面孔隙率ε等关键释放参数对VOC释放特性的影响也以图形的形式显示出来。结果表明，VOCs的舱内浓度在释放前期随α、ε和D的增大而增大，在后期随α和D的增大而降低，然而K_ma与D呈现相反的规律。总之，较高的α、ε、D、N和较低的K_ma均促进VOC的散发，进而降低室内污染程度，且时间分数阶导数参数α可以用来刻画VOCs在多孔人造板中的后期拖尾的反常慢扩散行为。本研究可为改善室内空气质量提供理论依据。

基金项目

国家自然科学基金项目(No.12102032)，北京市自然科学基金–北京市教委联合项目基金(KZ201810016018)。