不同温度层结城市街区内气流与污染物扩散的数值模拟研究

史学峰，郭栋鹏，王 冉，姚仁太

(1. 中国辐射防护研究院，太原 030002；2. 太原科技大学，太原 030024)

前 言

城市是人类工作和生活的聚居地，又是一个污染源密集的地区，城市区域污染问题日益突出，特别是有毒有害气体的释放对人体健康造成威胁，使得城市区域复杂下垫面的流场和污染物扩散问题越来越受到关注。目前，大气污染物扩散的风洞试验及数值模拟多数建立在中性层结下[1～5]，但是实际大气并非如此，近地面层的大气经常在稳定层结、不稳定层结以及中性层结之间变化，与中性层结下相比，温度层结下研究街道峡谷污染物流动与迁移的数值模拟研究较少[6～9]，城市地区，温度层结主要是由白天在建筑表面或地面上直接产生的太阳辐射加热附近的空气引起的，从而影响城市街区内污染物的流动与扩散。因此，为更加真实模拟城市街区内大气污染物扩散情况以及对环境的影响，需要研究在不同温度层结下的污染物的扩散。

关于城市街区中温度层结对污染物流动与扩散的风洞实验研究相对有限，主要是由于风洞在控制气流或表面温度存在一定技术难度[10]。Uehara K等人[11]利用层结风洞技术研究了不同大气稳定度(Rib=-0.21-0.79)对街道峡谷中流场影响，结果表明：当大气处于稳定条件下时，处于街区中的空腔涡旋趋于减弱，反之增强，大气处于稳定时，由于浮力作用导致街区内向下气流减弱，当Rib大于0.4时，街区内速度接近于0。Rafailidis[12]在英国萨里大学的EnFlo风洞中研究了温度层结对城市街道峡谷通风的影响。测量结果表明，稳定的分层条件导致污染物会滞留在峡谷内。Kovar-Panskus等人[13]使用风洞研究了太阳辐射诱导的壁面加热对2-D街道峡谷内流动的影响。通过比较有无壁面加热的情况发现，在极低的弗洛德数下，加热迎风面壁面，在靠近峡谷地面的位置产生了非常弱的二次流，但是没有说明壁面加热对污染物的流动与扩散的影响。Yassin等人[14]通过风洞实验和现场试验分别研究了不同稳定度条件下不同风向城市街道峡谷中污染物的扩散，结果表明稳定度和风向是影响大气扩散的影响因素，稳定层结下最大浓度高于中性和不稳定层结，该风洞和现场数据可以用于计算流体动力学(CFD)验证。

随着CFD 的不断发展, CFD模拟技术已逐渐用于估算各种建筑物环境下大气污染物的扩散规律，Tominaga等人[15]出版了基于稳态RANS模拟的“ AIJ建筑物周围行人风环境CFD实际应用的指南”。Blocken等人[16]将十步法应用于环境流体力学模拟，并通过两个回顾性案例分析表明，该方法适用于环境流体力学方面的研究。除了这些一般准则之外，还发布了一些非常具体的准则。其中包括(1)在计算域中对大气平衡边界层进行一致建模，(2)生成高质量的网格，(3)使用现场和实验室数据进行验证。前两个领域的大部分工作都集中在稳态RANS上，这些准则的建立为更准确和可靠的CFD模拟提供支持。Sini and Anquetin[17]使用CFD研究了空腔区内的机械浮力诱导流动，结果表明表面温度变化会影响空气流动和污染物交换，甚至不同的温差可以引起街区内湍流结构的转化，如由单涡向多个反向旋转涡流转换。Kim和Baik[8]使用2-D k-e模型研究了不同纵横比和街区底部加热对街区内流动的影响。结果模拟出不同形式的流动，表明底部加热对街区内流动起重要作用，但该研究并没有进行污染物扩散的研究。Xie等人[6,18-19]研究了在独立街道峡谷和不同横纵比的街区中由机械浮力诱导的流动，进一步表明了浮力在街道峡谷流动和污染物输送问题中的重要性。研究表明，温度层结下街道峡谷中的风特性是无量纲参数Gr/Re2的函数。Tan等人[7]用CFD模拟中性与不稳定层结下，城市街区对流场与污染物扩散的影响，结果表明，不稳定层结下，地表温度显著影响流场结构与污染物的扩散。随着理查森数的增大，近地面污染物的浓度逐渐减小，但其在近地面模拟误差较大，这可能由于使用二维模型所导致的。在这里需要说明的是，风洞实验中的街区是由三维方形建筑物排列组成，但在实际验证时大部分学者采用了简化的二维模型进行验证[6～8]，结果吻合存在一定的误差，一方面可能与用二维简化模型模拟三维模型有关[6]，另一方面可能是由于使用不同湍流模型的原因或者与使用不同的壁面有关[20]。

