据统计，人们每天90%以上的时间在室内环境中度过^［1］，因此室内空气质量的优劣与人体健康息息相关。室内污染物包括气态污染物和颗粒物，其中气溶胶颗粒物粒径通常小于10 μm，可以通过呼吸作用进入人体。人们长期接触并暴露在高浓度颗粒物污染物的环境下，会对自身健康造成不利影响，包括呼吸系统疾病、心血管疾病等^［2–5］。

研究表明^［7–8］，冬季大气环境中颗粒物污染物较其他季节更为严重。室内颗粒物污染物来源包括室外渗透、室内污染源散发和壁面沉积再悬浮等^［9］。对于室内污染源散发和室外渗透已有很多研究。Li等^［10］研究了5种通风方式下人体呼出的气溶胶颗粒物的扩散和去除率，结果表明地板通风时颗粒物污染物去除率最高。Cetin等^［11］通过实验和数值模拟研究了不同高度处点源释放的颗粒物在置换通风和地暖作用下的扩散和沉积特性，结果表明颗粒物扩散与出风口位置相关。Taheri等^［12］研究了在不同送风角度下地板通风对室内颗粒物沉积和去除的影响，结果表明进气角为90°时颗粒物去除率最高达到58%。Alotaibi等^［13］研究了地板再悬浮颗粒物和呼出气溶胶颗粒物在不同通风条件下的扩散和分布特性，并提出在置换通风和混合通风之间采取切换操作的建议。Dehghan等^［14］研究了不同采暖环境下颗粒物释放位置对其扩散和分布的影响，结果表明颗粒物大多分布在墙壁和天花板或者停留在近地面。Zhou等^［15］研究了不同风速和地板温度对颗粒物去除率的影响，结果表明地板温度在308 K时颗粒物清除时间可缩短15%。Zhang等^［16］研究了厨房烹饪产生的颗粒物在供暖通风房间内的扩散分布特性，结果表明在地暖供热条件下，室内颗粒物分布均匀。

此外，再悬浮颗粒物也是室内污染物的重要来源，其中粒径小于10 μm的颗粒物占再悬浮颗粒物的30%（质量分数）以上^［9］。颗粒物再悬浮过程分为两个阶段：第一阶段颗粒物从壁面分离进入边界层；第二阶段颗粒物离开边界层进入流场^［17］。而现有研究主要集中在第一阶段，包括在行走过程中不同粒径颗粒物的再悬浮特性，即再悬浮率及分布情况^{［18–20］}。Assaad等^［21］建立了再悬浮颗粒物扩散计算流体力学（CFD）模型，研究发现，扰动过程的前25 s，颗粒物再悬浮率快速增加，且不同材料的沉积表面对颗粒物再悬浮特性影响较大。

除了室内人员活动引起的颗粒物再悬浮，冬季辐射地暖系统也会引起地板处不同粒径颗粒物再悬浮^［13］，因此需要新风系统以保证室内环境质量^［22-23］。同时，还需兼顾能耗与热舒适性。混合通风可使室内空气充分混合，允许送风与室内空气间存在较大温差而不引起冷风感^［24］；辐射地暖系统可以提供理想的空气垂直温度梯度，使人体的手脚温差较小，从而提供较好的室内热舒适性^［25］。基于此，本文建立辐射地暖房间的混合通风模式，通过实验和数值模拟，探究不同阿基米德数Ar下再悬浮颗粒物沉积和去除特性。

1 研究方法 1.1 物理模型及相似理论

本文中房间物理模型 (3 m × 3 m × 3 m) 如图1所示。该模型耦合了混合通风和地板采暖模式，其通风口尺寸为0.6 m × 0.3 m。图2为采用缩尺比1∶3建立的环境实验舱实物。环境实验舱主体尺寸（长 × 宽 × 高）为1 m × 1 m × 1 m，舱体正面中心设有0.6 m × 0.6 m的可视化监视窗；舱体上方中间设有0.4 m × 0.4 m的可封闭开口。进风管和回风管长均为0.4 m，以保证进、出风口充分发展。通过变频风机送风和底部加热板供热。经过调试，管道内送风均匀，湍流度较小；舱体底部温差较小，控制在±5 ℃以内。距离进口侧和出口侧壁面0.2 m处设有2个直径为2 cm的圆形测量孔。利用TSI热线风速仪、多点温度测量仪和红外相机等对舱体内部速度和温度进行测量。测量线L₁、L₂及混合通风配置如图3所示。

