温度层结下建筑物群周围流场影响的数值模拟研究

郭栋鹏，王 冉，赵 鹏，李云鹏，姚仁太，刘 瑶

(1.太原科技大学，太原 030024；2.中国辐射防护研究院，太原 030006；3.山西省环境保护技术评估中心,太原 030024)

对于核设施厂址(如核电站等)，其存在一些阵列式的建筑物群(如图1)，由于其复杂下垫面与大气稳定度的影响，使其周围流场变得非常复杂，流线弯曲、速度剧烈地不连续及湍流分布不均匀。其中，大气稳定度对气载放射性物质扩散的影响尤为重要。因此，为更加真实模拟气载放射性物质扩散情况以及对环境的影响，需要研究不同温度层结下复杂建筑物群对流动与污染物扩散的影响。

图1 某核设施厂址建筑物群Fig.1 Buildings at a nuclear facility site

对于这类特殊工程厂址对环境的影响起初使用风洞模拟手段进行[1-4]，随着计算流体力学(CFD)的发展，CFD模拟技术已逐渐用于预测各种建筑物周围流动特性以及建筑物周围大气污染物的扩散规律[5-8]，但这些研究均是建立在中性层结的基础上进行的。关于温度层结的数值模拟研究比较有限。Zhang等[9]使用k-ε模型在中性和稳定层结条件下研究了标准体近场污染物的流动与扩散规律，但其使用的是均一来流的入口条件，未考虑来流风切变的影响。Santos等[10]使用k-ε模型，以立方体为研究对象，研究了不同温度层结对污染物扩散的影响，并与风洞实验进行比较。Ashrafi等[11]与Orkomi等[12]应用CFD技术研究了不同温度层结烟羽的抬升与扩散规律。在我国，部分学者应用CFD数值模拟技术在中性层结下对建筑物近场污染物流动与扩散规律进行研究[13-17]，而关于不同温度层结下建筑物群近场的流动模拟研究尚不多。

本文采用STAR-CD提供的RNGk-ε模型对不同大气稳定度条件下建筑物对附近流场的影响进行了模拟。首先应用CFD技术对Yassin[4]风洞实验进行模拟，并将模拟结果与风洞实验结果进行对比分析，检验其对温度层结的模拟能力，建立温度层结CFD数值模拟方法。并参照某核设施厂址，选取8×7建筑物群，对不同温度层结(不稳定、中性和稳定)条件下对规则建筑物矩阵中流动的影响规律进行初步研究，为评价核设施厂址复杂下垫面对流动影响及其近场气载放射性污染物扩散影响奠定基础。

1 数值模拟

基于Reynolds averaged Navier-Stokes (RANS)方法的三维CFD模型STAR-CD V3.26被设计用于平坦和复杂地形中大气流动和污染物扩散模拟。其使用修正标准k-ε湍流模型求解稳态流场，并通过有限体积法(FVM)离散流动和湍流的控制方程。本模拟使用Boussinesq近似简化计算，即假设计算域中气体为不可压缩流体，且密度定常。连续性和动量守恒方程如下：

(1)

(2)

式中，xi(x1=x,x2=y,x3=z)分别为纵向、横向和垂直方向；Ui为沿xi的平均速度分量，m/s；P为压力，Pa；μt为湍流涡黏性；μ为运动学黏性；ρ为密度，kg/m3；SB=gρβ(Tamb-T)，为浮力源项；T为温度，℃；Tamb为环境温度，℃；g为重力加速度，kg/s2；β=1/T，为热膨胀系数，1/℃。

湍流动能(k)以及耗散率(ε)传输方程如下：

(3)

(4)

(5)

(6)

式中，Pk为切应力相关的机械k的产生项，m2/s2；Pb为浮力相关的k的产生项，m2/s2；σk、σε、σh为湍流普朗特数；Cε1、Cε2、Cε3、Cμ为湍流模型经验常数。表1给出了计算中模型常数。

表1 湍流模型常数Tab.1 Turbulence model constants

2 模拟描述

2.1 几何描述

本文研究不同温度层结条件下规则建筑物矩阵中对流动的影响。模拟使用的三维矩形障碍物高度(H)为10 m，其大小为H×3H×H(长×宽×高)。研究布局是由56个障碍物组成的矩阵，其规模为：8个障碍物分布于流向上，7个障碍物分布于横向上。矩阵流向长度(Lx)和横向宽度(Ly)分别为36H和33H。建筑间隔在流向和横向上分别为4H和2H。坐标系原点(x/H=0和y/H=0)位于矩阵水平横截面(z/H=0)的中心[见图2(a)]。矩阵第一排建筑物迎风面位于x/H=-18，最后一排建筑物背风面位于x/H=18。建筑物矩阵几何示意图详见图2。

