自然通风由于具有控制室内污染物浓度水平、改善室内空气品质、调节室内热环境等功能，并且具有节能、环保等优点，越来越受到人们的重视。文献[1]数值分析了通风量、窗户位置、开窗面积以及室外风速之间的关系。文献[2]对深圳2个典型户型在风攻角变化时房间的自然通风进行了数值研究。文献[3]实验研究了窗墙比、建筑前后风压差、进出风口面积比以及进出风口位置等条件对通风量的影响。文献[4]利用CFD软件数值模拟了自然通风条件下不同排风口布置方式对热压通风的影响。文献[5]对教室的自然通风效果进行了能耗分析。文献[6]对意大利传统建筑的自然通风节能潜力进行了调查研究，发现在意大利3种气候区中，自然通风相对于空调系统分别使冷负荷降低了41%，46%和 52%。文献[7]调查分析了温带国家丹麦7月和8月建筑房屋的被动散热方式及自然通风的节能性能，得出了使用自然被动通风取代原先90%的机械通风时间仍满足室内热舒适性要求的结论。文献[8]对德国法兰克福办公楼的通风监测数据表明：冬季通风时间超过全天时间的10%~25%，会使冬季的能量消耗增加。文献[9]对印度复杂气候条件下自然通风效果进行了仿真模拟。以上国内外研究工作主要是针对夏季或者过渡季展开的。而寒冷地区冬季气温低、昼夜温度波动大，人们通常利用长时间关闭窗户来维持室内温度的恒定。这样，人体、地板、墙壁、家具等散发的污染物就会聚集而使得室内空气质量严重恶化。如何满足居住者对室内采暖温度的需求和保证符合室内空气质量标准要求的室内空气品质因此成为暖通界研究的热点问题。笔者就此问题展开研究，借助CFD软件、采用Realizable k-ε 湍流模型，以兰州某民用住宅为例，对室内流场、温度场、污染物浓度场、舒适性进行了数值分析，旨在探寻寒冷地区冬季自然通风的合理方式。

所研究住宅为3室2厅，位于8层建筑的第4层，结构尺寸为：x×y×z=10.5 m×13.2 m×2.9 m，平面示意图如图 1所示。其中，1#为客厅与餐厅，2#、3#、4#为卧室。外墙总传热系数K=0.46 W/(m²·K)。外窗传热系数K_c=2.2 W/(m²·K)。室内为地板供暖。通过改变窗户C1、C2的开度来控制进入室内的新风量，通过改变窗户C3、C4的开度来控制室内的排风量。冬季通风量过大会造成室温下降过快，这里取每一窗户的最小可开启面积为：宽×高=0.1 m×1.5 m(下文中简称为A)；最大可开启面积为：宽×高=0.2 m×1.5 m(下文中简称为2A)。

图1(Fig. 1)

图1 计算房间平面示意图 Fig. 1 Physical model of the calculating room

为简化计算，做如下假设：

1)室内空气不可压缩，密度的变化采用 Boussinesq假设，忽略流体粘性力作功而引起的耗散热；

2)流动为非稳态，流态为湍流；

3)室内空气假定为辐射透明介质，不参与辐射换热，忽略各内壁面之间辐射换热；

4)不考虑冷风渗透影响，室内气密性良好。

采用标准k-ε模型、RNG k-ε 模型以及Realizable k-ε 模型分别对文献[10]的通风过程进行了数值计算，并将数值结果与文献[10]的实验结果进行了比对。结果表明：Realizable k-ε模型计算结果和实验结果吻合度最好。所以，选用Realizable k-ε模型进行模拟。室内空气湍流流动与组分输运方程具有相同的形式^[11]，对于非稳态传热问题，通用控制方程形式均可表示为

式中：φ分别代表连续性方程及动量方程中的速度在x轴向的分量u、在y轴向的分量v、在z轴向的分量w、能量方程中的流体温度T、湍流动能方程中的湍流动能k、湍流动能耗散率方程中的湍流动能耗散率ε、组分输运方程中的污染物体积分数τ_p；Γ_φ分别代表相应方程中的有效扩散系数0、μ+μ_t、Pr+μ_t/Pr_t、μ+μ_t/σ_k、μ+μ_t/σ_ε、μ+μ_t/σ_τ；S_φ分别代表相应方程中的源项，组分输运方程中的源项为空气密度ρ，其他方程中源项具体表达式见文献[12]。μ为分子粘性系数；μ_t为湍流粘性系数，μ_t=c_μρk²/ε，c_μ为经验常数；Pr_t为能量方程中的湍流普朗特数；σ_k为湍流动能方程中的普朗特数；σ_ε为湍流动能耗散率方程中的普朗特数；σ_τ为组分输运方程中的普朗特数。