在我国，大多研究人员仅研究中性层结下复杂建筑物群周围的流场及其污染物的扩散规律进行了相关的风洞实验研究[21～24]，对不同温度层结下复杂街区对流动与污染物扩散的影响的研究尚不多见[25～26]。

本文主要应用STAR-CD提供的RNGk-ε模型对不同温度层结下地面污染源在城市街区内中间位置时，城市街区对附近流场以及污染物扩散的影响，建立的3-D CFD物理模型与Uehara K[11]风洞试验模型相同，采用布尼斯克(Bonssinesq)假设在动量方程添加了外部源相模拟温度层结，通过莫宁奥布霍夫长度(L)描述不同温度层结下大气的运动，并应用Pontiggia M.等人[21]方法对RNGk-ε模型k方程与ε方程进行修正，探索受复杂建筑物群影响不同温度层结CFD模拟技术，预期建立数值风洞。

1 风洞试验

Uehara K等人[11]试验是在日本国家环境研究所(Japanese National Institute for Environmental Studies)风洞试验室进行，该风洞试验段长24m，宽3m，高2m，实验段风速范围为0.2～10 m/s，应用来流空气温度(Ta)与风洞底板温度(Tf)控制试验过程中温度廓线，城市街区模型缩比为1∶300，所有模型长度、宽度、高度(H)均为100mm的立方体，风向轴线(X方向)上建筑物间隔为100mm，因此街谷的比例为1∶1(宽度∶高度)，横上(Y方向)建筑物间隔为50mm，风洞试验街区模型共有13列，主要研究第5与第6列建筑物形成城市街区内对流场影响,见图1。风洞试验共进行了7种温度层结的研究，整体理查森数(Rib)变化范围为-0.21-0.79，不同温度层结下主要参数见表1。为了进一步研究城市街区对污染物扩散的影响，本文在第5与第6列建筑物形成城市街区地面中心位置处设置地面源，地面源高度(Hs)为2mm(0.02H)内径为4mm。

2 数值模拟

2.1 计算域与网格设置

表1 Uehara K等人风洞试验主要参数Tab.1 Main experimental parameters in Ueharaet al

图1 模拟几何示意及建筑物周围网格划分图Fig.1 Simulated geometry and the computational domain mesh

CFD模拟过程，平衡大气边界层方法至关重要，因此，CFD数值模拟应首先对所构造的大气边界层的水平均匀性进行检验。为了检验边界层的保持性，本文先进行了空洞模拟，空洞计算域尺寸与本文数值模拟计算区域相同，检验结果如图2所示，由图2可知，不同距离处CFD数值模拟大气边界层较好的保持。

图2 不同位置处速度变化(Ld为计算域的长度) Fig.2 Vertical profiles of the u velocity (Ldis the length of the domain).

2.2 湍流模型

