烟雾扩散问题
1 研究背景简介
1.1 研究背景与意义
香烟烟雾中含有数千种化学物质，其中包括一氧化碳、尼古丁、苯等多种对人体有害的成分。这些物质在室内封闭环境中容易积聚，对人体健康构成显著威胁，尤其是二手烟和三手烟的危害已引起广泛关注。报告显示，2018年中国不吸烟者的二手烟暴露率为68.1%，其中家庭和工作场所是二手烟暴露的主要场所，这表明室内烟雾污染已成为影响公众健康的重要问题。因此，我们希望通过研究烟雾在室内环境中的扩散规律，为烟雾的控制、空气污染防治或室内空气质量管理提供参考，并帮助优化通风设计或制定健康防护措施。
烟雾的扩散过程受到室内温湿度、气流模式、房间结构及家具布置等多种因素的影响，具有较高的复杂性和随机性。深入研究烟雾扩散的传输机制，可以帮助揭示其污染范围与浓度分布，为优化室内空气质量管理提供理论支持。但在估计和预测吸入暴露风险时，由于与生物伦理学和动物保护相关的潜在限制，涉及人类志愿者和其他用于毒理学研究的哺乳动物替代模型的体内研究存在局限性。出于这样的考虑，我们希望通过理论分析建模与数值模拟的方式，得到烟雾扩散的一般规律，从而为提出可能的防控手段提供理论依据。
通过实验和仿真研究，能够验证不同通风方案对烟雾扩散的控制效果，为改善居住环境提供技术支持；利用数值模拟，可以精确再现不同环境下烟雾的扩散过程，从而更直观地评估各类干预措施的有效性。这种方法不仅降低了研究成本，还能为制定更加科学、精确的室内空气污染防控策略提供强有力的支持。
1.2 研究目标
本研究旨在探讨室内封闭环境下烟雾的扩散特性，基于扩散方程等物理原理与数学物理方法对烟雾的扩散情况进行解析求解与数值模拟，重点分析烟雾的流动路径、扩散速率和浓度分布，以优化空气管理和烟雾控制策略。
在具体的研究过程中，我们主要采用特定尺寸的密封容器来模拟室内封闭环境以便于实验开展与现象观察；同时，我们的研究主要关注扩散过程中烟雾的扩散效果、扩散速度以及不同时刻的浓度分布这三个要素，并通过定量计算测量与定性观察现象相结合的方式，得到烟雾扩散浓度空间分布随时间变化的的一般数学物理规律，从而为设计室内烟雾浓度控制系统提供理论依据。
2 理论分析求解
2.1 研究对象
在实际研究过程中，为方便实验的开展，我们选择了粒径分布在0.1~10微米范围内的植物甘油（丙三醇）气雾作为研究对象，这也是电子烟烟雾的主要成分之一。之所以使用植物甘油气雾来代替香烟烟雾进行研究，主要出于以下考量：
安全性：植物甘油为食品级材料，无健康与环境污染风险
可控性：植物甘油气雾生成过程可控，可进行定量化分析
易获取性：价格低廉，易于获取
可重复性：物质纯度高，化学性质稳定
具体而言，植物甘油气雾具有以下几条显著的理化特性：
无色无臭，有甜味，极易溶于水
密度（25℃）：1.26 g/cm^3
分子量： 92.09 g/mol
主要粒径属于典型气溶胶粒径范围
挥发性极低，以小颗粒悬浮于空气中
在研究过程中进行这样的材料替换，具有如下的合理性：
植物甘油气雾与香烟烟雾的气溶胶形成机制相同（蒸发-冷凝），二者在扩散和湍流行为上具有共性
二者粒径分布相似，具备相似的扩散动力学特征
二者在扩散过程中的视觉效果接近，便于追踪扩散行为
在本研究的理论分析部分，重点关注植物甘油气雾的浓度（即单位体积中该种物质颗粒的质量）ρ(x, y, z, t)这一物理量在研究区域（长为x_0 = 40cm，宽为y_0 = 50cm，高为z_0 = 60cm的长方体规则区域）空间范围内的分布情况及其随时间的变化情况。
2.2 数学物理方程的导出
针对烟雾扩散这一常见现象，有现成的实验规律——菲克定律（又称扩散定律）可以直接使用，其基本公式如下：
$$
\overrightarrow{q} = - D\nabla\rho
$$
其中：
- 式中负号代表扩散转移的方向（浓度减小的方向）和浓度梯度（浓度增大的方向）相反；
- D：扩散系数，与物质的种类以及环境温度有关；
- ∇ρ：浓度梯度，描述浓度不均匀分布的程度；
- 扩散流强度，即单位时间通过单位（横截）面积的质量，描述扩散运动的强度：
$$
\overrightarrow{q} = \frac{dm_{g}}{dVdt}
$$
在本研究中，针对常温环境（25℃）下植物甘油气雾在空气中的扩散系数，可进一步通过Chapman-Enskog方程
$$
D_{AB} = \frac{0.00143 \times T^{1.75} \times \sqrt{\frac{1}{M_{A}} + \frac{1}{M_{B}}}}{P \times \sigma_{AB}^{2} \times \Omega_{D}}
$$
进行估算，代入空气的平均分子量28.97 g/mol和估计的分子直径0.5nm可计算得到此时的扩散系数约为：
$$
D = 0.0962{cm}^{2}/s
$$
由于我们的研究区域是一个长方体的三维规则区域，故在研究坐标系的选择上，以长方体区域的底端顶点为原点，长、宽、高延伸方向作为$x、y、z$正方向建立空间直角坐标系，原点在长方体上的对角顶点坐标为$(x_{0},\ y_{0},z_{0})$。针对三个不同方向，还可以进一步地给出菲克定律的分量形式：
$$
q_{x} = - D\frac{\partial\rho}{\partial x},\ \ q_{y} = - D\frac{\partial\rho}{\partial y},\ \ q_{z} = - D\frac{\partial\rho}{\partial z}
$$
接下来将基于菲克定律和质量守恒定律给出三维扩散方程，即本研究涉及的数学物理方程：
对于空间中任一点(x, y, z)附近的无穷小空间而言：该空间内浓度变化取决于穿过其表面的扩散流。
考虑单位时间内x方向的扩散流：
左表面流入流量
$$
\left. \ q_{x} \right|{x}dydz
$$
右表面流出流量
$$
\left. \ q{x} \right|{x + dx}dydz
$$
故在单位时间内x方向净流入流量为
$$
\left. \ q{x} \right|{x}dydz - \left. \ q{x} \right|_{x + dx}dydz
$$
由于取无穷小空间，dx足够小，于是有
$$
\frac{\partial q_{x}}{\partial x} = \frac{\left. \ q_{x} \right|{x + dx} - \left. \ q{x} \right|{x}}{dx} 即 \left. \ q{x} \right|{x + dx} - \left. \ q{x} \right|{x} = \frac{\partial q{x}}{\partial x}dx
$$
代入化简后可以得到单位时间内x方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial x^{2}}dxdydz
$$
同理可得：单位时间内y方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial y^{2}}dxdydz
$$
单位时间内z方向净流入流量为
$$
D\frac{\partial^{2}\rho}{\partial z^{2}}dxdydz
$$
在本研究中，出于解析求解可行性的考虑，在理论分析部分暂时仅考虑常压下有源无汇的理想场景（即仅在某固定坐标点(x_0, y_0, z_0)处有一恒定扩散源，其强度（单位时间内单位体积产生的气体质量）为F(x,y, z,t)），而在后续数值求解与仿真模拟时才会把负压源视作汇的产生原因从而进一步补充对流的相关条件（对流项）。
在此条件下，基于质量守恒定律，单位时间在该无穷小空间内增加的气体质量等于单位时间内净流入的气体质量与扩散源产生的气体质量，其中扩散源产生的气体质量F(x,y, z,t)仅在(x_0, y_0, z_0)处有恒定值F_0而在其余空间各处均为0。于是列出并化简得到如下三维扩散方程（相比通用控制方程暂时不考虑对流项）：
$$
\frac{\partial\rho}{\partial t} - D\left( \frac{\partial^{2}\rho}{\partial x^{2}} + \frac{\partial^{2}\rho}{\partial y^{2}} + \frac{\partial^{2}\rho}{\partial z^{2}} \right) = \rho_{t} - D\nabla^{2}\rho = 0, 0 < x < x_{0}, 0 < y < y_{0}, 0 < z < z_{0}
$$
其中：
浓度的时间增长率（关于时间t的偏导数），即单位时间该无穷小空间内增加的气体质量:
$$
\frac{\partial\rho}{\partial t}
$$扩散系数D仅与物质种类和温度有关，若研究某种特定气体在理想的均匀介质（温度在空间各处保持一致且不随时间变化）条件下的扩散，该系数值为常数；反之则其也应表示为关于空间坐标与时间的函数D(x, y, z, t)（在确定气体类型时仅与温度的分布函数T(x, y, z, t)有关）。
2.3 初始条件与边界条件
该问题的初始条件是显然的：研究区域空间内在初始状态（t = 0）没有该种烟雾气体分布，且在某固定点处有一恒定扩散源，始终以恒定的扩散强度（单位时间内单位体积产生的气体质量）产生气体并向区域空间内扩散。将这样的初始条件写成数学表达式的形式：
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
式中F_0表示在固定坐标端点
$$
\left( \frac{x_{0}}{2},\ \frac{y_{0}}{2},\ 0 \right)
$$
处的恒定扩散源的扩散强度。结合实验时采用的实验装置，参考其技术文档可知，其值约为5g/cm^3。
关于边界条件，针对我们所研究的扩散问题，通常采用Dirichlet边界条件：
$$
\frac{\partial\rho}{\partial n} = 0
$$
显然这是一个齐次边界条件，也是一个典型的第二类边界条件。为方便后续求解处理，针对我们所建立的空间直角坐标系，可以将其进一步写成分量形式：
$$
\left{ \begin{array}{r}
\left. \ \frac{\partial\rho}{\partial x} \right|{x = 0} = \left. \ \frac{\partial\rho}{\partial x} \right|{x = x_{0}} = 0 \
\left. \ \frac{\partial\rho}{\partial y} \right|{y = 0} = \left. \ \frac{\partial\rho}{\partial y} \right|{y = y_{0}} = 0 \
\left. \ \frac{\partial\rho}{\partial z} \right|{z = 0} = \left. \ \frac{\partial\rho}{\partial z} \right|{z = z_{0}} = 0
\end{array} \right.\
$$
2.4 数学物理方程求解
针对2.2节提出的数学物理方程，可以通过分离变数法，结合2.3节给出的边界条件与初始条件，得到本征值问题的解和本征值的取值，并将多个本征值问题的解依次相乘得到本征解，再对所有可能的本征值求和，即可得到级数解。更进一步的，通过2.3节给出的初始条件，还可以对于级数解中的常数系数进行进一步确定，从而得到该数学物理方程完整的级数解解析表达式。
2.4.1 变量分离
对于扩散的烟雾而言，虽然其浓度会在空间中的不同位置随时间而变化，但时间与空间这两者在浓度瞬态变化的过程中实际上可以看作是相互独立的，即浓度随时间变化的速率ρ_t与浓度在空间分布中的变化速率∇ρ无关。
基于这样的假设，可以对烟雾浓度函数ρ(x, y, z, t)做如下的变量分离：
$$
\rho(x, y, z, t) = v(x,y,z) \bullet T(t)
$$
空间部分v(x,y,z)：描述烟雾浓度在三维空间中的分布，可以认为它与时间无关，或者说在某个时刻，密度的空间分布是静态的；
时间部分T(t)：描述烟雾密度随时间的变化，通常假设它是一个仅依赖于时间的函数，反映了扩散过程中气体密度随时间的衰减或增长。
在该式中，自变数x, y, z只出现于v之中，自变数t只出现于T之中，烟雾浓度的一般表示式具有分离变数的形式。
更进一步的，由于在空间中x、y、 z方向上的烟雾浓度变化情况也彼此独立，还可以对于空间部分v(x,y,z)做进一步的变量分离，即：
$$
v(x,y,z) = X\left( x) \bullet Y(y) \bullet Z(z \right)
$$
其中，X(x)为烟雾浓度在空间中x方向的分布情况，Y(y)为烟雾浓度在空间中y方向的分布情况，Z(z)为烟雾浓度在空间中z方向的分布情况。
2.4.2 常微分方程与本征值问题
显然三维扩散方程
$$
\rho_{t} - D\nabla^{2}\rho = 0
$$
为齐次方程，对于该齐次方程，可以将分离变量后的烟雾浓度代入三维扩散方程中，化简得到如下式子：
$$
T_{t}(t)v(x,y,z) = DT(t)\nabla^{2}v(x,y,z)，即\frac{T_{t}(t)}{DT(t)} = \frac{\nabla^{2}v(x,y,z)}{v(x,y,z)} = \frac{X^{‘’}(x)}{X(x)} = \frac{Y^{‘’}(y)}{Y(y)} = \frac{Z^{‘’}(z)}{Z(z)}
$$
其中，左边均为关于时间t的函数，与坐标位置x, y, z无关；右边则是关于坐标位置x, y, z的函数，与时间t无关。两边相等是显然不可能的，除非两边实际上是同一个常数，不妨将其记为-λ，于是有：
$$
\frac{T_{t}}{DT} = \frac{X^{‘’}}{X} = \frac{Y^{‘’}}{Y} = \frac{Z^{‘’}}{Z} = - \lambda
$$
将该式关于空间部分X(x)、Y(y)、Z(z)与时间部分T(t)分别分离，可以得到关于X、Y、Z以及关于T的一系列常微分方程：
$$
\left{ \begin{array}{r}
T^{‘} + \lambda DT = 0 \
\nabla^{2}v + \lambda v = 0
\end{array} \right.\ 即 \left{ \begin{array}{r}
T^{‘} + \lambda DT = 0 \
X^{‘’} + \lambda X = 0 \
Y^{‘’} + \lambda Y = 0 \
Z^{‘’} + \lambda Z = 0
\end{array} \right.\
$$
考虑2.3节给出的齐次边界条件即Dirichlet边界条件：边界处
$$
\frac{\partial\rho}{\partial n} = 0
$$
有：
$$
\left{ \begin{array}{r}
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0 \
Y^{‘}(0) = Y^{‘}\left( y_{0} \right) = 0 \
Z^{‘}(0) = Z^{‘}\left( z_{0} \right) = 0
\end{array} \right.\
$$
以上常微分方程与齐次边界条件共同构成了本征值问题。
2.4.3 本征值问题求解
首先，对于第一个方程
$$
\mathbf{T}^{‘} + \mathbf{\lambda DT} = \mathbf{0}
$$
而言，该方程为一阶常微分方程，可直接进行求解，结合初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
可得到：
$$
\mathbf{T}\left( \mathbf{t} \right) = \mathbf{F}_{\mathbf{0}}(\mathbf{1} - \mathbb{e}^{- \mathbf{\lambda Dt}})（C为某常数）
$$
其物理意义上的合理性如下：
初始时刻t = 0：
$$
T(0) = F_{0}\left( 1 - \mathbb{e}^{- 0} \right) = 0
$$
即空间内在初始状态下（t = 0）没有该种烟雾气体分布，与初始条件吻合；稳定时刻t->∞：
$$
T(\infty) = F_{0}\left( 1 - \mathbb{e}^{- \infty} \right) = F_{0}
$$
即浓度达到稳定值F_0，与扩散物理过程实际现象一致；时间项指数因子-λDt决定了扩散速度，其中D为扩散系数，越大扩散越快，λ为空间模态本征值，模态越高衰减越快，这也符合实际物理扩散过程。
而针对后面三个方程，与对应的齐次边界条件形成了三组本征值问题，由于其形式完全一致，下面将以其中一组本征值问题
$$
\left{ \begin{array}{r}
X^{‘’} + \lambda X = 0 \
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0
\end{array} \right.\
$$
为例，利用幂级数法求得其本征解：
选定点x_1 = 0，显然常微分方程
$$
X^{‘’} + \lambda X = 0
$$
的系数函数在该点领域内是解析的，因此该点为方程
$$
X^{‘’} + \lambda X = 0
$$
的常点。针对常点x_1，可以给出X(x)在该点邻域上的泰勒级数形式：
$$
X(x) = \sum_{k = 0}^{\infty}{a_{k}{(x - x_{1})}^{k}}
$$
将其代入常微分方程可得：
$$
\sum_{k = 0}^{\infty}{\lbrack(k + 2)(k + 1)a_{k + 2} + \lambda}a_{k}\rbrack{(x - x_{1})}^{k} = 0
$$
于是有：
$$
(k + 2)(k + 1)a_{k + 2} + {\lambda a}{k} = 0,\ k = 0,1,2\ldots
$$
可进一步得到系数递推公式：
$$
a{k + 2} = \frac{- \lambda}{(k + 2)(k + 1)}a_{k}
$$
推得：
$$
\left{ \begin{array}{r}
a_{2k} = \frac{\lambda^{k}{( - 1)}^{k}}{(2k)!}a_{0} \
a_{2k + 1} = \frac{\lambda^{k}{( - 1)}^{k}}{(2k + 1)!}a_{1}
\end{array} \right.\ ,\ k = 1,2,3\ldots
$$
因此可得到常微分方程
$$
X^{‘’} + \lambda X = 0
$$
的解：
$$
X(x) = a_{0}X_{0}(x) + a_{1}X_{1}(x)
$$
其中：
$$
X_{0}(x) = 1 + \frac{- \lambda}{2!}x^{2} + \frac{\lambda^{2}}{4!}x^{4} + \ldots + \frac{\lambda^{k}( - 1)^{k}}{(2k)!}x^{2k} + \ldots = \sum_{k = 0}^{\infty}{\frac{\lambda^{k}( - 1)^{k}}{(2k)!}x^{2k}} = \cos(\sqrt{\lambda}x)
$$
$$
X_{1}(x) = x + \frac{- \lambda}{3!}x^{3} + \frac{\lambda^{2}}{5!}x^{5} + \ldots + \frac{\lambda^{k}( - 1)^{k}}{(2k + 1)!}x^{2k + 1} + \ldots = \sum_{k = 0}^{\infty}{\frac{\lambda^{k}( - 1)^{k}}{(2k + 1)!}x^{2k + 1}} = \frac{\sin(\sqrt{\lambda}x)}{\sqrt{\lambda}}
$$
此时的收敛半径为
$$
R = \lim_{n \rightarrow \infty}\left| \frac{a_{n}}{a_{n + 2}} \right| = \lim_{n \rightarrow \infty}\left| \frac{(n + 2)(n + 1)}{- \lambda} \right| = \infty
$$
即级数解X_0(x)$和X_1(x)在复数域上始终收敛。X_0(x)仅含x的偶次幂，为偶函数；X_1(x)仅含x的奇次幂，为奇函数。同时，式中含有sqrt（λ）项，这要求λ≥0。
于是有：
$$
X^{‘}(x) = - a_{0}\sqrt{\lambda}\sin(\sqrt{\lambda}x) + a_{1}\cos(\sqrt{\lambda}x)
$$
将齐次边界条件
$$
X^{‘}(0) = X^{‘}\left( x_{0} \right) = 0
$$
代入可得:
$$
\left{ \begin{array}{r}
a_{1} = 0 \a_{0}\sqrt{\lambda}\sin(\sqrt{\lambda}x_{0}) + a_{1}\cos(\sqrt{\lambda}x_{0}) = 0
\end{array} \right.\
$$
显然有a_1= 0，但对于a_0而言，如果a_0也为0，这样的解是没有意义的；而要使得a_0不为零，就要求
$$
\lambda = 0或\sin\left( \sqrt{\lambda}x_{0} \right) = 0
$$
恒成立，这意味着
$$
\sqrt{\lambda}x_{0} = k\pi
$$
由于
$$
x_{0} > 0、\sqrt{\lambda} \geq 0
$$
于是k的可取值为
$$
k = 0,1,2,\ldots
$$
在这样的条件下，由于x_0与π均为定值，本征值λ则有常数k唯一确定，此时不妨将常数k视为新的本征值。于是可得到该本征值问题的本征解：
$$
X(x) = a_{0}\cos\left( \sqrt{\lambda}x \right) = a_{0x}\cos\left( \frac{k\pi}{x_{0}}x \right)，k = 0,1,2,\ldots
$$
与此类似的，另外两个本征值问题的本征解为：
$$
Y(y) = a_{0y}\cos\left( \frac{m\pi}{y_{0}}y \right)，m = 0,1,2,\ldots
$$
$$
Z(z) = a_{0z}\cos\left( \frac{n\pi}{z_{0}}z \right)，n = 0,1,2,\ldots
$$
值得注意的是，三个本征值问题中出现了三个新的本征值k、m、n，他们分别与时间部分涉及到的本征值λ有定量关系，但是当他们取不同值时，原本征值λ也会有不同的取值，即本质上有三个不同的本征值k、m、n取代了原先的本征值λ。
基于以上本征值问题的解和本征值的取值，可以将其各部分相乘得到本征解，再对本征解按照所有可能的本征值求和，即可得到烟雾浓度的级数解：
$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{a_{0x}a_{0y}a_{0z}F_{0}\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
由于式中a_0x、a_0y、a_0z、F_0均为常数（可由初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
进一步确定），不妨将其替换为某一常数
$$
A = a_{0x}a_{0y}a_{0z}F_{0}
$$
以便于后续计算，即此时可得到级数解为：
$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
2.4.4 级数解常数系数的确定
需要指出，函数
$$
\varnothing_{kmn} = \cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)
$$
是定义在长方体区域
$$
x \in \left\lbrack 0,x_{0} \right\rbrack,\ \ y \in \left\lbrack 0,y_{0} \right\rbrack,z \in \lbrack 0,z_{0}\rbrack
$$
上的一组正交基函数，其满足如下正交性与归一化条件：
$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\varnothing_{kmn}\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}} = \left{ \begin{array}{r}
x_{0}y_{0}z_{0},\ (k,m,n) = (k^{‘},m^{‘},n^{‘}) \
0,(k,m,n) \neq (k^{‘},m^{‘},n^{‘})
\end{array} \right.\
$$
对于级数解
$$
\rho(x,\ y,\ z,\ t) = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}\mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} = \sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\varnothing_{kmn}}}}{(1 - \mathbb{e}}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
而言，将两侧同时与Φ_k’m’n’相乘，并在区域内积分：
$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}} = \int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{(\sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{A\varnothing_{kmn}}}}\mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})\varnothing_{k^{‘}m^{‘}n^{‘}}dxdydz}}}
$$
利用正交性条件，上述式子右边仅保留(k,m,n) = (k’, m’, n’)的项，其余项积分为0：
$$
\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{kmn}dxdydz}}} = A(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})x_{0}y_{0}z_{0}
$$
代入初始条件
$$
\rho(x,y,\ z,0) = F_{0}\delta(x - \frac{x_{0}}{2})\delta(y - \ \frac{y_{0}}{2})\delta(z - 0)
$$
即可在t = 0处解出系数
$$
A = \frac{\int_{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x,\ y,\ z,\ t)\varnothing_{kmn}dxdydz}}}}{\left( 1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} \right)x_{0}y_{0}z_{0}}
$$
由于分子分母在t = 0处均为0，故利用洛必达法则上下对时间$t$求导得：
$$
\rho_{t} = T^{‘} = \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)Dt}（忽略系数）
$$
$$
A = \frac{\left. \rho_{t} \right|{t = 0}\int{0}^{z_{0}}{\int_{0}^{y_{0}}{\int_{0}^{x_{0}}{\rho(x, y, z, 0)\varnothing_{kmn}dxdydz}}}}{\left( {\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}\mathbb{e}}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t} \right)x_{0}y_{0}z_{0}}
= \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\varnothing_{kmn}\left( \frac{x_{0}}{2},\frac{y_{0}}{2},0 \right)
= \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\cos\left( \frac{k\pi}{x_{0}}\frac{x_{0}}{2} \right)\cos\left( \frac{m\pi}{y_{0}}\frac{y_{0}}{2} \right)\cos\left( \frac{n\pi}{z_{0}}0 \right) = \frac{F_{0}}{D\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}x_{0}y_{0}z_{0}}\cos\left( \frac{k\pi}{2} \right)\cos\left( \frac{m\pi}{2} \right)，其值随k,m,n的变化而变化
$$
最后，将其代入2.4.3节解得的级数解，可以最终得到植物甘油气雾浓度的解析表达式为：
$$
\rho(x,\ y,\ z,\ t) = \frac{F_{0}}{Dx_{0}y_{0}z_{0}}\sum_{k = 0}^{\infty}{\sum_{m = 0}^{\infty}{\sum_{n = 0}^{\infty}{\frac{1}{\left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)\pi^{2}}\cos\left( \frac{k\pi}{2} \right)\cos\left( \frac{m\pi}{2} \right)\cos\left( \frac{k\pi}{x_{0}}x \right)\cos\left( \frac{m\pi}{y_{0}}y \right)\cos\left( \frac{n\pi}{z_{0}}z \right)}}}(1 - \mathbb{e}^{- \left( \frac{k^{2}}{x_{0}^{2}} + \frac{m^{2}}{y_{0}^{2}} + \frac{n^{2}}{z_{0}^{2}} \right)D\pi^{2}t})
$$
2.5 基于MATLAB的级数解可视化
基于上述求解过程得到的烟雾浓度解析表达式，编写了对应的MATLAB程序，将不同时刻烟雾浓度在空间内各点处的分布情况进行计算与可视化。值得注意的是，由于级数解本质是无限项求和，而在MATLAB程序实际计算中需要指定求和项数，在此综合考虑计算精度与算力需求，选取$k,m,n$三个本征值各10项进行计算与可视化；除此之外，为与实际实验相对应，计算的扩散时间设置为30秒，同时为直观观察扩散过程中烟雾浓度随时间的变化，设置计算步长为1秒。
代入相应的数值并运行程序可以得到如下的可视化结果（代码详见附件1）：
可以看到，随着时间的推移，烟雾的可见范围越来越大（设定可见浓度为0.3mg/cm^3），且浓度从扩散中心点向四周呈递减趋势，这与我们的生活经验高度一致；在扩散时间为0s（即扩散未开始）时，空间区域内烟雾浓度均为零，这也与我们设定的初始条件相一致。以上事实初步验证了理论推导结果的正确性。
除此之外，为与后续的实验过程相对应，还记录了特定测量点(x_0/2, 0, z_0)处的烟雾浓度随时间的变化情况，并根据换算公式
$$
ppm = \frac{C \times 10^{6}}{\rho}
$$
将原有的气雾质量浓度C（单位：g/cm^3）换算为体积比ppm（百万分之一），得到的结果如下：
可以看到，由于该处离扩散中心点较远，此处烟雾浓度大概在扩散时间20s后达到基本稳定，稳定时的烟雾浓度值约为380ppm。
3 数值仿真模拟
3.1 数值分析方法：有限元法
由于严格的解析求解对于方程的复杂度以及初始条件与边界条件等都具有一定的要求，第2章所涉及到的理论分析也仅仅是围绕简化后的烟雾扩散问题展开。相比于真实的场景，这样的理论分析忽略了很多现实环境中的复杂变量（比如对流）。为使得我们的研究更具有实用价值，我们需要对于更加泛化的数学物理模型，采用数值求解的方法去模拟真实的扩散情况。
在这一部分，我们主要考虑到，2.3节提到的初始条件中烟雾源处视为单点并采用δ函数进行数学表示的做法过于理想化。与实验相对应的，我们采用的烟雾发生器的烟雾出口也不仅仅是一个单点，而是一个半径约为0.6cm的圆面。因此，我们希望通过有限元这一数值求解方法，对于初始条件迭代后的问题进行数值模拟分析。
图3.1 烟雾发生器烟雾出口孔径测量
基于有限元法的基本思想，我们主要按照如下的方法步骤进行数值分析：
首先进行区域剖分，将长方体研究区域Ω划分为有限个小单元，用于构建离散有限元模型。在划分单元时，采用长方体网格剖分方式，将区域划分为n_x * n_y * n_z个小立方体。每个小立方体单元有8个顶点，因此共有N = (n_x+1) * (n_y+1) * (n_z+1)个节点。对这些节点按列优先顺序进行编号，依次为1、2、…、N，同时对于各单元按体积分块顺序编号，编号为1、2、…、E。
接下来需要选择合适的单元基函数以表示单元内的浓度分布。在此我们选取基函数Φ_ i(x,y,z)近似表示浓度：
$$
\rho(x,\ y,\ z,\ t) \approx \sum_{i = 1}^{N}{c_{i}(t)\phi_{i}(x,y,z)}（c_{i}(t)表示节点i处的浓度值）
$$
并采用线性插值基函数（梯形单元），即单元内任意一点的浓度由8个顶点的浓度值线性插值（符号根据顶点位置的坐标取正或负）：
$$
\phi_{i}(x,y,z) = \frac{1}{8}(1 \pm x)(1 \pm y)(1 \pm z),\ i = 1,2,\ldots,8
$$
值得注意的是，在单元内，基函数在对应节点取值为1，其余节点取值为0，同时基函数也满足局部支撑性，即在单元外为零。
基于划分好的单元和选取的单元基函数，可以写出对应的单元积分表达式。通过对单元内的质量矩阵、刚度矩阵和源项向量进行积分，单元积分表达式应使得方程余量最小、本质边界条件余量最小、自然边界条件余量最小。将扩散方程乘以任意测试函数Φ_j(x,y,z)后，对区域Ω积分即可得到单元积分表达式：
$$
\int_{\Omega}^{}{\phi_{j}\frac{\partial\rho}{\partial t}\mathbb{d}\Omega} = D\int_{\Omega}^{}{\phi_{j}\nabla^{2}\rho\mathbb{d}\Omega} + \int_{\Omega}^{}{\phi_{j}F\mathbb{d}\Omega}
$$
该表达式中共有三项：
时间导数项（差分离散）：
$$
\int_{\Omega}^{}{\phi_{j}\frac{\partial\rho}{\partial t}\mathbb{d}\Omega} \approx \int_{\Omega}^{}{\phi_{j}\frac{\rho^{n + 1} - \rho^{n}}{\mathrm{\Delta}t}\mathbb{d}\Omega}
$$空间扩散项：通过分部积分，将扩散项转化为：
$$
\int_{\Omega}^{}{\phi_{j}\nabla^{2}\rho\mathbb{d}\Omega} = - \int_{\Omega}^{}{\nabla\phi_{j}\nabla\rho d\Omega} + \int_{\partial\Omega}^{}{\phi_{j}\nabla\rho \bullet \overrightarrow{n}dS}
$$
在齐次边界条件下，边界项为0。
- 离散形式：代入插值函数近似：
$$
\int_{\Omega}^{}{\phi_{j}\frac{\rho^{n + 1} - \rho^{n}}{\mathrm{\Delta}t}\mathbb{d}\Omega} = D\int_{\Omega}^{}{\nabla\phi_{j}\nabla(\sum_{i = 1}^{N}{C_{i}^{n + 1}\phi_{i}})\mathbb{d}\Omega} + \int_{\Omega}^{}{\phi_{j}F\mathbb{d}\Omega}
$$
基于列出的单元积分表达式，可以对各个单元进行分析，计算其对应的质量矩阵、刚度矩阵与源项向量，从而建立有限元方程。将离散化方程写为矩阵形式：
$$
M\frac{c^{n + 1} - c^{n}}{\Delta t} + Kc^{n + 1} = F
$$
其中：
质量矩阵：
$$
M_{\mathbb{i}j} = \int_{\Omega}^{}{\phi_{i}\phi_{j}d\Omega}
$$刚度矩阵：
$$
K_{\mathbb{i}j} = D\int_{\Omega}^{}{\nabla\phi_{i}\nabla\phi_{j}d\Omega}
$$源项向量：
$$
F_{\mathbb{i}} = \int_{\Omega}^{}{F\phi_{i}d\Omega}
$$
进而可以将单元积分表达式更新为有限元方程：
$$
(M + \Delta tk)c^{n + 1} = Mc^{n} + \Delta tF
$$
将每个单元映射到标准立方体单元（边长为2，中心为原点），可以根据基函数梯度
$$
\nabla\phi_{i} = (\frac{\partial\phi_{i}}{\partial x},\ \frac{\partial\phi_{i}}{\partial y},\ \frac{\partial\phi_{i}}{\partial z})
$$
在这些标准单元上使用高斯积分进行质量矩阵和刚度矩阵的计算：
$$
M_{\mathbb{i}j}^{(e)} = \int_{\Omega_{e}}^{}{\phi_{i}\phi_{j}d\Omega}
$$
$$
K_{\mathbb{i}j}^{(e)} = D\int_{\Omega_{e}}^{}{\nabla\phi_{i}\nabla\phi_{j}d\Omega}
$$
将源项向量
$$
F_{\mathbb{i}}^{(e)} = \int_{\Omega_{e}}^{}{F\phi_{i}d\Omega}
$$
离散化为节点的平均值后，也可以通过插值基函数对其进行积分。
完成单元分析后，可以将单元矩阵的局部节点编号1、2、…、8映射到总体矩阵的全局节点编号，并在每个单元中将其局部质量矩阵M^e、刚度矩阵K^e和源项向量F^e加入到总体矩阵M、K和向量F中，进行总体合成。
在对于总体有限元方程进行求解之前，还需要在总体矩阵方程中施加边界条件。在此我们采用消行修正法：对于齐次边界条件，将其带入总体有限元方程，并将边界上的节点浓度固定置零，同时修改总体矩阵对应行列为单位矩阵，右端项置零；对于非齐次边界条件（即初始条件），在边界节点区域（以扩散源中心（x_0/2, y_0/2, 0）为圆心，半径为0.6cm的圆面）中给定固定浓度值F_0，并对应调整右端向量F（已在单元积分中满足）：
$$
F_{i} \Leftarrow F_{i} - k_{ij}F_{0}
$$
最后，对于总体有限元线性方程
$$
(M + \Delta tk)c^{n + 1} = Mc^{n} + \Delta tF
$$
可以采用迭代方法进行求解：
初始条件：
$$
c^{0} = 0
$$时间步进：对于每个时间步，首先计算右端向量
$$
Mc^{n} + \Delta tF
$$
再采用稀疏矩阵求解器计算下一步的浓度
$$
c^{n + 1}
$$迭代终止：当时间达到设定终止时间T时，终止求解过程。
3.2 基于Python的有限元数值仿真
基于3.1节提出的有限元分析方法，利用Python中的MeshPy依赖库进行网格生成，并编写对应Python代码进行有限元数值仿真的计算与可视化。出于程序运行的稳定性等方面考虑，数值计算时仅选取原长方体区域的四分之一（高度不变，长宽各取一半，坐标原点即为扩散源中心）作为扩散区域。
运行程序（代码详见附录2），得到以下数值仿真结果：
可以看到，有限元仿真Python程序的运行结果能够较好地反映扩散过程中烟雾浓度分布的变化情况，基本与理论分析可视化中展现的趋势一致，在此不再赘述。相较于理论分析的可视化结果而言，有限元仿真的结果更加精细也更加丰富，特别是在截面的浓度分布上，很好地反映了扩散过程的基本特点，但在空间中的浓度分布可能受到了边界条件的过分影响，导致其明显呈现沿容器边缘扩散的趋势，而在容器内部的扩散范围较为有限。
3.3 基于COMSOL仿真软件的有限元数值仿真
尽管对于初始条件进行了一定的修正，但上述的有限元数值模拟仍然缺乏对于对流项的考虑，这可能是导致其结果受边界条件影响较大的原因之一。除此之外，上述的理论分析与数值计算也仅仅对于烟雾扩散过程中的浓度分布变化进行计算与模拟，而缺乏了烟雾扩散速度方面的计算分析，这也受到了解析求解方法的限制。因此，我们希望借助成熟的商业有限元仿真软件COMSOL进行更进一步的数值仿真，以分别模拟常压扩散与负压吸附两种情况下的植物甘油气雾颗粒扩散情况，并得到常压扩散情况下的扩散速度场与压力分布情况，帮助我们进一步完善烟雾扩散的基本模型。
为更加真实地反映实际扩散过程中的复杂物理环境，我们需要在COMSOL软件中指定相应的物理场模型：
湍流模型：帮助模拟湍流流动对气体扩散过程的影响，增加流体混合、热量传递和物质扩散等过程的复杂性；
RANS模型：通过对Navier-Stokes方程进行时间平均，简化了直接求解瞬时湍流方程的复杂度，计算成本相对较低，同时可以较为准确地描述气体扩散过程中湍流对流场的影响，尤其是在稳态流动情况下；
低雷诺数k-ω模型：采用湍动能（k）和湍流频率（ω）作为主要变量，能够较好地捕捉近壁区和低雷诺数流动特性，在复杂的流动区域中能更好地描述气体的湍流行为，提供更准确的边界层模拟，从而提高气体扩散的模拟精度。
除此之外，我们还需要使用温度场模型来模拟扩散过程中粒子的运动情况。对于烟雾产生端，采用650~750K的的随机梯度划分以达到不同速度的效果。
为了模拟出烟雾粒子扩散过程中的运动状态，我们还需要模拟出颗粒的流动。由于实际实验中采用的植物甘油气雾本质上是一种固体烟雾颗粒，所以在仿真中，我们选取遵循一阶牛顿规则的粒子，并给予每个颗粒一定的质量（与植物甘油气雾的理化性质相一致，设置其密度为1.26g/cm^3）。
处于仿真计算的复杂度与运行时间考虑，我们仅在x = x_0/2处的垂直截面区域内（50cm*60cm）进行仿真，并设置仿真时间为6秒、时间步长为0.1秒、网格单元大小为1cm。同时与实验环境保持一致，温度设置为25℃（298.15K）。
考虑到需要对于常压扩散与负压吸附两种情形分别进行仿真，还需要对于仿真边界（壁）进行一定的设置：
常压扩散：将各边设置为壁，壁条件设置为满散射；
负压吸附：将上边界设置为出口，边界条件设置为静压力（-0.02MPa）
完成了以上的所有设置流程，可以点击”计算”运行有限元仿真，仿真结果如下所示：
可以看到，通过COMSOL软件进行有限元数值仿真，可以有效划分三角网格并准确模拟出常压扩散情况下的扩散速度场及其强度分布情况，同时仿真得到的粒子运动轨迹大体趋势也与先前的解析求解与有限元数值计算较为吻合。仿真结果成功验证了先前理论模型推导的有效性，并在扩散速度方面进行了有力补充。
4 实验现象验证
4.1 预实验：常压封闭环境下的水雾扩散
在进行理论模型的推导与演算前，为更准确地把握初始条件、边界条件以及部分环境因素的设置，并对于烟雾扩散的具体现象进行初步了解，我们首先基于手头已有的实验材料进行了定性观察的预实验。在预实验中，我们使用香薰机作为烟雾发生装置，但由于装置限制，仅采用纯净水作为烟雾发生原料以产生水雾，通过观察水雾扩散现象来初步把握扩散过程的基本特征。
{width=”1.466816491688539in”
height=”1.9565212160979877in”}
基于以上预实验装置，搭建如下的实验装置场景：
启动香薰机，开始产生水雾，并通过长方体容器底部开孔向容器内部区域空间进行烟雾扩散。设置扩散时间为10分钟，拍摄视频记录扩散实验现象。
可以看到，随着扩散时间的推移，容器内聚集的水雾越来越多，且越靠近底部、越靠近扩散源中心，水雾的浓度也就越高（从现象来看，可以认为水雾团越密集的地方浓度越高），且当扩散时间达到一定程度时，容器内的水雾逐渐达到饱和，开始向容器底部沉降，并在容器底部呈聚集态势。
预实验的定性观察帮助我们对于烟雾扩散的基本现象有了宏观上的认知，并通过实物实验的方式让我们对于研究对象以及研究区域有了更加具象的认识，这对于我们后续理论模型推导与数值计算仿真的有序推进打下了良好的基础。但定性实验只能帮助我们认识现象，要对于理论模型的研究成果进行更加严谨的实证评估，还需要进行定量化的实验与测量。
4.2 正式实验：常压封闭环境下的植物甘油气雾扩散
在完成了理论模型的推导以及解析求解与数值仿真计算之后，我们还希望通过实际定量实验的测量结果来验证理论模型的合理性。为进行定量化实验，我们将预实验中的部分实验装置进行了替换，以严格定量控制烟雾的产生过程并实时测量特定测量点位处的烟雾浓度，从而与理论计算与仿真结果进行比照。
对于烟雾发生器而言，为达到定量产生植物甘油气雾的效果，选取了Selens手持烟雾发生器作为新的实验器材，该款烟雾发生器可以产生固定强度的植物甘油气雾，查阅其产品技术文档可知，其固定产生烟雾的强度约为F_0 ≈ 5g/cm^3（与先前理论推导采用相同符号）；同时为定量测量定点(x_0/2,0,z_0)处（与2.5节一致）的烟雾浓度，选取了MQ-2烟雾浓度传感器进行测量，通过STM32F1开发板对传感器进行供电，同时将数据传输至PC端并实时在屏幕上显示测量到的烟雾浓度值（单位：ppm，与g/cm^3的换算关系在2.5节已提及）。
需要注意的是，由于该款传感器读取烟雾浓度的过程是渐变而非突变，因此在使用前需要一定的预热（初始化）时间，且初始化完成后，在还未开始烟雾扩散实验时，读取到的烟雾浓度也并不为零，因此在后续数据处理时还需要移除对应的偏移量，实测平均偏移量在180ppm左右。
基于以上实验装置，搭建如下的实验装置场景：
对于该款手持烟雾发生器而言，需要手动长按开关才可持续发生烟雾，且该产品具有安全保护机制，最多只能连续产生30s的植物甘油气雾。启动装置后即开始进行烟雾扩散实验，烟雾通过长方体容器底部开孔向容器内部区域空间不断进行扩散。拍摄视频记录扩散实验现象与扩散过程中屏幕上显示的烟雾浓度的实时变化数值及曲线。
可以发现，随着扩散时间的增加，容器内聚集的植物甘油气雾也越来越多，且扩散的过程与理论解析求解可视化以及有限元数值仿真展示的烟雾扩散趋势基本一致。同时，在常压状态下，相较于水雾扩散而言，容器内聚集的植物甘油气雾在停止发生一段时间后仍然没有明显消散，这是由于植物甘油独特的理化性质所导致的，也比水雾扩散实验的结果更加接近香烟烟雾扩散的实际情况。
将烟雾传感器读取到的烟雾浓度数值实时传输到PC端保存，利用MATLAB进行可视化并与理论解析求解计算的定点测量结果进行比照，结果如下：
可以看到，通过烟雾传感器测量得到的烟雾浓度同样在20s左右达到基本稳定，且两者达到稳定后的浓度数值基本接近（在5%的误差允许范围内，约为380ppm），这说明我们的理论模型能够较好地解释实验测量结果；但比较而言，相比起实际测量的数据，理论求解得到的浓度在达到稳定前均高于实际测量结果，这一方面可能由于理论计算时进行了部分近似假设造成了偏差，另一方面也可能是实际测量时在传感器测量与数据读取传输过程中产生了一定的延迟。
通过理论求解计算结果与实验测量观察结果的比照可以发现，我们建立的理论模型对于烟雾扩散情况的模拟效果无论从现象（可视化效果）上还是从浓度数值上都能对于实际的烟雾扩散现象进行较好的模拟，这充分证明了我们理论解析求解模型研究的有效性。
4.3 正式实验：封闭环境中负压作用下的植物甘油气雾扩散
在正式实验阶段，除了在常压条件下进行烟雾扩散实验以验证烟雾扩散理论模型的准确性之外，为了后续进一步地提出可能的烟雾吸附方案并验证其有效性，我们还在负压环境下进行了对应的植物甘油气雾扩散实验，为后续开发提供参考。
在负压源的选取上，我们采用了涵道风机来营造负压环境。该款涵道风机可以通过旋钮调节产生不同档位强度的负压，通过查阅该产品技术文档可知，我们在实验时次啊用的最低档位负压值为P_0≈-0.02MPa，在通过COMSOL有限元仿真对于负压环境下的烟雾扩散进行数值模拟时也采用了相同的参数。其余的实验装置保持不变。
基于常压实验中搭建好的实验装置，在原烟雾发生处的开口对侧开出一个相同大小（半径3cm）的圆孔，并将涵道风机进风口对准该圆孔，完成负压吸附环境下植物甘油气雾扩散实验环境的搭建：
通过STM32F1开发板驱动涵道风机（供电）并通过旋钮将其调至最低档位，并采用与常压环境下实验相同的实验步骤与实验条件，拍摄视频记录扩散实验现象与扩散过程中屏幕上显示的烟雾浓度的实时变化数值及曲线。
可以发现，随着扩散时间的增加，容器内聚集的植物甘油气雾越来越多，但没有像常压状态下一样出现大量烟雾的堆积与沉降现象，且整体扩散方向受负压源牵引作用明显，这点与COMSOL有限元仿真得到的烟雾粒子运动轨迹基本一致；同时可以注意到，同样是在停止烟雾发生30秒后（期间持续进行负压吸附），在负压吸附的持续作用下容器内聚集的植物甘油气雾迅速地被排出容器外，而不像常压环境下的扩散那样会长时间堆积在容器内难以消散，这也充分验证了负压吸附这一解决方案的可行性与有效性，为后续的开发应用提供了良好的实验基础。
将烟雾传感器读取到的烟雾浓度数值实时传输到PC端保存，利用MATLAB进行可视化并与理论解析求解计算的定点测量结果进行比照，结果如下：
可以看到，相较于常压环境中的烟雾扩散，在负压作用下测量点处的植物甘油气雾浓度达到基本稳定的速度明显变快（大概在10s左右），且稳定后测得的烟雾浓度数值（100ppm左右）相较于常压下的稳定测量结果（380ppm左右）有明显下降，这样的数值结果在实验现象的基础上更进一步地验证了负压吸附这一解决方案的可行性与有效性，同时为后续的开发应用提供了良好的实验数据基础。
5 研究结论应用
5.1 研究结论
基于以上对于植物甘油气雾扩散问题在解析求解、数值模拟以及实验测量方面的深入研究，针对1.2节提出的研究目标，可以给出如下三点主要的研究结论：
在扩散过程中不同时刻的浓度分布方面，针对常压下密闭环境中的烟雾扩散现象，通过理论计算得到了不同时刻空间内各点处的浓度分布，并结合实验测试结果，印证了理论分析的高精度，在20s后空间内浓度分布基本稳定；
在扩散过程中烟雾的扩散速度方面，通过COMSOL有限元仿真软件对常压下密闭环境中的烟雾扩散进行数值模拟，得到扩散处于基本稳定状态下空间截面内速度场的分布，为分析浓度分布变化趋势提供参考；
在扩散过程中的烟雾扩散效果与除烟效果方面，将常压下和负压吸附下的COMSOL有限元仿真结果与实际实验现象进行对照，可以发现理论分析、数值求解以及实际实验三者的扩散现象结果基本一致，且在引入负压对流后，稳定后的烟雾浓度整体较常压时明显降低，且在烟雾发生停止后能够迅速清除残余的烟雾，具有良好的除烟效果。
除此之外，在这三点研究结论的基础上，我们还通过理论、仿真与实验的有机结合，得到了植物甘油气雾扩散的一般规律与数学物理模型，可以通过该理论模型对于实际的烟雾扩散情况进行定量化的数值模拟分析，从而对于烟雾扩散过程进行更精准的把控与必要的干预管理，并通过多种手段对该理论模型的有效性进行了充分验证，这也为我们后续应用产品的开发打下了良好的理论与实验基础。
5.2 研究成果应用：室内烟雾浓度控制系统
基于数值分析结果与研究结论，我们开发了一套简易的闭环控制系统作为原型样机，通过实时检测监测点处的烟雾浓度，实现对于风机吸附功率的灵活动态调控，以优化用户体验。
在该系统的开发方面，采用STM32F1开发板作为系统主控，通过ADC实时读取MQ-2传感器返回的烟雾浓度数值，并通过PID闭环控制实时调控给涵道风机供电的PWM占空比，从而实现对于风机吸附功率的实施调控。
通过运行该闭环控制系统并测试其烟雾浓度控制效果，可以发现该系统对于烟雾有着较高的灵敏度，并在检测到高浓度烟雾后能够迅速调动涵道电机作出反应，快速吸附容器内的植物甘油气雾并始终维持空间内的烟雾浓度低于150ppm，说明样机能够初步满足用户的基本除烟需求。
6 课程收获感悟
通过本次数学物理方法课程的学习，以及烟雾扩散课程项目的开展，我们通过在实际项目中运用课程中学习到的理论知识，在扩散过程中烟雾浓度分布的定量化求解过程中，对于数学物理方程的导出、边界与初始条件的确定以及基于变量分离与本征值方法求解级数解的过程有了更加深刻的认识与理解，同时通过课程中介绍的数值求解有限元方法对于更加泛化的情况进行了数值模拟与讨论。在解析与数值求解的理论分析基础上，我们还利用成熟的COMSOL有限元仿真软件对于烟雾扩散的粒子运动轨迹以及速度场分布进行数值仿真以作为对我们理论模型的补充，并通过定量化的实验测量进一步验证了理论模型的有效性，为我们后续进行实际的产品开发与应用提供了良好的理论与实验基础。
在本次课程项目报告的最后，再次感谢老师的耐心指导以及提出的宝贵意见，也希望我们能够将数学物理方法教给我们的定量化分析思维带入到未来的项目研究与学习工作中，对我们未来的项目开展与人生成长起到帮助。
参考文献
[1] Gao, Naiping and Jianlei Niu. “Modeling particle dispersion and deposition in indoor environments.” Atmospheric Environment (Oxford, England : 1994) 41 (2007): 3862 - 3876.
[2] Hoegg U. R. (1972). Cigarette smoke in closed spaces. Environmental health perspectives, 2, 117–128.
[3] Kuga K, Ito K, Yoo S-J, et al. First- and second-hand smoke dispersion analysis from e-cigarettes using a computer-simulated person with a respiratory tract model. Indoor and Built Environment. 2018;27(7):898-916.
[4] Holmberg, S., & Li, Y. (1998). Modelling of the Indoor Environment – Particle Dispersion and Deposition. Indoor Air, 8, 113-122.
[5] Al-sarraf, A.A., Yassin, M.F. & Bouhamra, W. Experimental and computational study of particulate matter of secondhand smoke in indoor environment. Int. J. Environ. Sci. Technol. 12, 73–86 (2015).
[6] 陈占秀,陈冠益,王艳,等.丙三醇与1,6-己二醇混合物降温凝固过程的分子动力学模拟[J].化工学报,2013,64(07):2316-2321.
[7] 沈玉峰,孟强龙.基于物理模型的烟雾扩散模拟[J].电子技术,2011,38(03):61-63.
[8] 李强,普小云.用毛细管成像法测量液相扩散系数——等折射率薄层测量方法[J].物理学报,2013,62(09):185-191.
[9] 杜青青.室内烟雾扩散模拟技术研究与实现[D].北京工业大学,2018.
[10] 梁昆淼编.数学物理方法[M].北京：高等教育出版社,1978:593页.
附录
1 MATLAB解析求解可视化代码
1 | % 参数定义 |
2 Python有限元仿真代码
1 | import numpy as np |