方程中各系数取值^[12]：c_μ=0.09,σ_k=1.0,σ_ε=1.2,Pr_t=1.0,σ_τ=1.0,c₂=1.9，c₁=max{0.43,$\tilde \eta $/(5+$\tilde \eta $)}，$\tilde \eta $=Sk/ε，S=(2S_i,jS_i,j)^1/2，S_i,j=($\partial $u_i/$\partial $x_j+$\partial $u_j/$\partial $x_i)/2。

寒冷地区冬季气温低，一年内较冷的天数所占比重大。以黄河上游的内陆城市兰州为例，属于中温带大陆性季风气候，全年降水少，太阳辐射潜力大，全年日照时数约为2 600 h。1月平均气温为-6.9 ℃，极端最低温为-21.7 ℃。计算中，根据室外逐时温度和太阳辐射强度得到室外综合温度，再以气固耦合的方式并考虑围护结构的衰减和延迟效应，获得通过外墙的热流密度，以此热流密度作为数值模拟室内热环境的热边界条件。所以，兰州冬季室外的低温条件决定了建筑自然通风模式与其他温和地区自然通风模式必然存在显著差异。

室内污染物的种类繁多。其中，C₇H₈(甲苯)是引发“建筑综合症”以及人体癌变的诱因之一，对人体健康极为不利。国内外对室内C₇H₈污染非常重视，故笔者选取C₇H₈为室内的代表性污染物。

1)入口边界：综合IWEC(美国国家气象数据中心)和CSWD(清华大学和中国气象局)所提供兰州冬季气象参数，兰州在冬季以北风为主导风向，故以北外窗C1、C2为计算区域的流动入口边界；11月的月平均风速最小，两种统计结果均为0.2 m/s，3月份风速最大，2种统计结果均为1.2 m/s，故风速取1.2 m/s和0.4 m/s2个值，对应湍流强度分别为4.92%和5.64%(开口面积为A，水力直径为0.187 5 m)以及4.54%和5.21%(开口面积为2A，水力直径为0.352 9 m)。进风温度取0 ℃(即273.15 K)。进风口空气成分及体积分数为：O₂为21.0%，CO₂为0.03%，H₂O为0.04%，其余均为N₂。

2)出口边界：排风口与进风口距离较远，通风路径长，可认为气体流动已充分发展，故排风口设为自由出流边界(Outflow)^{[13, 14, 15]}。

3)壁面边界：气固交界面的空气流动取速度无滑移条件。

4)热边界设置：采用非稳定传热方法得到不同时刻通过围护结构的热流密度，以通风初始时刻11:00时的热流密度为依据进行参数设置。具体数值为：东墙-8.26 W/m²，西墙-8.44 W/m²，南墙-8.00 W/m²，北墙-8.65 W/m²，南窗-20.68 W/m²，北窗-33.29 W/m²，负号表示室外向室内传递的热流为负值。室内采暖热负荷均由地板提供，设为定热流边界条件，根据文献[16]设定地板的热流密度为45 W/m²。内墙、天花板均设为绝热边界。

5)污染源条件设置：来自地板的污染物C₇H₈是沿地板外法线方向向上均匀散发的。根据文献[17]对室内C₇H₈允许浓度规定，结合散发时长以及房间体积，可求得散发速率为5×10^-9kg/s。选用层流模型模拟了污染物的扩散过程(限于篇幅，另文讨论)。对计算房间密闭情况下污染物散发1、6、12、15 h后的浓度分布进行分析并结合人们冬季开关窗的习惯，确定以密闭15 h后的室内污染物浓度场为自然通风数值模拟的初始条件。