大气环境中用来描写流体运动的基本方程是Navier-Stoke方程，空气流动速度一般不大，可将空气当作不可压缩流体。近地层气流的运动一般为复杂的湍流运动，目前常采用雷诺时均法对基本控制方程进行处理，空气温度、密度变化不大，流动符合Bonssinesq假设，在动量方程添加了外部源相(SB=ρβ(Tamb-T)g, 其中ρ 为空气的密度, Tamb为环境温度, β= 1/T 为膨胀系数)。由于RNG k-ε湍流模型能够较好的模拟绕流问题[15]，本文选用RNG k-ε模型封闭N-S方程进行计算。为了在CFD模拟中实现不同温度层结，本文通过添加外部源相对选用的湍流模型进行修正，该外部源相添加方法由Pontiggia M.等人提出，并将中性与稳定条件下CFD模拟结果与现场试验结果进行对比分析[27]。guo等人应用Pontiggia M.等人提出方法对温度层结下建筑物对流场与污染物扩散的影响进行了研究[28]。

2.3 边界条件

本文CFD数值模拟入口边界风廓线与温度的取值及其变化规律均与Uehara K等人应用风洞试验研究不同温度层结下城市街区对周围流场结构影响的模型相同，计算时入口边界条件风廓线见图3。污染源排放速度为0.08m/s，通过调整入口边界大气温度(Ta)以及计算区域底部温度(Tf)实现不同稳定度的模拟。不同高度处湍流动能廓线与湍流耗散廓线分别根据式(1)，(2)计算。

图3 入口边界条件风廓线Fig.3 The vertical velocity profiles of the approaching flow.

(1)

(2)

出口边界条件按自由出口设定。计算域两侧及顶部采用滑移壁面条件。建筑物表面和地面设置一定的摩擦速度(u*)与粗糙度(z*)，采用无滑移的壁面条件并由壁面函数法确定壁面附近流动，使用标准壁面函数处理近壁面的流动。示踪气体的扩散采用被动扩散，示踪气体温度与环境温度相同。不同稳定度条件下主要参数见表2。

采用压力-速度修正算法 SIMPLE[30]联立求解各离散方程，采用精度较高QUICK差分格式。

表2 不同稳定度条件下主要参数Tab.2 Main parameters for different stratified conditions

3 结果分析与比较

本文主要研究不同温度层结下，地面源位于距城市模型前边界1m第五与第六排建筑物街区内中间位置，第五与第六排建筑物街区对流场、温度场与污染物扩散的影响。主要分析沿中心线(y=0)街区内, 街区中间位置(x/H= 0.0)流场结构、温度场以及街区背风侧(x/H= - 0.5)、迎风侧(x/H= 0.5)与街区内浓度场分布特征。

3.1 流场结构

研究采用RNG k-ε湍流模型计算了城市街区对周围流场的影响，并采用无量纲的归一化速度来消除风速不同引起的差异。归一化速度(u/U700)或(w/U700)为局地纵向平均速度(u)或垂直方向上平均速度(w)与边界层顶部δ=700 mm处纵向平均速度(U700)之比。

3.1.1 纵向平均速度分布，u/U700

不稳定、中性与稳定层第五与第六排建筑物街区中心 x/H= 0.0 位置纵向上归一化速度数值模拟与Uehara K等人[11]风洞试验结果比较见图4。

图4 x/H= 0.0 位置纵向上归一化速度变化，u/U700Fig.4 Vertical profiles of mean velocity component in the longitudinal direction, u/U700