基于相似理论，对于非等温流动，缩尺模型与全尺度模型需满足雷诺数Re、格拉晓夫数Gr和阿基米德数Ar相似才能确保两模型相似。由分析可知，平均Re超过4 000时，两模型流场结构基本相似，即与Re无关；Gr大于10⁹时，靠近加热表面的流动为完全湍流^［26］。因此，只需保证Ar相似即可。Ar定义为

$ \qquad {A}{r}=\frac{g{\alpha }_{\mathrm{v}}\Delta tL}{{u}^{2}} $

(1)

式中：g为重力加速度，m/s²；α_v为空气体积膨胀系数，1/K；Δt为房间地面温度与送风温度的差值，K；L为房间特征尺度，m；u为送风速度，m/s。

1.2 数值模型

采用雷诺平均Navier−Stokes方程对室内通风及热浮力引起的湍流与扩散进行描述^［27］。湍流模型采用RNG k−ε模型^［28］，其具有更好地处理高应变率和可承受流线弯曲程度较大流动的能力。控制方程为

连续性方程

$ \qquad \frac{{\partial \left( {\rho {u_j}} \right)}}{{\partial {x_j}}} = 0 $

(2)

动量方程

$ \qquad \begin{split} \frac{{\partial \left( {\rho {u_i}{u_j}} \right)}}{{\partial {x_j}}} =& - \frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial {u_i}}}{{\partial {x_j}}} - \rho \overline {{u_i}{u_j}} } \right) +\\& \rho \beta {g_i}\left( {{T_{{\text{ref}}}} - T} \right) + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) \end{split} $

(3)

能量方程