图2 规则建筑群几何布局Fig.2 Geometrical configuration of the regular buildings

2.2 计算域和网格

模拟中计算域在流向和横向长度均为60H，计算域高度为25H，底部中心位于坐标系原点。使用212×208×58个笛卡尔网格对计算域进行离散化。网格在垂直方向、流向和横向上分布是不均匀的。在矩阵中，水平方向上网格尺寸为0.2H，水平网格分辨率在矩阵上风向和下风向区域逐渐减小。在2H高度以下网格保持均匀最小的垂直间距，最小垂直网格尺寸为0.1H，垂直网格间距在2H高度以上逐渐增大。使用一套精细的网格280×268×62进行网格独立性分析，图3比较了矩阵中心处使用精细网格和粗糙网格得到的流向速度垂直分布。散点图中数据点大致分布在1∶1线上，网格独立性分析结果表明两套网格中流向速度差异是很小的。由于使用精细网格并不会造成研究结论的改变，所以本文选用粗糙网格进行模拟。

图3 网格无关性验证Fig.3 The grid independence verification

2.3 边界条件

(1)入流和出流边界

模拟中计算域侧面边界类型需要根据入流风向确定，即上游边界设为入流边界，下游边界设为出流边界。本研究考虑垂直于建筑物矩阵的入流风向(平行于x轴)并基于莫宁-奥布霍夫相似理论来描述表面层。入口处平均风速和温度垂直廓线由下列公式定义：

(7)

T(z)=θ(zr)+

(8)

(9)

θ*=-q0/u*

(10)

q0=Qh/(ρCp)

(11)

式中，u为速度，m/s；z为高度，m；zr为参考高度，m；z0为地面粗糙度，m；θ为位温，℃；u*为摩擦速度，m/s；L为莫宁-奥布霍夫长度，m；θ*为温度尺度；λadia=-0.009 766，绝热直减率；κ为冯卡门常数，0.41；Qh为地面感热通量，W/m2；Cp为定压比热，kJ/(kg·K)；Ψ,ψ分别为稳定度修正函数，在中性条件下Ψ和ψ均等于0。

稳定条件：

(12)

(13)

不稳定条件：

(14)

(15)

(16)

(2)对称与壁面边界

计算域顶部及两侧面(平行于x轴)应用对称边界条件，这些边界上所有变量正向梯度均等于0。底部边界定义为无滑移壁面条件并设置粗糙度z0等于0.02 m，该壁面上速度分量设为0并使用标准壁面函数计算其它变量。壁面动量由如下对数律法则描述：

(17)

式中，E为壁面粗糙度函数；u+=Up/u*，是无量纲速度，其中Up是平行于壁面的流体速度；y+=ρu*y/μ，是壁面至网格中心的无量纲距离；y为网格中心到壁面的距离。

不同模拟边界层中参考高度处(zr=H)风速(UH)和温度(TH)分别设为2 m/s和298 K。表2给出了不同类型模拟边界层的微气象参数。

表2 模拟边界层主要微气象参数Tab.2 Micro-meteorology parameters of the simulated boundary layers

3 模型有效性分析

3.1 风洞实验描述

Yassin[4]风洞实验主要研究不同温度层结建筑物对流场及其污染物扩散规律的影响，该风洞实验段长16.0 m，宽1.2 m，高1.0 m。试验以1∶300制作模型，模型长度(L)、宽度(W)、高度(H)均为100 mm的建筑物。为了与风洞实验结果比较，数值模拟的计算区域设为16.0 m×1.2 m×1.0 m(长×宽×高)，风洞实验模型见图4，风洞实验主要参数见表3，不同温度层结下温度层结下速度廓线、温度廓线见图5。

3.2 模拟设置

为了与风洞实验结果比较，数值模拟的计算区域与风洞实验相同，数值模拟的计算区域设为16.0 m×1.2 m×1.0 m(长×宽×高)，模型尺寸与Yassin[4]风洞实验模型相同(见图4)，计算域中网格总数为180万。网格结构采用具有良好拓扑结构的六面体网格，计算区域内最大网格尺寸为30 mm，建筑物表面网格尺寸为5 mm。计算域中上风向设为入流边界，下风边界则为出流边界。计算域顶部应用对称边界条件，底部使用无滑移壁面条件。入流速度、温度和湍流廓线的描述均与之前(2.3节)保持一致，不同温度层结下，拟合入流速度，温度廓线见图5，并依据表3给出的数据设置输入参数。

图4 建筑物及其测量位置图Fig.4 Building and survey location map