由图4可知，在街区内部与上部(0.1

3.1.2 街区内部流场结构

由图5可知，不同温度层结下，街区内部均形成一个稳定的涡，街区内部涡的强度受温度层结影响较明显，稳定层结下，街区内涡流较弱，不稳定层结下，街区内涡流较强，中性层结位于中间。

不同温度层结条件下街区内垂直方向上气流w/U700分布见图6。由图6 可知，在街区背风面出现了向上气流，街区迎风面出现了较强的向下气流，主要由于气流遇到下游建筑物迎风面时进入街区内部，到达地面，然后沿上游建筑物背风面上升，在街区顶部作为来流的补充[11]。稳定层结下向上与向下气流相对较弱，主要由于稳定层结抑制了垂直方向上气流的运动。稳定层结下迎风面最强w/U700= -0.14，背风面最强w/U700= 0.04；不稳定层结下向上与向下气流相对较强，迎风面最强w/U700= -0.28，背风面最强w/U700= 0.14；中性层结下，迎风面最强w/U700= -0.22，背风面最强w/U700= 0.1。

图6 街区内垂直方向上气流w/U700分布Fig.6 Contour ofw/U700 within the street canyon

3.2 湍流结构

图7 湍流动能随高度(z/H)的变化Fig.7 Contours of the y/H= 0 plane in the street canyon

3.3 温度场结构

稳定与不稳定层结条件下街区内中心轴线温度与归一化温度分布见图8。由图8可知，不稳定条件下,数值模拟结果比Uehara K等人风洞试验结果略高，稳定条件下,数值模拟结果比Uehara K等人风洞试验结果略低，但是温度梯度变化规律相似，均反映了街区内温度随高度的变化规律主要由于Uehara K等人风洞试验没有给出来流温度廓线，仅给出来流空气温度与地板温度，这增加了数值模拟的不确定性；近地面，数值模拟结果与Uehara K等人风洞试验结果的差异，主要由于数值模拟壁面函数对近地面热转化不准确，从而改变了近地面附近湍流混合[20]。不稳定层结与稳定层结在0.8

图8 街区内中心温度变化Fig.8 Profiles at the central part of the street canyon

3.4 浓度场分析

图9显示了在不稳定、中性和稳定层结条件下，街区背风侧墙面(x/H= - 0.5)、迎风侧墙面(x/H= 0.5)浓度场分布特征，并与Chang 和 Meroney[31]风洞试验进行对比分析，Chang和Meroney风洞试验主要研究中性层结下，地面源位于街区中间位置，不同城市街区高宽比以及污染源上风向不同列建筑物对污染物扩散的影响，本文污染源上风向有5列建筑物，扩散规律与Chang和Meroney风洞试验污染源上风向有8列建筑物的试验基本相同，均属于相同城市粗造度，因此本文街区对污染物扩散规律的影响主要与Chang和Meroney风洞试验 H/W=1并且污染物上风向有8列建筑物结果进行对比分析。

图9 背风面与迎风面浓度变化，K Fig.9 Concentrations on the upwind and downwind walls K

由图9可知，中性层结下，数值模拟结果与Chang和Meroney风洞试验结果吻合较好，风洞试验由于测量点位置的影响，背风侧墙上(x/H=- 0.5)最大值出现在z/H=0.2处，归一化浓度最大值 K=26.2；数值模拟最大值出现在z/H=0.05处，归一化浓度最大值 K=36.6，在z/H=0.2处，归一化浓度最大值 K=28.3。迎风侧墙(x/H=0.5)风洞试验最大值出现在z/H=0.45处，归一化浓度最大值 K=2.18，数值模拟最大值出现在z/H=0.47处，归一化浓度最大值 K=2.12。

由图9可以观测到热稳定分层对污染物扩散的影响，街区内部不同稳定度变化产生的浓度波动。不同温度层结，背风侧墙上(x/H=- 0.5)与迎风侧墙(x/H=0.5)数值模拟结果均表现为近地面稳定层结归一化浓度值最大，而不稳定层结归一化浓度值最小，中性层结位于中间，归一化浓度最大值均出现在近地面，主要由于稳定层结抑制了垂直方向上气流的运动。背风侧墙上(x/H=- 0.5)稳定条件下归一化浓度最大值 K=49.7，不稳定条件下归一化浓度最大值 K=28.2，中性条件下归一化浓度最大值 K=36.6；迎风侧墙上(x/H=0.5)稳定条件下归一化浓度最大值 K=2.75，不稳定条件下归一化浓度最大值 K=1.55，中性条件下归一化浓度最大值 K=2.12。不同温度层结下背风侧墙上(x/H=- 0.5)污染物的浓度均大于与迎风侧墙(x/H=0.5)浓度。主要由于在街区内部形成一个循环涡，顺时针气流携带污染物撞击到背风面墙上，因此背风侧墙上污染物归一化浓度值大于迎风面[32]。

图10为不同稳定度条件下沿风向轴线(y=0处)城市街区内污染物扩散的垂直剖面图，由图10可知，不同温度层结下，较高归一化污染物浓度均出现在污染源的上风向，稳定条件下归一化浓度值最大 K=412.6，不稳定条件下归一化浓度值最大 K=296.3，中性条件下归一化浓度值最大 K=361.6。主要由于在街区内部形成一个循环涡，顺时针气流携带污染物向上风向运动，并撞击到背风面墙上，从而在污染源的上风向形成高浓度区。稳定条件下,由于气流垂直运动相对较弱，街区内涡流的强度较小，污染物主要表现为在街区背风面归一化浓度值最大，不稳定条件下,由于气流垂直运动相对较强，街区内涡流的强度较大，污染物主要表现为街区背风面归一化浓度值相对较小。而中性层结位于中间。总之,街区内气流的交换主要受空腔区与来流的湍流混合及其混合强度的影响[26]。

图10 不同稳定度条件下沿风向轴线污染物扩散的垂直剖面图Fig.10 Vertical profile of pollutant diffusion along wind direction axis under different stability conditions

4 结 论

本文采用3-D计算流体力学(CFD)模拟不同温度层结下城市街区内气流与污染物扩散规律，并与Uehara K等人风洞试验结果进行对比分析，通过莫宁奥布霍夫长度描述不同温度层结下大气的运动，并应用Pontiggia M.等人方法对RNG k-ε模型k方程与ε方程进行修正，建立数值风洞模拟计术，结果表明。

4.1 Pontiggia M.等人对RNG k-ε模型k方程与ε方程进行修正方法可以应用在数值风洞中模拟大气边界层气体的流动与污染物扩散规律。

4.2 稳定层结下，街区内涡流较弱，不稳定层结下，街区内涡流较强；街区背风面出现了向上气流，街区迎风面出现了较强的向下气流，不稳定条件下背风面上升气流较强，稳定条件下背风面上升气流较弱。

4.4 不同温度层结下，数值模拟结果与Uehara K等人风洞试验结果温度梯度变化规律相似，均反映了街区内温度随高度的变化规律，近地面，数值模拟结果与Uehara K等人风洞试验结果的差异，主要由于地表粗糙度估算不准确，从而改变了近地面附近湍流混合。

4.5 由于热稳定分层对污染物扩散的影响，不同温度层结街区内部变化产生的浓度波动。不同温度层结下，较高归一化污染物浓度均出现在污染源的上风向，主要由于在街区内部形成一个循环涡，顺时针气流携带污染物向上风向运动，并撞击到背风面墙上，从而在污染源的上风向形成高浓度区，因此不同温度层结下背风侧墙上(x/H=- 0.5)污染物的浓度均大于与迎风侧墙(x/H=0.5)浓度，归一化浓度最大值均出现在近地面。稳定层结归一化浓度值最大，主要由于稳定层结抑制了垂直方向上气流的运动。

分析表明: 不同温度层结下，由于街区内部与环境大气温度存在温度差，影响了街区内流场结构与湍流特征，从而影响了街区内污染物的扩散特征。尽管本文提出的新模拟方法能够提高预测精度，但其只是阐明街区中温度层结对污染物流动与扩散影响机理的垫脚石。未来仍需对本研究中的方法进行严格的测试和改进，特别是在强稳定层结条件下流场结构及其污染物扩散规律的研究是下一步研究的重点。