随着工业化和城市化进程逐渐加快，人口密集的城市地区的空气质量备受关注。影响城市地区空气质量的主要污染源有扬尘、机动车尾气、恶臭和其他有毒有害污染物^［1］，其中，某些建筑（如住宅、医院等）内产生的废气污染物，往往会通过屋顶烟囱直接排放到大气环境中。在不利气象条件下，烟囱排放污染物可能会通过“下洗效应”返回到近地面，或者通过建筑通风系统进入室内，从而对近地面室内外空气质量和人体健康造成危害^［2-3］。

近年来，国内外学者对城市建筑环境间流动及污染物扩散进行了较为广泛的研究，主要采用数值模拟和风洞试验等方法^［4-6］。随着计算能力的提升以及各种算法的不断改进，数值模拟方法以其成本低、效率高、直观性强等优点，已经成为研究流动、传热与传质过程的重要方法之一^［7］。研究表明，污染物扩散主要受气象条件（如风向、风速、太阳辐射等）、建筑几何形状（高宽比、屋顶结构、架空结构等）、烟囱设置（尺寸、位置、排放速率）等影响^［8-10］。多位学者研究了孤立建筑物在不同风向角下的污染物扩散，结果显示，方位角等的改变会显著影响污染物的浓度分布^［11-13］。但是相比于研究孤立建筑物，更具有现实意义的是研究在城市环境中相邻或周围建筑对屋顶烟囱排放污染物扩散的影响。

本文采用standard k−ε两方程模型，并结合三维稳态对流扩散方程，利用商业软件FLUENT对不同风向（风向角α = 0° ～ 180°）下，下台阶建筑屋顶烟囱释放污染物进行数值模拟，研究污染物的扩散分布特性，分析建筑周围流场分布，并着重讨论不同建筑壁面尤其是背风面上和建筑间街谷内行人呼吸高度处的污染物分布情况。

1 研究方法 1.1 物理模型

本文选取的建筑物理模型如图1所示。图中以1∶100的缩尺建立下台阶式两栋建筑模型，并以矮层建筑（记为“A”）的高度作为参考高度H（0.2 m），建筑A尺寸为2H × 1.5H × H，高层建筑（记为“B”）尺寸为H × 1.5H × 2H；烟囱污染源位于建筑A屋顶正中央，直径为0.05H，高度为0.08H。按照来风向与建筑中心轴线的夹角设定不同的风向角α。图2中以俯视视角展示了0°、30°、60°、90°、120°、150°和180°七种风向角情况。

1.2 数值模型 1.2.1 控制方程

本研究涉及的所有工况中都将气流运动视为不可压缩流动，且不考虑热效应。气流运动的控制方程为连续性方程和动量方程（RANS）。

不可压缩流动的连续性方程

$\qquad \dfrac{{\partial {u_j}}}{{\partial {x_j}}} = 0 $

(1)

动量方程

$\quad \begin{aligned} \dfrac{{\partial ({u_j}{u_i})}}{{\partial {x_j}}} = & \dfrac{\partial }{{\partial {x_j}}}\left( {(\mu + {\mu _{\rm{t}}})\dfrac{{\partial {u_i}}}{{\partial {x_j}}}} \right) - \dfrac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}}{\text{ + }}{{g}_{{i}}} \end{aligned} $

(2)

式中：x_i、x_j为笛卡尔坐标；u_i、u_j为x_i、x_j方向上的平均速度分量；μ为分子黏度；μ_t为湍流黏度；ρ为密度；p为压力；g_i为重力加速度的分量。

研究表明，对于预测城市湍流流动，standard k−ε模型^［14］的计算结果与实验数据相比^［15-17］，具有较好的数值精确度。因此，本文选择standard k−ε模型，具体方程为

k方程

$\qquad \dfrac{{\partial \left( {\rho k{u_{_j}}} \right)}}{{\partial {x_j}}} = \dfrac{\partial }{{\partial {x_j}}}\left( {\dfrac{{{\mu _t}}}{{{{\text{α } _k}}}}\dfrac{{\partial k}}{{\partial {x_j}}}} \right) + {G_k} - \rho \varepsilon $