图5 不同温度层结下速度廓线、温度廓线Fig.5 Vertical distributions of mean velocity and Temperature in the simulated boundary layer under thermal stability

表3 不同稳定度条件下主要参数Tab.3 Main parameters under different stability conditions

3.3 流场有效性分析

本文分别在沿建筑物中心线x/H=0.75、1.0、1.5和2.0(如图4)四个不同位置进行流场特征的有效性分析，在大气环境模拟领域，通常采用归一化速度消除不同模拟风速引起的建筑物对流场结构的影响差异。归一化速度(u/UH)为局地纵向平均速度(u)与来流建筑物顶部纵向平均速度(UH)之比。下风向不同距离处建筑物对周围流场影响的归一化速度比较结果见图6。

图6 不同位置处不同模型归一化速度随高度的变化Fig.6 Vertical profiles of mean velocity component in the longitudinal direction

由图6可知，不同温度层结下本文CFD数值模拟结果与Yassin风洞实验结果吻合较好，变化趋势基本一致。在建筑物背风向近地面速度显著减小，特别是当大气处于稳定状态时。随着下风距离的增大，气流混合逐渐均匀，风速逐渐恢复到来流状态，但是大气处于稳定层结时，风速恢复相对较慢，不稳定层结风速恢复相对较快。

当大气处于稳定状态时，在建筑物高度以下(z/H<1.0)近地面层速度减小略高于中性与不稳定状态，特别是在x/H=1.5、2.0处，而在x/H=0.75、1.0处，由于受建筑物机械扰动的影响，大气的温度层结对流场影响不明显。

4 建筑物群对流场影响

4.1 模型的建立

通过标准体对模拟方法的验证，将本文建立的温度层结数值模拟技术应用于核设施厂址理想建筑物群流场结构的模拟。本文使用56个高度(H)为10 m的矩形建筑物组成的矩阵来模拟核设施厂址理想建筑物群，每个障碍物大小 (x×y×z)为：H×3H×H。研究区域由8×7个建筑物组成，纵向(x轴)上每排分布有8个，横向上每列分布有7个。建筑物矩阵纵向距离和横向距离分别为36H和33H，建筑物之间纵向和横向间隔分别为4H和2H。笛卡尔坐标系原点位于障碍物矩阵中心，即x/H=0和y/H=0两条线在z/H=0平面上的交汇处(见图2a)。第一排建筑物迎风面位于x/H=-18，最后一排建筑物背风面位于x/H=18。该建筑物矩阵的几何布局见图2。

4.2 结果分析

(1)纵向速度，u/UH

图7显示了归一化流向速度u/uref在建筑物矩阵中心面y/H=0上的分布。不同稳定度条件下，u/UH均出现了明显的垂直梯度。街谷内部出现了大范围的负向速度区域，且在街谷底部负向速度值较大，u/UH<-0.2的区域出现在街谷底部区域(z/H<0.5)。这说明不同稳定度条件下，街谷底部出现了强烈的回流。不同稳定度条件下，第一个街谷中u/UH<-0.2的区域明显大于后续街谷，第二排建筑物之后各个街谷中u/UH的分布大体上是相似的，且相对于第一个街谷，不同稳定度条件下u/UH<-0.2的区域均发生减小。这说明第一个街谷中涡旋的强度是非常强烈的，因为气流通过第一排建筑物时产生的强烈的机械扰动对涡旋起主导作用。第二个及后续的街谷中，u/UH<-0.2的区域在不稳定条件下最大，中性条件下次之，稳定条件下最小，该区域在稳定条件下减小最为明显。这是因为第二排建筑物后，温度层结对街谷涡旋的影响变得明显，稳定的大气条件明显减弱了涡旋强度。不同稳定度条件下街谷底部上游部分|u/UH|值大于下游部分。

图7 不同稳定度条件下，y/H=0平面上归一化速度u/UH的等值线云图Fig.7 Contours of u/UH on the plane y/H=0 under different stability conditions

图8给出了街谷4中不同位置处u/UH的垂直分布情况。不同位置在z/H<0.5范围内u/UH值均为负，这证明了街谷内1/2建筑高度以下范围内以回流为主。z/H<0.5范围内，水平位置x/H=-1过渡至1，|u/UH|值逐渐减小，表明回流速度从街谷内上游至下游部分发生减小。在x/H=-1和1处，不稳定条件下的|u/UH|值略大于中性和稳定条件的，这表明不稳定条件下回流速度较大。街谷中，不同稳定度条件下，u/UH随高度增加逐渐趋于入流状态，而稳定条件下u/UH恢复至入流水平的高度略大于中性和不稳定条件。这表明稳定条件下建筑矩阵对较高处风速的衰减影响较大。