在三维直角坐标系下采用6面体结构网格对计算空间进行离散。考虑进、排风口处流动参数的高梯度变化，为获得这些区域流动参数的准确信息，对其进行了局部加密。利用本文构造的3套网格(网格数分别为520 590、672 138和868 434)进行了数值试验，所得(x=8.98 m，y=3 m)处的速度和温度计算结果均吻合得很好，说明所得解是网格独立性的。考虑计算的经济性，决定选用672 138作为后续计算的网格数。

采用有限体积法对控制方程(1)进行离散；应用SIMPLE算法求解速度/压力耦合问题^[18]；梯度项方程采用Green-Gauss Cell Based格式离散；动量方程、能量方程、湍流动能方程、湍流动能耗散率方程、组分输运方程均采用二阶迎风格式离散。

进风温度为0 ℃，进风速度为1.2 m/s和0.4 m/s，与不同进风口开启面积、排风口开启面积组合，形成12个工况，如表 1所示。

表1(Table 1)

表1 自然通风模拟工况Table 1 Simulation cases of natural ventilation 工况 进风温度 /℃ 进风速度/ (m·s-1) 进风窗口编号及开启面积 排风窗口编号及开启面积 1 2 3 4 5 6 7 8 9 10 11 12 0 0 0 0 0 0 0 0 0 0 0 0 1.2 1.2 1.2 1.2 1.2 1.2 0.4 0.4 0.4 0.4 0.4 0.4 C 1:A C 1:A C 1:2 A C 1:A C 2:A C 1:A C 2:A C 1:A C 2:A C 1:A C 1:A C 1:2 A C 1:A C 2:A C 1:A C 2:A C 1:A C 2:A C 3:A C 3:A C 4:A C 3:A C 3:A C 4:A C 3:A C 4:2 A C 3:A C 3:A C 3:A C 4:A C 3:A C 3:A C 4:A C 3:A C 4:2 A C 3:A

表1 自然通风模拟工况 Table 1 Simulation cases of natural ventilation

为了确定合理的通风工况，规定同时满足：通风1 h后距地面1.1 m高度水平面(规定为工作区高度)上温度的平均值不低于16.0 ℃和该平面上污染物浓度不高于文献[17]的要求值2.17×10^-9 kmol/m³,则认为该工况是符合通风要求的。据此，当进风温度为0 ℃、进风速度为1.2 m/s时，工况1~2大约在通风22 min时温度下降到室内采暖要求的最低温度16.0 ℃^[17]，而C₇H₈要降到室内空气质量标准要求上限值2.17×10^-9 kmol/m³以下则需要29 min，这说明在通风时间上不能同时满足最低采暖温度和最高污染物允许浓度的要求。工况3~6更不能满足自然通风的要求(此处不再给出计算结果)。

当进风温度为0 ℃、进风速度为0.4 m/s时，距地面1.1 m高度水平面上温度的平均值和该平面上污染物浓度平均值随时间变化趋势如图 2和图 3所示。由图 2可以看出，不同通风工况下室内温度随时间的下降速度不同。但通风60 min，工况7~12均能满足室内采暖温度需要。图 3中工况9~12对应的C₇H₈平均浓度变化趋势非常接近，并且在51 min时接近文献[17]规定的上限值。综合室内采暖温度要求及污染物浓度上限值两个因素，确定工况9~12的通风时长为51 min。而工况7~8的污染物浓度随时间下降趋势基本一致但过于缓慢，不利于健康。

图2(Fig. 2)

图2 平均温度变化曲线 Fig. 2 Average temperature profile

图3(Fig. 3)

图3 C7H8平均浓度变化曲线 Fig. 3 Average concentration profile of C7H8

综上，进风温度为0 ℃、进风速度为0.4 m/s时的工况9~12为可行的通风工况。但根据对室内竖向温差的统计结果，工况9中地面上方0.1 m到1.1 m间的竖向温差大于3.0 ℃，不满足ISO7730对热舒适的要求。所以，后续分析只列举工况10~12的计算结果。

计算中，同时满足以下条件，认为计算已收敛：1)连续性方程及动量方程残差设为10^-3；2)能量方程求解残差设为10^-6；3)监视位置气流参数不再波动；4)进出口流体质量守恒。

通风效率表示外界进风的排污或排热能力。能维持稳态散发情况下较低的室内污染物浓度或能较快地降低室内污染物浓度，其通风效率就高。评价室内污染物浓度变化一般利用工作区的通风效率来反映污染物浓度水平，其定义式如式(2)所示。