(3)

ε方程

$ \begin{array}{c}\dfrac{\partial \left(\rho \varepsilon {u}_{j}\right)}{\partial {x}_{j}} = \dfrac{\partial }{\partial {x}_{j}}\left(\dfrac{{\mu }_{{\rm{t}}}}{{\text{α}}_{\varepsilon }}\dfrac{\partial \varepsilon }{\partial {x}_{j}}\right) + \dfrac{\rho \varepsilon }{k}\left( {{{C}}_{1\epsilon }{G}_{k} - {C}_{\text{2}\varepsilon }\varepsilon } \right)\end{array} $

(4)

$\qquad {\mu _{\rm{t}}} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon } $

(5)

式中：k为湍动能；ε为湍动能耗散率；C_1ɛ、C_2ɛ和C_μ等均采用FLUENT软件中的默认参数，值分别为1.44、1.92和0.09；α_k和α_ɛ分别为与k和ε对应的湍流普朗特数，分别取为1.0、1.3；G_k为由平均速度梯度引起的湍动能。

$\qquad {G_k}= -{\rho {\overline {{{u'}_ i}{{u'}_ j}}}}\dfrac{{\partial {u_i}}}{{\partial {x_j}}} $

(6)

式中，$u_i^{'}$、$u_j^{'} $分别为x_i、x_j方向上的脉动速度。

采用Boussinesq假设的方法评估G_k，有

$\qquad {G_k} = {\mu _{\rm{t}}}{S^2} $

(7)

式中，S为平均应变率张量的模量。

$\qquad S = \sqrt {2{S_ {ij}}{S_ {ij}}} $

(8)

式中，S_ij为应变率张量。

对于烟囱释放污染物过程采用三维稳态对流扩散方程^［18］，即

$\qquad {u_i}\dfrac{{\partial {C^\beta }}}{{\partial {x_i}}} = \dfrac{\partial }{{\partial {x_j}}}\left( {\left( {D_{\rm{m}}^\beta + \dfrac{{{\nu _{\rm{t}}}}}{{{{S c_ {\rm{t}}}}}}} \right)\dfrac{{\partial {C^\beta }}}{{\partial {x_j}}}} \right) + {S_ {\rm{p}}} $

(9)

式中：C^β为第β种组分的质量分数；$ {D}_{\text{m}}^{\beta } $为分子扩散系数；v_t为湍流涡黏性；Sc_t为湍流施密特数，值为0.7；S_p为污染源项。

1.2.2 计算域

采用计算流体力学（CFD）软件ANSYS FLUENT 14.5进行数值模拟，并按照前处理—分析—后处理的流程对所研究问题进行求解。图3以风向角α =0°为例，展示了数值模拟时的计算域设置，其中，x – y – z直角坐标系的原点(x = 0, y = 0, z = 0)位于建筑A、B中心处，计算域大小的选取原则严格遵循AIJ^［19］的指导指南，入口平面在建筑A迎风墙上游6H处，出口平面在建筑B背风墙下游15H处，顶部平面距离建筑B屋顶6H，侧面平面距离建筑侧壁6H。

1.2.3 求解计算

入口采用速度入口，其风速U(z)、湍动能k以及湍动能耗散率ɛ分别定义为^［20］

$\qquad U(z) = {U_{{\text{ref}}}}{\left( {\dfrac{{\textit{z}}}{H}} \right)^a} $

(10)

$\qquad k(z) = \dfrac{3}{2}{\left[ {I({\textit{z}})U({\textit{z}})} \right]^2} $

(11)

$\qquad \varepsilon ({\textit{z}}) = 0.187{{{C}}_\mu }^{\tfrac{1}{2}}k({\textit{z}})\dfrac{{{U_{{\text{ref}}}}}}{H}{\left( {\dfrac{{\textit{z}}}{H}} \right)^{a - 1}} $

(12)

$\qquad I({\textit{z}}) = 0.1{\left( {\dfrac{{\textit{z}}}{H}} \right)^{ - a - 0.05}} $

(13)

式中：U_ref为选取的参考风速，值为3 m·s⁻¹；a为风剖指数，值为0.18；H为0.2 m；I(z)为湍流强度。