$ \qquad \frac{{\partial \left( {\rho h{u_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left( {\frac{\lambda }{{{c_p}}}\frac{{\partial h}}{{\partial {x_j}}} - \rho \overline {{u_j}h} } \right) + {S_{\mathrm{h}}} $

(4)

RNG k−ε湍流模型

$ \frac{\partial }{{\partial {x_j}}}\left( {\rho k{u_j}} \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _k}{\mu _{{\text{eff}}}}\frac{{\partial k}}{{\partial {x_j}}}} \right) + {G_{\mathrm{k}}} + {G_{\mathrm{b}}} - \rho \varepsilon - {Y_{\mathrm{M}}} + {S_{{k}}} $

(5)

$ \begin{split}\frac{\partial }{\partial {x}_{j}}\left(\rho \varepsilon {u}_{j}\right)=&\frac{\partial }{\partial {x}_{j}}\left({\alpha }_{\varepsilon }{\mu }_{\text{eff}}\frac{\partial \varepsilon }{\partial {x}_{j}}\right) + {C}_{1\varepsilon}\frac{\varepsilon }{k}({G}_{{\mathrm{k}}} + \\&{C}_{3\varepsilon }{G}_{{\mathrm{b}}})- {C}_{2\varepsilon }\rho \frac{{\varepsilon }^{2}}{k}-{R}_{\varepsilon } + {S}_{\varepsilon } \end{split} $

(6)

式中：x_i、x_j为笛卡尔坐标；u_i、u_j分别为待求解时均速度和流体时均速度，m/s；ρ为空气密度，kg/m³；p为空气压力，Pa；μ为空气动力黏度，kg/（m·s）；β为空气热膨胀系数，1/K；T_ref为参考温度，K；T为空气温度，K；g_i为i方向的重力加速度，m/s²；h为空气定压比焓，J/kg；S_h为热源源项，W/（m ³·s）；λ为空气热导率，W/（m·K）；c_p为空气定压比热容，J/（kg·K）；k为湍动能；ε为耗散率；G_k、G_b分别为受平均速度梯度、浮力影响产生的湍动能；Y_M为在可压缩流动中脉动膨胀对总耗散率的影响；R_ε为ε方程中的附加项；C_1ε、C_2ε、C_3ε均为常量；α_k和α_ε分别为k和ε的有效普朗特数的倒数；S_k和S_ε分别为自定义湍动能和湍动耗散率的源项；$ \rho \overline{u_i u_j} $、$ \rho \overline{u_jh} $分别为雷诺应力和雷诺传热量；$ \mu_{_{{\mathrm{eff}}}} $为RNG模型的修正系数。

本文选取适用范围最广的离散坐标（DO）辐射模型^［29］，并采用离散相模型（DPM）追踪单个颗粒物在室内的运动轨迹。拉格朗日跟踪粒子轨迹方程为

$ \qquad \frac{{{\mathrm{d}}{u_{\mathrm{p}}}}}{{{\mathrm{d}}t}} = {F_{\mathrm{D}}}\left( {{u_j} - {u_{\mathrm{p}}}} \right) + \frac{{g\left( {{\rho _{\mathrm{p}}} - \rho } \right)}}{{{\rho _{\mathrm{p}}}}} + {F_{\mathrm{a}}} $

(7)

$ \qquad{F_{\mathrm{D}}} = \frac{{18\mu }\;{{{C_{\mathrm{D}}}Re}}}{{{24\rho _{\mathrm{p}}}d_{\mathrm{p}}^2}}$

(8)

式中：u_p为颗粒物速度，m/s；ρ_p为颗粒物密度，kg/m³；d_p为颗粒物直径，μm；C_D为阻力系数；$ F_{\mathrm{D}} $为弛豫时间的倒数，s⁻¹；F_a为附加力，m/s²；t为粒子在该网格中的停留时间，s；式（7）等号右边的项分别表示曳力、浮力和附加力。

本文中附加力只考虑萨夫曼升力、热泳力和布朗力。计算时空气速度为平均速度，因此采用随机游走扩散（DRW）模型模拟湍流脉动速度对颗粒物运动的影响。该模型假设脉动速度服从高斯概率分布，且各向同性。脉动速度分量$ u_i^\prime $为

$\qquad u_i^\prime = {\zeta _i}\sqrt {\bar u_i^{\prime 2}} = {\zeta _i}\sqrt {\frac{2k}{3}} $

(9)

式中：ξ_i为正态分布随机数；$ \bar{u}_i^\prime $为脉动平均速度分量。

采用质点网格法计算污染物质量浓度。质点网格法以计算网格为基础，将质量浓度和轨迹相关联，使用多个控制网格对包含颗粒物的空间进行离散化，最后根据颗粒物在控制空间内的停留时间计算颗粒物质量浓度^[30]。

在地板上方0.01 m平面处释放颗粒物，模拟地板加热引起的热对流所导致的进入流场的再悬浮颗粒物。2 500个释放点均匀分布于该平面，两点之间距离为0.02 m，每个释放点释放40个颗粒物，共计10⁵个颗粒物。再悬浮颗粒物密度为870 kg/m^3[31］，直径为2.5 μm，初始速度均为0，温度为300 K，释放率为2.5 μg/s^［32］。

边界条件设置如表1所示。

采用有限体积法对方程进行离散，压力−速度耦合方程采用SIMPLE算法进行求解，控制方程离散格式均采用二阶迎风格式，收敛条件为各变量的相对计算残差小于10⁻⁶。

1.3 工况设置及网格敏感性分析

混合通风是通过房间上部送风、下部回风来稀释室内空气，以改善空气质量，因此设置通风口、回风口位置分别为同侧和异侧，如图3所示。表2为混合通风方式下的参数，其中：T_in为送风温度，K；T_f为地板供热温度，K；ΔT为地板供热温度和送风温度的差值，K；$ Pr $为普朗特数。针对混合通风方式，研究6种Ar下室内流场及颗粒物的分布规律。

异侧通风工况下送风速度为0.5 m/s^［15］。采用六面体结构化网格进行区域离散，选择标准壁面函数来处理近壁面层。第一层网格选取在距离壁面0.005 m处。网格敏感性分析时采用网格精度不同的4套网格，网格数量分别为25万、65万、158万和398万。选取L₁处速度分布进行敏感性分析。图4为网格敏感性分析及网格设置。由图4(a)可知，网格数量为25万和65万时的计算结果相比于158万和398万时的偏大，而网格数量为158万和398万时的计算结果几乎一致。这说明计算结果在网格数量大于158万后趋于一致。158万网格划分如图4(b)所示。

采用速度均方差$ {\varepsilon _{\text{v}}} $作为网格独立性的定量标准^［33］，即

$ \qquad {\varepsilon _{\text{v}}} = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {\left(\frac{{{V_{{{i,{{n}} - 1}}}} - {V_{i,{n}}}}}{u}\right)} } $

(10)

式中：V_i,n−1、V_i,n分别为编号为n−1、n的网格的速度，m/s；N为检查采样点的数量。

随着网格不断加密，$ {\varepsilon _{\text{v}}} $逐渐减小，且网格数量由158万增加至398万时其变化很小。因此，本文选用数量为158万的网格进行计算。

2 结果与讨论 2.1 模型验证

本文中环境实验舱工况分为等温工况和非等温工况。等温工况指实验过程中未对地板进行加热，地板与送风气流的温差为0；非等温工况指在实验过程中对地板进行加热，且地板与送风气流保持稳定的温差。

图5(a)、(b)分别为Ar = 0及异侧通风工况下，L₁和L₂处速度的模拟结果和实验结果对比。由图可知，实验结果和模拟结果一致性较好。这表明所采用的数值模型可以较好地预测通风房间内速度分布。

图5(c)、(d)分别为Ar = 8.04及异侧通风工况下，L₁处速度和无量纲温度θ的模拟结果和实验结果对比。θ的计算式为

$\qquad \theta = \frac{{{T_{\text{c}}} - {T_{{\text{in}}}}}}{{{T_{\text{f}}} - {T_{{\text{in}}}}}} $

(11)

式中，T_c为测点温度，K。

从图5可知，非等温工况下速度最大值出现的位置相较于等温工况下出现的位置低，且与无量纲温度最小值出现的位置（高度为0.7 ～ 0.8 m）一致。这表明通过数值方法可以模拟出房间内在速度高值区影响下形成的温度低值区，即“空气湖”效应^［34］。此外，L₁处实验测得的速度和无量纲温度与模拟结果总体一致性较好。这表明该模型可较精确地预测室内混合对流情况。综上，将RNG k−ε湍流模型应用于研究地暖通风房间内的流动与扩散是可靠的。

2.2 流场随Ar的变化

图6为异侧和同侧通风工况、不同Ar下轴面流线图及无量纲速度云图，其中：V/u为无量纲速度；V为室内风速，m/s。等温工况下，由于康达效应^［31］和壁面阻挡，入射气流在房间上方形成一个逆时针旋涡。回风口位置会影响室内流场：同侧通风工况下房间中部形成顺时针旋涡，而异侧通风工况下气流经回风口排出，无旋涡生成。

加热地板后，送风气流由于密度差向下偏折，在受到阻挡后向四周扩散，因此送风落点从壁面转移到地板，临界点为Ar = 6.03。Ar的增加导致等温工况下的旋涡消失并形成新的旋涡：同侧通风工况下在房间右侧形成逆时针旋涡，异侧通风工况下在房间左侧形成顺时针旋涡。

2.3 温度场随Ar的变化

图7为异侧和同侧通风工况、不同Ar下无量纲温度云图。Ar增加时，低温区从壁面向地板迁移。这是由于冷气流因地暖作用产生了密度差。回风口位置对室内温度分布影响较大。异侧通风工况下，Ar < 6.03时，近地面温度梯度明显，房间中上部温度均衡；Ar ≥ 6.03时，室内温度基本一致，且由于送风气流的分割作用，温度低值区迁移到回风口远处。Ar = 2.01时，旋涡影响温度分布，从而形成稳定的等温区域；Ar > 2.01时，旋涡对室内温度影响较小，靠近送风气流的区域温度偏低，靠近地暖的区域温度偏高，其他区域温度接近。

图8为异侧和同侧通风工况下L₁和L₂无量纲温度梯度分布。由图可知：低温点随Ar增加向下偏移，这与送风气流向下偏折有关；同侧通风工况下，Ar > 4.02时各Ar下的温度分布几乎重合，表明该工况下较异侧通风工况下更能维持室内温度稳定。

2.4 颗粒物质量浓度分布随Ar的变化

图9为异侧和同侧通风工况、不同Ar下混合通风地暖房间2.5 μm再悬浮颗粒物质量浓度云图，其中C为颗粒质量浓度，μg/m³。Ar为0时，异侧通风工况下颗粒物质量浓度低值区主要在房间上方，其余区域为颗粒物质量浓度高值区。这是受气流输送的稀释和回风口过早排出气流的影响所致。同侧通风工况下，再悬浮颗粒物主要集中在地板附近，其余区域颗粒物质量浓度较低。这是受回风口位置影响形成的大旋涡稀释了颗粒物质量浓度所致。

当加热地板时，浮升力和送风气流导致颗粒物质量浓度分布均匀，且Ar的增加使得回风口位置对室内颗粒物质量浓度分布无影响。异侧通风工况下，Ar ≥ 8.04时，颗粒物质量浓度分布均匀；Ar较小时，房间左侧存在颗粒物质量浓度高值区，这是旋涡效应所致。同侧通风工况下，颗粒物质量浓度高值区在地板附近，且Ar增大时，颗粒物质量浓度高值区缩小，这是因为送风落点逐渐向地板转移，从而稀释了颗粒物质量浓度；Ar为4.02时，旋涡消失造成气流流向室内，因此房间左侧出现颗粒物质量浓度高值区。

异侧和同侧通风工况下不同平面再悬浮颗粒物沉积率如图10所示。沉积率的定义是在某表面上的颗粒物数与从边界层释放的颗粒物总数之比。Ar为0时，地板上颗粒物沉积率最高，因此影响颗粒物沉积的主要因素是其自身重力；而随着Ar的增加，由于热浮升力的影响，地板上的颗粒物沉积率呈下降趋势，地板和其他壁面的颗粒物沉积率呈上升趋势。

再悬浮颗粒物去除率如图11所示。去除率的定义为从计算域逸出的颗粒物数与从边界层释放的颗粒物总数之比。Ar为6.03是临界点，Ar小于6.03时，同侧通风工况下去除率高于异侧通风工况下的值，而Ar超过6.03时则反之。同侧通风工况下去除率随Ar的增加而降低，这与流场在Ar为6.03时送风落点迁移至地面相对应，此时房间右侧形成逆时针旋涡并将颗粒物卷吸入，使其难以排出。异侧通风工况下，受热浮升力和旋涡的共同影响，在Ar为2.01时颗粒物去除率最大，为31%。

行人呼吸高度处再悬浮颗粒物平均质量浓度如图12所示。行人呼吸高度的定义是高度为1.5 m（对应模型：0.5 m）。由图可知，不同Ar下，异侧通风工况下行人呼吸高度处颗粒物质量浓度高于同侧通风工况下的值。虽然受热浮升力的影响质量浓度逐渐降低，但最低质量浓度仍高于同侧通风工况下的值，且最低值在Ar为2.01时，为21.56 μg/m³。

综上，选择同侧通风工况且Ar为2.01时，室内环境处于较适宜的状态。

3 结论

采用经实验验证的RNG k−ε湍流模型，研究混合通风和地板供热条件下室内再悬浮颗粒物沉积和去除特性。主要结论为：

(1) 混合通风下随着Ar的增加，送风气流都会向下弯折而改变室内流场。Ar为6.03时，送风落点从壁面迁移到地板。混合通风下，受Ar的影响，分别会在房间右侧形成逆时针旋涡，在房间左侧形成顺时针旋涡。

(2) 混合通风下室内温度场会受Ar的影响。异侧通风工况下，靠近送风气流的区域室内温度较低，远离送风气流的区域温度则较高，且存在明显的温度梯度；而同侧通风工况下，Ar大于4.02时，室内温度均匀分布，无明显的温度梯度。

(3) 加热地板后，颗粒物质量浓度分布逐渐变得均匀，无明显的高值区，因此Ar增加后回风口位置对室内再悬浮颗粒物质量浓度分布无影响。同侧通风工况下随Ar增加再悬浮颗粒物去除率降低，异侧通风工况下去除率最高值在Ar为2.01时。异侧通风工况下随Ar增加行人呼吸高度处颗粒物质量浓度降低，但最低质量浓度仍高于同侧通风工况下的值，且最低值在Ar为2.01时，为21.56 μg/m³。综上，选择同侧通风工况且Ar为2.01时，室内空气质量较为适宜。