式中：η为通风效率；C_g为工作区平均浓度，kmol/m³；C_P为排风口浓度，kmol/m³；C_S为进风口浓度，kmol/m³。

对工况10~12，统计不同通风时长Z=1.1 m平面上的C₇H₈平均浓度作为工作区污染物平均浓度，同样的方法得到排风口对应时长的污染物平均浓度C_P，进入室内的新风中认为不含有C₇H₈，所以，取C_S=0。据此得到3种工况下的通风效率，如图 4所示。

图4(Fig. 4)

图4 通风效率变化曲线(Z=1.1 m) Fig. 4 Ventilation efficiency profile (Z=1.1 m)

可以看出，3种工况的通风效率均在0.95~1.0之间波动。通风15 min后通风效率均呈下降趋势。这是因为刚开始通风时，房间内C₇H₈分布较均匀，随着通风的持续进行，C₇H₈浓度出现分层现象，Z=1.1 m水平面上的C₇H₈面平均浓度下降趋势因此变缓，从而导致排风口的C₇H₈平均浓度下降速率比Z=1.1 m水平面上的C₇H₈面平均浓度下降速率高。但总体上，工况12的通风效率要高于工况10和工况11的通风效率，这与不同的开窗方式密切相关。

由以上分析可知，“双进口、单出口”的开窗方式对于污染物的排除效率较高。这也是自然通风中流场、温度场与污染物浓度场耦合作用的结果。限于篇幅，笔者只给出速度场和温度场的模拟结果，浓度场变化另文讨论。

为了同时保证室内采暖温度要求和污染物浓度不超标，2.3.1 节所确定的工况10~12通风时长为51 min。下面给出这3个工况通风51 min时的模拟结果。

图 5为不同工况时客厅进风窗口中轴线(X=1.15 m)纵截面上速度场分布。可以看出，3个工况表现出了相似的特征：温度较低的室外新风进入房间后，沿窗台下内墙壁形成下降气流，与地板附近的热空气相遇后便朝远离窗口的方向运动，与被加热后向北运动的热气流相遇，形成了扁平状的漩涡。在靠近南墙附近区域，流动变得更加紊乱。但工况12靠近窗户附近的回流速度平均水平要比工况10和工况11的低。

图5(Fig. 5)

图5 X=1.15m截面速度矢量 Fig. 5 Velocity vectors in X=1.15 m section