出口设置为流出边界条件（outflow），侧面以及顶面设置为对称边界条件（symmetry），建筑壁面以及地面设置为无滑移边界条件且近壁面区域采用标准壁面函数（standard wall functions, SWFs）。控制方程采用有限体积法（FVM）离散，压力与速度采用压力耦合方程组的半隐式方法（SIMPLE）进行耦合，对流离散项采用二阶迎风格式，各变量相对残差小于10⁻⁶，并且流场无变化即视为收敛。

无量纲浓度K定义为

$\qquad K = \dfrac{{{C_{\rm{r}}}{U_{{\text{ref}}}}{H^2}}}{{{C_{\rm{e}}}Q}} = 10\;695.2{C_{\rm{r}}} $

(14)

式中：C_r为各监测点或平面处污染物摩尔分数，本文选取六氟化硫（SF₆）作为污染物；C_e为烟囱排放口处污染物摩尔分数，值为1/7；Q为烟囱出口处污染物排放率，值为7.85 × 10⁻⁵ m³·s⁻¹。

1.2.4 网格独立性分析

数值仿真计算时既要做到不失精度，又不浪费计算资源，因此选取合适的节点数网格非常重要。本文在烟囱释放污染物工况下α =0 °时设计了3种不同节点数网格，分别为网格1（91万）、网格2（172万）和网格3（327万），且计算域完全对称，为了节约计算时间，所以仅划分了一半的网格。提取街道峡谷中心轴线处无量纲浓度进行对比分析。网格独立性分析如图4所示。经比较，网格1和网格2的无量纲浓度相差较大，最大差值大于10 %，而网格2和网格3的无量纲浓度相差较小，最大差值小于5%。说明在网格2的基础上再增加节点数对计算结果几乎无影响，因此本文选取网格2进行数值仿真计算。图5为建筑周围以及烟囱附近网格，其建筑壁面最近一层网格高度为0.025H。

2 结果与分析 2.1 流场以及污染物浓度分布随风向角α的变化

图6为不同风向角α下流场分布以及行人呼吸高度处无量纲速度u/U_ref分布，图7为不同风向角α下烟囱释放的污染物浓度分布，并且为了便于观察，将建筑物做半透明处理。经对比分析可知，如图6（a）所示，当α = 0°时，经过建筑A屋顶的来流，由于受建筑B的阻挡，大部分气流流向街道峡谷，并在街道峡谷内形成逆时针旋涡；街道峡谷内行人呼吸高度处无量纲速度为负值，表明街道峡谷存在旋涡。在这种气流运动影响下，如图7（a）所示，烟囱释放的污染物大部分跟随气流进入街道峡谷内。此外，由于建筑A两侧也存在气流产生的漩涡，污染物被气流从街道峡谷内带动至建筑A的两侧壁面处，而建筑B四周壁面则被经过烟囱的来流包裹，污染物也流动至建筑B的四周壁面处。如图6（b）～（f）所示，在α = 30°、60°、90°、120° 、150°五种工况下，当来流方向与建筑中心轴线成一定角度后，经过烟囱的大部分气流未进入街道峡谷内，因此污染物扩散对街道峡谷内行人呼吸高度处平面影响相对较小。如图7（e）、（f）所示，当α = 120°、150°时，污染物均未扩散至街道峡谷内，但由于气流在建筑背风区后形成旋涡，并且在行人呼吸高度处平面无量纲速度也均为负值，说明气流产生了漩涡，且对污染物扩散至背风壁面处有一定的影响。如图7（b）、（c）所示，当α = 30°、60°时，污染物扩散至建筑B的背风壁面处；如图6（d ）～ （f）所示，当α = 90°、120°和150°时，经过烟囱的气流会在建筑A背风区形成漩涡。如图7（d）～（f）所示，污染物会扩散至建筑A的背风面，但对建筑B周围壁面几乎无影响；如图6（g）所示，当α = 180°时，两建筑街道峡谷下方有明显的逆时针漩涡，在其上方则形成了一个顺时针的漩涡，这与黄晓天等^［21］的研究结果一致，如图7（g）所示，烟囱释放的污染物扩散至建筑B背风壁面以及建筑之间的街道峡谷内。