图8 街谷4中x/H=-1、0和1处u/UH的垂直分布Fig.8 Vertical profiles of u/UH at x/H=-1,0 and 1 in the canyon 4

(2)垂直速度，w/UH

图9显示了归一化垂直速度w/UH在建筑物矩阵中心面y/H=0上的分布。在第一排建筑迎风角附近，强烈的风切变作用造成了w/UH强烈的垂直分布。街谷中，w/UH在中部和迎风侧附近为负值，在背风侧附近为正值。这表明气流在街谷背风侧和迎风侧附近的垂直运动分别是向上和向下的。街谷内部，背风侧向上的运动属于上游涡旋的一部分，而迎风侧向下的运动则属于下游涡旋。从图中可以看出，不同稳定度条件下，背风侧附近|w/UH|值要略大于迎风侧，说明了背风侧向上的气流运动要强于迎风侧向上的运动。这反映出街谷内部上游涡旋流动要强于下游涡旋。街谷背风侧附近，w/UH值在不稳定条件下要略大于中性和稳定条件，且w/UH>0.2的区域在不稳定条件下是最大的，在稳定条件下w/UH均小于0.15。这是因为稳定的大气条件抑制了垂直方向的气流运动，而在不稳定条件下，这种运动是更强烈的。

图9 不同稳定度条件下，y/H=0平面上归一化速度w/UH的等值线云图Fig.9 Contours of w/UH on the plane y/H=0 under different stability conditions

图10 不同稳定度条件下，y/H=0平面上归一化湍流动能的等值线云图Fig.10 Contours of on the plane y/H=0 under different stability conditions

图11 街谷4中x/H=-1、0和1处的垂直分布Fig.11 Vertical profiles of at x/H=-1，0 and 1 in the canyon 4

(4)流场结构

图12显示了前3个街谷中y/H=0平面上速度场流线图。不同稳定度条件下，所有街谷内部可以观察到两个涡旋的出现，这流动模式即OKE[18]提出的尾迹干扰流动。不同稳定度条件下，街谷1中的涡旋结构略不同于后续街谷(街谷2和3)，特别是在稳定条件下，这种差异更为明显。稳定条件下街谷1中上游涡旋规模较小且并不完整。在街谷2和3中流动趋于稳定，不同稳定度条件下街谷中上游涡旋的规模明显大于下游涡旋。相较于街谷1，街谷2和3中上游涡旋规模略微减小，上游涡旋中心距街谷背风侧的水平距离略也较短，而下游涡旋规模出现略微增大。不稳定条件下，上游涡旋的中心高度和尺度均大于中性和稳定条件。下游涡旋在稳定条件下也略小于不稳定和中性条件。

图12 前3个街谷中y/H=0垂直剖面的流场结构Fig.12 Streamlines of the velocity field on the plane y/H=0

5 结论

通过对大气处于不同温度层结时建筑物群对周围流场结构影响的数值模拟研究，本文首先应用Yassin[4]风洞实验结果对数值模拟有效性进行分析比较，其次研究不同温度层结时理想建筑物群对流场结构的影响。结果表明：

(1) CFD 模拟结果与Yassin[4]风洞实验结果能较好地吻合，并且研究发现稳定层结建筑尾流区范围内速度和湍流动能减小。

(2) 不同温度层结下，街谷内流向速度受建筑物影响均出现很大程度的减弱。在街谷底部均有强烈的回流出现，并且不稳定条件下街谷底部回流速度大于中性和稳定条件。在第2街谷之后，不同温度层结下街谷内u/UH分布开始达到平衡。不稳定条件下，建筑物对纵向通道内u/UH的影响较小，纵向速度分布在与入流状态基本一致。而中性和稳定条件下，通道内u/UH随距离的增大逐渐出现明显衰减，并且z/H<2.0范围内u/UH分布逐渐趋于线性。稳定条件下建筑物矩阵对纵向速度的最大影响高度大于不稳定和中性条件。

(3) 不同温度层结下街谷内上游和下游部分均有涡旋出现。并且温度层结对涡旋大小和中心位置有着显著的影响。不稳定条件下，街谷涡旋强度大于中性和稳定条件。稳定条件下由于大气层结的抑制，街谷涡旋强度是较弱的。

不同温度层结下，由于建筑物群内部与环境大气温度存在温度差，从而影响了流场结构与湍流特征，本文通过与风洞实验结果的对比，建立了温度层结数值模拟技术，为下一步应用于评价不同温度层结下核设施建筑物群近场气载放射性污染物的流动与扩散奠定了基础。