图 6为不同工况时X=7.93 m(2#房间和3#房间两房门中轴线)纵截面上的速度矢量分布。可以看出，3个工况表现出了相似的特征：温度较低的室外新风进入房间后，地面附近的冷气湖自北朝南运动，流经房门时由于过流断面缩小而加速，遇到3#房间南墙时在浮升力作用下形成上升气流，在房间半高度处迅速降落并朝北向回流，再一次流经房门时受下部冷气流的顶托作用，加速后朝2#房间顶部运动，形成了热羽流，工况12的热羽流影响范围和强度要比工况10和工况11的弱。这样，就在2#房间上部形成了一个范围不同的热滞留区，这是污染物容易聚集的位置。

图6(Fig. 6)

图6 X=7.93 m截面速度矢量 Fig. 6 Velocity vectors in X=7.93 m section

图 7为不同工况时客厅进风窗口中轴线(X=1.15 m)纵截面上温度场分布。可以看出，3种工况下客厅温度场分布表现出了较一致的特征：地面附近温度梯较大，沿地面外法线方向呈现出典型的温度分层现象。但温度高于294 K的范围，工况12的最大，工况11的最小。这是因为工况11的开窗面积最大，室内气流速度平均水平较高，冷热气流的热交换更充分，从而形成了较厚的冷气湖。

图7(Fig. 7)

图7 X=1.15 m截面温度场 Fig. 7 Temperature profile in X=1.15 m section

图 8为不同工况时客厅进风窗口中轴线(X=7.93 m)纵截面上温度场分布。可以看出，3种工况下，2#房间地面附近均存在较厚的冷气湖。到了2#房门时，由于过流断面的缩小，冷气湖厚度被大幅压缩。3种工况下，在2#房门和3#房门之间的顶部区域，形成了不同范围的热滞留区，这与该区域漩涡的生成有关。与图 7不同的是，该截面上温度高于294 K的范围，工况12的最大，工况11的最小。这是因为工况11中窗户4的开度要大于其他2个工况的开度，减小了室内气流运动的阻力，原来流经4#房间外通道并流向3#房间的部分气流此时会“短路”进入4#房间，进入3#房间的冷气流相对减少，温度整体上会偏高。

图8(Fig. 8)

图8 X=7.93 m截面温度场 Fig. 8 Temperature profile in X=7.93 m section

图 9为工况12的4个水平面上的温度分布。图 9(a)显示了冷气流自北朝南运动过程中对脚踝高度(Z=0.1 m)温度的影响特征，每一房间的不同位置温度差异显著。4#房间温度整体上偏高，这是因为未开窗增大了流入气流的运动阻力，而3#房间作为回流气流唯一通过的房间，温度偏低。图 9(b)为膝盖部位(Z=0.5 m)温度分布。与图 9(a)对照，发现各房间温度更趋均匀，温度水平有所提升，温度梯度有较大幅度减小。除了2#房间温度低于采暖最低要求温度外，其余房间绝大部分范围的温度均在采暖最低要求温度以上。图 9(c)显示了Z=1.1 m高度各房间水平面上的温度分布。各房间温度均匀性较图 9(b)有所改善，所有房间的温度均高于采暖最低要求温度值。另外，在靠近客厅南部位置出现了局部热点，这可能是因为在这一位置有旋涡存在，形成了热滞留区。图 9(d)显示了Z=1.7 m高度各房间水平面上的温度分布。各房间温度值及均匀性比较低的几个面上的值均有所提高，只有4#房间可能由于冷热气流的交汇对房间这一高度处的温度均匀性产生了一定影响，但对温度值的影响并不显著。

图9(Fig. 9)

图9 水平面上温度场(工况12) Fig. 9 Temperature profile in horizontal sections (case 12)

图 10给出了工况11~12沿房间高度不同水平面上的舒适性评价结果。可以看出，3种工况的舒适性表现出了相似的变化规律：PMV值沿房间高度方向逐渐增大，但均小于0。根据PMV-PPD的热舒适评价指标及ISO7730标准规定：以PPD≤10%为设计依据^[19]，即90%以上的人感到满意的热环境为热舒适环境，此时对应的PMV= -0.5～+0.5。但考虑中国的实际经济状况，一般计算时都将该值放宽至PPD≤20%^[20]，此时对应的PMV= -0.75～+0.75。所以，根据中国对舒适性要求的PMV值，对于工况10，在0.65 m以上的空间可以满足舒适性要求；对于工况11，在0.72 m以上的空间可以满足舒适性要求；对于工况12，在0.42 m以上的空间可以满足舒适性要求。这3个工况对应的高度均低于一般建筑室内确定PMV值的水平面高度1.1 m，但室内仍偏凉。另一方面，在1.1 m水平面上，工况10的PMV值接近-0.5，工况12的PMV值为-0.3，而工况11的PMV值近似为-0.6。所以，只有工况10和工况12能满足国际上对舒适性要求的PMV值要求，但房间也偏凉。为了克服由于偏凉而引起的不适感，居住者可以通过改变衣着来适应环境。实际上，冬季室内温度偏低，反而有利于抑制污染物的释放^[21]。

图10(Fig. 10)

图10 PMV平均值沿房间高度变化曲线 Fig. 10 Average value of PMV in vertical direction

综上所述，本文所得的3种通风模式均能满足舒适性要求，工况12的舒适度最好，工况10次之，工况11稍差一些。

以兰州某民用住宅为研究对象，采用Realizable k-ε 模型对冬季自然通风时室内热环境进行了数值分析，得到了如下主要结论：

1)寒冷地区无空调设施建筑冬季通过人为调节的自然通风是可行的，但并不是单纯的开关窗行为，具体方案应与户型、建筑布局、当地气象资料、大气环境状况以及人们的生活习惯相结合；

2)室内气流组织与维持室温、降低室内污染物浓度密切相关，合理选择通风时段、组织气流，对于提高通风效率关系重大；

3)所得的3种通风工况均能满足舒适性要求，但存在室内偏凉现象，这有利于抑制室内污染物的释放。