2.2 烟囱污染物浓度分布特性

图8为建筑四周壁面以及行人呼吸高度处平面对应编号。图9为不同风向角α下建筑四周壁面以及行人呼吸高度处平面上平均无量纲浓度$ {\overline K} $［${\overline K} = {{{\overline C}_{\rm{r}}}{U_{{\rm{ref}}}}{H^2}/\left( {{C_{\rm{e}}}Q} \right)}$，其中${\bar C_{\rm{r}}} $表示壁面处污染物的平均摩尔分数］比较。当α = 0°时，壁面A−4与行人呼吸高度处平面P上平均无量纲浓度$ \overline {K} $均处于一个相对较高的水平。这是由于气流直接经过烟囱后，会携带污染物进入建筑间的街道峡谷内，但壁面B−1上的平均无量纲浓度较小。这是因为形成的逆时针漩涡致使污染物大部分聚集至壁面A−4以及行人呼吸高度处平面P上。建筑侧壁面也会受到污染物扩散的影响，结合图6、图7可知，壁面A−3主要是气流在建筑A下游区域形成回流区导致的。当α = 30°、60°时，平面P上基本无污染物，壁面A−3上有少许的污染物，且2种风向角下污染物在壁面上的平均无量纲浓度很接近，当α = 30°时相比于α = 60°时，烟囱释放的污染物会更多地影响壁面B−3和B−4；当α = 90°时，仅建筑A的背风面A−3上平均无量纲浓度较高。这是因为建筑B在该风向角下基本不受经过建筑A屋顶烟囱的气流影响，此时流动结构相对简单，气流在建筑A背风区形成的漩涡致使污染物少量聚集于壁面A−3附近。当α = 120°、150°时，烟囱释放的污染物仅对建筑A的背风面有影响；当α = 180°时，流动结构变得相对复杂，建筑四周壁面以及平面P均有污染物扩散分布，且由于气流在建筑B背风区产生漩涡，将污染物卷吸至建筑间的街道峡谷内，且壁面B−1的无量纲污染物相对较大。随着风向角逐渐变大，建筑A屋顶中心烟囱释放的污染物对壁面A−1的影响逐渐变大，对壁面A−3的影响则是先变大后变小，而只有α = 0°和180°时污染物才会扩散至壁面A−2和A−4，尤其是当α = 0°时处于街道峡谷内的壁面A−4上污染物的平均无量纲浓度最大。建筑B四周壁面以及平面P只有当α = 0°、180°时才有污染物到达；不同的是当α = 180°时建筑B背风面的尾流漩涡将部分污染物卷吸至平面P处，该平面平均无量纲浓度比α = 0°时的减少80%。

3 结　论

本文采用CFD数值模拟方法，研究不同风向角下两建筑模型中矮层建筑屋顶烟囱排放污染物的扩散分布特性，得到以下结论：

（1）气流来流处是否有建筑对流动结构影响较大，来流方向无建筑时，仅会在烟囱所在建筑背风面形成旋涡，污染物也只会扩散至建筑物背风区域以及下游区域；当α = 180°时，来流处存在建筑且高度相对较高时，气流会产生回流区，将污染物带回街道峡谷区域。

（2）来流风向角的变化对烟囱释放污染物扩散影响较大，会对建筑不同的壁面产生显著影响。壁面A−1、A−2、B−1在α = 180°时污染物的平均无量纲浓度达到最大，壁面A−3在α = 90°时污染物平均无量纲浓度均达最大，壁面A−4、B−2、B−3、B−4以及街道峡谷内行人呼吸高度处平面P在α = 0°时污染物的平均无量纲浓度最大；壁面A−1在α = 0° ～ 60°时基本不受污染物影响，壁面A−2、A−4、B−1、B−2以及平面P在α = 30° ～ 150°时基本不受污染物影响，壁面A−3在α = 150°时污染物的平均无量纲浓度最小，壁面B−3、B−4在α = 90°、120°、150°时基本不受污染物影响。