沈阳建筑大学 冯国会 毕扬 蔡易霖 曹广宇
重庆大学 张亿先
摘 要：病毒性微生物气溶胶污染着室内的环境，危害着人们的身体健康。运用不同通风手段排除室内污染物是近些年众多研究人员努力的方向。本文运用大涡模拟，考虑人的呼吸过程与热羽流，研究了不同通风速度对气溶胶颗粒物扩散的影响。结果表明，通风速度的增加会使气溶胶颗粒开始排出的时间提前，但同时也会增加其扩散范围，使其扩散的更加均匀。但通风速度从0.75m/s增加至1m/s后，排除颗粒物数量的比例却从18.904%下降到了17.464%。说明通风量增加到一定程度后，气溶胶颗粒的排除效果不会明显的提升。
关键词：通风；气溶胶；大涡模拟
基金项目： 国家自然科学基金项目（51378138）；辽宁省教育厅基础研究项目（LZ2014030）。
0 引言
近年来，随着生活水平的提高，人们对室内空气品质也越来越重视。人体呼出微生物病菌在室内传播，严重影响室内空气品质的同时，危害着人们的身体健康。比如2003年非典和2007年H1N1。通过通风排除室内微生物病菌已被证明是最有效的手段[1]。很多研究人员都已经对气溶胶颗粒物的扩散进行了研究[2-4]。由于此过程影响因素过多，所以一部分研究人员对其影响因素进行了单独的研究。贺启滨等人[5]对不同粒径气溶胶颗粒物的受力大小进行了研究分析，结果表明，不同的粒径所受到的曳力、布朗力、saffman提升力、热泳力都相差较大，小粒径颗粒物的运动受曳力为主导，大粒径颗粒物的运动受重力主导。赵秀国等人[6]对人体呼吸道内气溶胶的沉积进行了研究，结果表明人体上呼吸道气流运动会导致吸入气溶胶在不同位置沉积。C. Voelker等人[7]研究了人体热羽流的气流特征。Chen 等人[8]对蒸发对人体呼出不同粒径的气溶胶的影响进行了研究，结果表明粒径在100μm至200 μm范围内的气溶胶的蒸发对其扩散运动影响较大。
人体呼吸产生的气溶胶扩散到被排除的过程，首先要通过鼻孔呼出的气流带动排出人体，紧接着被人体热羽流夹带着向上运动[9]，这时一部分气溶胶颗粒冲出热羽流的保护，而另一部分随着热羽流向上运动[10]。冲出热羽流的颗粒会被人的肢体行为或者通风所扰动，而向上的颗粒物会扩散的更加剧烈[11]。气溶胶颗粒物向上运动后通风对其扩散占主导地位，比如向下的0.45m/s的通风会完全破坏热羽流进而使气溶胶颗粒物向下运动[12-13]。所以在通风的房间内，不同的通风速度与热羽流对颗粒物的扩散有着怎样的影响，对使用通风控制颗粒物的扩散有着重要的意义。
自计算机技术应用于流体力学，不同的湍流模型被应用于各种场景进行模拟，并得到了较好的结果。随着科学技术的进步和计算机技术的而发展，大涡模拟越来越被广泛的运用到各种场景，并都得到了比传统雷诺时均模型更好的结果[14-16]。而Mathieu Moreau又针对大涡模拟的亚格子模型应用与颗粒物的结果进行了对比。结果表明混合模型 SGS要优于纯湍流粘度模型[17]。
本文采用大涡模拟，针对小型室内的不同通风速度对人体呼气产生的气溶胶颗粒进行了模拟分析，采用了精确的人体模型（由丹麦科技大学通过扫描人体制作），使用了UDF设置了正弦形式的呼吸速度，精确了热羽流等因素。将通风速度设为0.25m/s、0.5m/s、0.75m/s和1m/s四种。
1 模型建立
如图1所示，模拟场所为2m×1m×2.5m（X×Y×Z）的小型房间，采用最常用的侧送方式送风，风速分别为0.25m/s、0.5m/s、0.75m/s和1m/s。风口高度为2m，湍流强度为6.7%，湍流特征长度为3.3×10-2m。由于热羽流是由于空气在不同温度下的密度不同引起的，所以对空气密度采用采用boussinesq假设，膨胀系数为0.0036，单位1/℃。人体使用鼻孔呼吸，呼吸速度采用UDF将其设定为对时间的正弦函数[18]。公式为vt=2.3×sin（1.04×t），其中的正弦周期、速度最大值都是根据文献中设定，鼻孔面积0.56cm2。鼻孔呼出的角度是根据[18]设定。通过对网格数量和时间不长的独立性检验，最终采用105万网格和0.05s的时间步长。
图1 模拟场景
2 控制方程
2.1 大涡湍流模型
湍流流场由不同大小的涡旋组成，而这些涡旋大小的不同也导致了他们在流场中发挥着不同的作用，大的涡旋对平均流动的影响巨大，比如湍流的扩散以及热量、质量、动量的交换，雷诺应力的产生都是通过大尺度涡旋产生的；而小尺度涡旋主要起到的是耗散作用，通过耗散脉动来影响各种变量。大涡模拟将大涡与小涡分开处理，大涡采用直接求解N-S方程的方式，而小涡则使用了通用的亚格子模型。整体控制方程如下：
在求解计算前需要使用滤波函数区分大涡小涡，采用盒式滤波函数，方程如下：
其中Δ为截至宽度，可用下式表示：
对于不可压缩流体，无粘性动量源项的条件下，经过滤波后的LES连续性方程如下：
式中：ui为的是速度的矢量形式；ρ为的是密度（kg/m3）；t为时间。滤波后的动量方程表示如下：
其中：ui为速度在x轴的标量；uj为速度在y轴的标量，压力场p，其中：
而滤波后的能量方程如下：
其中：hs和λ分别为显热焓和热导率。
公式（7）中的子网格焓通量项使用梯度假设进行近似：
其中：μSGS为亚格子粘度，PrSGS为亚格子普朗特数，经常设为常量0.85。
亚格子模型的选择由于许多文献已经证明，WALE（局部涡粘度的壁面自适应）模型在离散相模型上有着更准确的结果[19~20]，所以我们使用WALE亚格子模型，其控制方程如下[21]：
在WALE模型中的涡流粘度通过以下方式建模：
其中Ls和Sdij在wale模型里分别为如下：
其中WALE常数Cw的默认值是0.325。
2.2 离散相模型
离散相模型通过在拉格朗日参考系中写入粒子上的力平衡来预测离散相粒子（或液滴或气泡）的运动轨迹。这个力的平衡等同于粒子的惯性和作用在粒子上的力，由下式表达：
其中F：为附加加速度项FD(u-up)为单位颗粒质量所受到的曳力。FD可由下式表示：
等式（12）中还可以添加其他的附加力，本文所考虑的其他力主要有热泳力，布朗力，saffman提升力。
悬浮在具有温度梯度的气体中的小颗粒在与梯度方向相反的方向上受到的力。这种现象被称为热泳。在公式（12）中，可以在附加加速度项中选择性地包含对粒子的热泳效应：
其中DT, p是热泳力系数，我们将其定义为Talbot所建议的形式如下[22]
其中，Kn=knudsen number =2λ/dp；K =k/kp，其中k为基于平移能量的流体导热率，表示为15/4μR；而kp是颗粒的热导率；Cs 、Ct 、Cm为常数，分别为1.17、2.18、1.14；mp为颗粒物的质量；T是局部流体温度；μ为流体粘度。该表达式默认颗粒物为球形，流体为理想气体。
3 结果与讨论
3.1 呼吸流场
提取不同通风速度的50s时刻的流场绘制如图2所示（a、b、c、d分别代表通风速度为0.25m/s、0.5m/s、0.75m/s和1m/s时的工况）。可以看出气流在人体面前产生了一个漩涡，这是由于四周墙壁对气流的影响。随着气流速度的增加，这个漩涡有了更大的速度，和更明显的轮廓。值得注意的是，当气流速度为0.2m/s时，由于通风速度与人体产生的热羽流速度大小相近，所以通风产生的漩涡中心位于更靠近排风口的偏上的位置，但是与人体周围热羽流产生的速度相比，这个速度还是很小。当气流速度为0.75m/s时，排风口的速度虽然很大，但根据流线我们可以看出，流出的气流大多数来源于顶部，而出风口流出的气流明显向下运动，这是由于出风温度比室内温度小密度大导致的。当通风速度为1m/s时，漩涡周围的速度远远大于漩涡中心的速度。
图2 流场速度与流线图
图3为不同通风风速下鼻孔呼出气流的流线图（a、b、c、d分别代表通风速度为0.25m/s、0.5m/s、0.75m/s和1m/s时的工况），可以看出，当通风速度为0.25m/s时，鼻孔产生的气流路径大部分一致且路程较短，这也是因为通风速度很小导致的。而当通风速度增加到0.5m/s时，流线路程远远大于0.25m/s时的路程，值得注意的是当通风速度增加至1m/s时流线更加分散，与0.5m/s和0.75m/s通风速度下大部分在人体前面的流线不同，会有更多出现在气流死角。
图3 鼻孔出流三维流线图
3.2 颗粒物
图4表示了颗粒物的扩散结果（a、b、c、d分别代表通风速度为0.25m/s、0.5m/s、0.75m/s和1 m/s时的工况），其中不同的颜色代表不同时间产生的气溶胶颗粒，（a）中的颗粒物初始时间大部分向上运动，对比图2（a）我们可以知道这是由于气流速度与向上的热羽流速度相当，而存在时间较久的颗粒物大部分出现在人体周围，刚产生的颗粒物大部分出现在人体上方。0.5m/s的工况下存在时间较旧的颗粒物确相对比较均匀，而新产生的颗粒物再向上运动的过程中受到通风气流的影响逐渐向前运动。值得注意的是0.25m/s、0.5m/s工况中气溶胶颗粒物存在浓度非常低的区域，0.25m/s的情况下出现在人体下方周围，0.5m/s的工况下出现在人体后方，这种现象在其他两个工况中几乎没有，而0.75m/s的工况中明显漩涡中心浓度较小，1m/s的工况中分散的很均匀，几乎没有颗粒物流动的死角。0.75m/s 工况中这样的现象是因为漩涡中心的速度很小，根据图2（c）也可以看出，漩涡中心几乎没有气流经过。
图5 气溶胶颗粒物数量随时间变化图
通过逐时追踪粒子数量，我们可以从图5中得知，当通风速度越大时，粒子开始从排风口排出所需要的时间越短，0.25m/s时大约在35s开始有气溶胶颗粒物排除，而1m/s的通风速度会使气溶胶粒子在20s左右就开始排除，但是就排除效率来说，1m/s的通风速度下，50s时刻残留的气溶胶粒子数却比0.75m/s的通风工况要多。这是可能是由于粒子过于分散造成的，这一点上我们可以看出增加通风量并不一定能提高污染物的排除效率。我们计算50s时刻排除的气溶胶颗粒物的百分比，四种通风风速下分别为为3.08%；12.408%；18.904%和17.464%。
4 结论
通过对四种通风风速对人体呼出气溶胶颗粒物的影响的模拟研究，得出以下结论：
（1）在小室内的通风会使气流形成漩涡，这个漩涡中心的速度通常会很小，会产生气溶胶颗粒物的滞留。
（2）就气溶胶颗粒物的扩散结果来说，通风风速越大，气溶胶的扩散越快也越均匀，这在某些程度上来讲是不利于控制排除他们的；通风风速越小，气溶胶的浓度越高。
（3）通风风速最好不应低于0.5m/s，因为较低的风速不足以主导颗粒物的运动，而是被热羽流也影响着。
（4）就排除效果来说，通风风速增加至1m/s时，效果的提升不再明显，只是在开始排除的时间上有所提前，而在50秒时排除气溶胶颗粒物的百分比反而下降了一点。
就以上结论建议通风时的风速控制在0.75m/s左右，既能增加对气溶胶颗粒物的排除速度，又能最有效率的提高污染物的排除效率。
参考文献
[1] 施珊珊, 纪文静, 赵彬. 不同通风形式下住宅内细颗粒物质量浓度及室内暴露量的模拟及比较[J]. 暖通空调, 2013(12):34–38.
[2] 贺启滨, 朱彤, 董昆,等. 气溶胶在全混通风小室中浓度衰减的试验研究[J]. 安全与环境学报, 2011, 11(3):109–111.
[3] 高乃平, 贺启滨, 李晓萍,等. 人工气候室内呼出气溶胶颗粒物分布的实验研究[J]. 同济大学学报(自然科学版), 2012, 40(11):1680–1685.
[4] 高乃平, 贺启滨, 李晓萍,等. 人工气候室内呼出气溶胶颗粒物分布的实验研究[J]. 同济大学学报(自然科学版), 2012, 40(11):1680–1685.
[5] 贺启滨, 高乃平, 朱彤,等. 人体呼出气溶胶在通风房间中运动的受力分析[J]. 安全与环境学报, 2011, 11(1):242–245.
[6] 赵秀国, 徐新喜, 孙栋,等. 人体上呼吸道气流结构与气溶胶沉积数值模拟[J]. 系统仿真学报, 2012, 24(8):1582–1587.
[7] Voelker C, Maempel S, Kornadt O. Measuring the human body’s microclimate using a thermal manikin.[J]. Indoor Air, 2014, 24(6):567–79.
[8] Chen C, Zhao B. Some questions on dispersion of human exhaled droplets in ventilation room: answers from numerical investigation.[J]. Indoor Air, 2010, 20(2):95–111.
[9] Brohus H, Nielsen P V. Personal exposure in displacement ventilated rooms[J]. Indoor Air, 1996, 6(3): 157–167.
[10] Yang C, Yang X, Zhao B. Person to person droplets transmission characteristics in unidirectional ventilated protective isolation room: The impact of initial droplet size[J]. Building Simulation, 2016, 9(5):597–606.
[11] Salmanzadeh M, Zahedi G, Ahmadi G, et al. Computational modeling of effects of thermal plume adjacent to the body on the indoor airflow and particle transport[J]. Journal of Aerosolence, 2012, 53(2):29–39.
[12] Qian H, Li Y, Nielsen P V, et al. Dispersion of exhaled droplet nuclei in a two‐bed hospital ward with three different ventilation systems[J]. Indoor air, 2006, 16(2): 111–128.
[13] Licina D, Melikov A, Sekhar C, et al. Human convective boundary layer and its interaction with room ventilation flow[J]. Indoor Air, 2015, 25(1):21.
[14] Tian Z F, Tu J Y, Yeoh G H, et al. On the numerical study of contaminant particle concentration in indoor airflow[J]. Building & Environment, 2006, 41(11):1504–1514.
[15] Hu C H, Ohba M, Yoshie R. CFD modelling of unsteady cross ventilation flows using LES[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2008, 96(10–11):1692–1706.
[16] Riber E, Moureau V, García M, et al. Evaluation of numerical strategies for large eddy simulation of particulate two-phase recirculating flows[J]. Journal of Computational Physics, 2009, 228(2):539–564.
[17] Moreau M, Simonin O, Bédat B. Development of Gas-Particle Euler-Euler LES Approach: A Priori Analysis of Particle Sub-Grid Models in Homogeneous Isotropic Turbulence[J]. Flow Turbulence & Combustion, 2010, 84(2):295.
[18] GUPTA J K, Lin C H, Chen Q. Characterizing exhaled airflow from breathing and talking[J]. Indoor air, 2010, 20(1): 31–39.
[19] 王兵, 张会强, 王希麟,等. 不同亚格子模式在后台阶湍流流动大涡模拟中的应用[J]. 工程热物理学报, 2003, 24(1):157–160.
[20] 俞建阳, 王若玉, 陈浮,等. 不同亚格子模型的对比分析及其运用[J]. 工程热物理学报, 2016, V37(11):2311–2318.
[21] Nicoud F, Ducros F. Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[J]. Flow Turbulence & Combustion, 1999, 62(3):183–200.
[22] Talbot L, Cheng R K, Schefer R W, et al. Thermophoresis of particles in a heated boundary layer[J]. Journal of Fluid Mechanics, 1980, 101(101):737–758.
备注：本文收录于《建筑环境与能源》2019年1-2月刊总第18期。
版权归论文作者所有,任何形式转载请联系作者。
