随着人类生产活动的迅速发展，大气中气溶胶的含量逐年增长，长期悬浮在空气中的气溶胶粒子严重危害人类健康[1-2]，对气候产生负面影响[3]。悬浮在大气中的气溶胶是通过粒子随气流输运时发生自然沉降被去除的。因此气溶胶粒子扩散运动的研究一直以来都受到了环境科学和大气科学领域研究者的广泛关注[4]。在气溶胶粒子输运研究中，粒子扩散速度是一个重要参数，它可以更好地量化这种化学形态微粒物质和表面大气能量、动量交换的行为，也是连接粒子输运和沉积的关键因素[5-7]。
风洞试验是模拟气溶胶粒子扩散实验常用实验手段，它不仅可以克服现场实验条件的差异性、不确定性等因素，也可以为数据模拟实验的仿真性结果提供实践支撑[8]。
PING et al[9]使用粒子图像测速法(PIV)研究了风洞中沙粒气溶胶扩散速度的变化规律。研究表明，粒子平均纵向速度随风速的增大而增大，随粒径的增大而减小。BO et al[10]用PIV在风洞中研究了沙层上方沙粒粒径对粒子速度概率分布的影响，发现粒径大小对其速度分布并无显著影响。LI et al[11]采用粒子动态分析技术(PDA)研究了风驱动下气溶胶粒子的垂向速度随高度的变化，发现平均垂向速度随扩散高度变化主要受粒径影响，大粒径粒子的速度随扩散高度的增加而增加，而小粒子则呈现先减后增的趋势。JIANG et al[12]在风洞中用粒子跟踪测速法(PTV)测定了气溶胶粒子平均速度的变化规律，发现在大多数情况下，粒子平均垂向速度远小于平均纵向速度。
气溶胶迁移运动发生在对流层中最靠近地面的气层——大气边界层[13]。根据大气边界层垂直方向的高度，可将其分为冠层(从地面到平均障碍物高度的2倍)、粗糙子层(障碍物平均高度的2～5倍)和惯性子层(高于粗糙子层)[14-15]。到目前为止，文献中所报道的气溶胶粒子的扩散风洞研究大多集中在冠层，对粗糙子层的研究较少。然而，研究粗糙子层的粒子运动特性可以侧面提供与近地表上的物质和能量交换机制的参考，因此，研究粗糙子层中大气气溶胶粒子的运动特性是十分必要的。
本文以实际某电厂厂址周围丘陵地形模型为下垫面，在直流吹气式风洞中，采用PDA技术模拟并测定了电厂点源释放的气溶胶粒子在典型地形上粗糙子层高度的扩散行为。这为气溶胶在实际电厂周围地形上的迁移扩散提供了理论依据，也为气溶胶在实际复杂地形粗糙子层的扩散模拟提供了参考。
1 实验内容与方法
本实验是基于中国辐射防护研究院(山西太原)大气边界层二号风洞实验室完成的，该风洞目前是亚洲横截面尺寸最大的环境风洞。风洞为直流吹气式风洞，包括动力段、扩散段、稳定段、收缩段和试验段，图1为风洞的概念图、外部图和三维移测机构图。洞体总长71.1 m，试验段尺寸长28 m、宽5 m、高3.5 m，可模拟风速范围为0.2～20 m/s.
图1 风洞实验室
Fig.1 Wind tunnel laboratory
风洞模拟实验的基本原理应满足以下原则：
1)运动相似，即模型与原型有相似的大气湍流流动，这可以调节风洞入口不同数量尖劈和地面的粗糙元的摆放来实现，通过风洞湍流特性结果可得，风速廓线遵循指数规律，且二者相应高度风速比约为1∶10.
2)动力相似，即根据湍流雷诺数相等准则(即根本茂准则)，在相似的大气湍流流动中，模型与原型的来流速度之比值等于它们的模型几何缩尺比的1/3次方。结合上文，可以确定本实验模型的缩尺比为1∶1 000.
3)几何相似，即本实验的典型地形是某电厂厂址地形按1∶1 000比例缩放的简化模型。为保证再现原型厂址地形该有的细节，地形及部分构筑物采用真石漆制作，保证经多次干磨和水磨仍保持表面垂直粗糙度分辨率为1～40 μm；模型制作结果与建筑物及地形等原型缩比值的偏差不超过±2 mm.
4)扩散相似，即污染物扩散相似包括对烟气排放压力、排放口高度和排放口直径的控制来实现。
为了方便，将释放源位置坐标定为(0 mm，0 mm)，在距释放口300～800 mm处有一高90 mm的小峰，在距释放口3 200～4 300 mm处有一高308 mm的大峰(为方便下文分别简化为峰1、峰2).为研究气溶胶在粗糙子层的迁移特征，分别在两座峰的峰前、峰顶和峰后选取合适测量点。测量点选取在山体的轴线上，分布在山体前方、顶部和后方的不同位置，有代表性的位置共计6个选点，示意图如图2所示，A=(500 mm，320 mm)，B=(600 mm，320 mm)，C=(800 mm，320 mm)，D=(3 500 mm，470 mm)，E=(4 050 mm，470 mm)，F=(4 300 mm，470 mm).
图2 气溶胶扩散风洞实验测量示意图
Fig.2 Schematic diagram of aerosol diffusion experiment in wind tunnel
本试验模拟的风速是1 m/s和2 m/s，准备工作就绪后，释放气溶胶粒子进行测量。图3为实验采用的液体粒子发生器，通过其释放气溶胶DEHS(癸二酸二异辛酯，粒径范围为0.1～6 μm).对于实验中丘陵下垫面的流场特征研究测量用的是丹麦DANTEC公司生产的热线风速仪(CTA).粒子扩散测量研究采用的是丹麦DANTEC公司生产的粒子动态分析仪(PDA)，其测量原理如图4.风洞的试验段中，有三维的移动测量机构，该装置可以架设测量探头以测量数据。
图3 气溶胶粒子发生器(10E03)
Fig.3 Picture of aerosol particle source
图4 粒子动态分析技术(PDA)原理图
Fig.4 Principle schematic of particle dynamic analyzer (PDA)
2 气溶胶在粗糙子层的迁移实验结果及分析
2.1 湍流特性
为考察该风洞大气模拟动态相似条件，采用热线风速仪(CTA)研究了实验条件下大气边界层中的湍流情况。根据以下公式，可得归一化速度廓线和纵向湍强。
归一化速度廓线：
u/u0.
(1)
纵向湍强：
(2)
式中：u为局地平均速度，m/s；u0为来流上层风速度，m/s；u′为纵向脉动速度，m/s.
图5给出的是各测量位置的归一化速度廓线。选取x=100 mm处为参考点，该位置基本不受地形的影响，可以认为是平坦表面的风速廓线。由图5可以看出，在不同风速及地形条件下，参考点的风速度剖面均服从指数规律。这说明，本风洞实验符合大气模拟实验的动态相似条件。从图中还可以看出，与参考点相比，各测量点的风廓线均发生了不同程度的偏移，且随着风速的增强、峰高的增加，偏移程度加大。在粗糙子层的高度范围，风速度廓线也表现为上述变化规律，只是与冠层高度范围相比，其变化规律略弱一些。
图5 不同风速不同峰位置的归一化速度廓线
Fig.5 Profiles of mean velocity over different peaks at different wind speeds
湍流强度(Ix)是湍流流动的另一个重要参数，是脉动风速的相对强度，反映了大气湍流的运动。图6给出的是不同来流风速下的湍流强度分布。从图中可以看出，在粗糙子层的高度，来流风速及地形复杂程度对湍流强度可产生一定程度的影响。
图6 不同风速不同峰位置的湍流强度
Fig.6 Turbulence intensity over different peaks at different wind speeds
2.2 粒子纵向和垂向速度的概率分布
图7-图10中的散点表示的是，不同风速(u=1 m/s和2 m/s)时，不同粒径(0.1～2 μm、2～4 m/s和4～6 μm)的气溶胶粒子在不同峰位置的粗糙子层高度内的纵向速度及垂向速度的典型概率分布。通过对图中散点进行拟合可以看出，粗糙子层范围内粒子的纵向、垂向速度概率分布均符合高斯分布。
图7 峰1位置粒子纵向速度概率分布
Fig.7 Probability distributions of longitudinal particles velocity at peak1
从图7可以看出，对于峰1上3个测量点，均符合风速越大、拟合相关系数R2越大的规律。这可能是由于，风速越大，烟羽中夹带的粒子数越多，其拟合效果越好。另一方面，当风速为1 m/s时，概率分布曲线的峰值随粒径增加(从d=0～2 μm到d=4～6 μm)而增大。这可能是由于较大粒径的粒子较重，不易受到周围其它作用力的影响，其运动受限。当风速增大到2 m/s时，较大粒径的粒子在强风作用下运动的无序性增加，因此所有粒径段的粒子纵向速度概率分布较一致。另一方面，从概率分布数据中发现，当来流风速由u=1 m/s增加到u=2 m/s时，粒子纵向速度值范围从0.5～1.3 m/s增大到1～2.5 m/s，频率的范围从0.12～0.17减小到0.05～0.09.这表明，随着风速的增大，气溶胶粒子的动力扩散性增强[16]。
对于峰2(图8)，与峰1相比，其纵向速度概率分布曲线更高更窄，这可能是由于，峰2的位置(x=3 500、4 050、4 300 mm)距离点源较远，释放的气溶胶粒子在来流风驱动下经过长距离运输才能到达，粒子运动相对平稳。
图8 峰2位置粒子纵向速度概率分布
Fig.8 Probability distributions of longitudinal particles velocity at peak2
从图9可以看出，当u=1 m/s时，峰1上方粒子的垂向速度值几乎都是负值，即峰1上方粗糙子层高度范围中的气溶胶粒子总是向下运动的。据文献报道[12]，由于冠层高度粒子受到湍流引起的上升力作用较大，其垂向运动总是包括向上和向下两部分。然而，本文研究的高度是粗糙子层，海拔位置较高，粒子受湍流影响较小，主要在重力作用下向下运动[17]。当来风速度增加到2 m/s时，垂向粒子速度有正值出现，即表明部分粒子在较高风速下可向上运动。这一点也证实了上文的论述：风速增大，湍流特性增强，部分粒子受湍流上升力的作用而向上运动。
图9 峰1位置粒子垂向速度概率分布
Fig.9 Probability distributions of vertical particles velocity at peak1
与峰1处的概率分布曲线相比，峰2处的曲线(图10)右移，这意味着向上运动的粒子数比例增加。这是由于，湍流特性在较高的峰2上方粗糙子层继续增强，与上述讨论一致。另一方面，据文献报道，冠层粒子的垂向速度通常比纵向速度小一到两个数量级[18]。然而我们通过对比，发现在粗糙子层高度范围内，粒子的纵向速度只有垂向速度的2～5倍。
图10 峰2位置粒子垂向速度概率分布
Fig.10 Probability distributions of vertical particles velocity at peak2
2.3 粗糙子层粒子速度随高度的分布
为了研究气溶胶粒子在粗糙子层不同高度下呈现的不同运动状态，对0.1～2 μm粒子的平均纵向速度随高度的变化进行了描述，如图11所示。图中的实心点代表的是典型地形下不同下风向距离，粗糙子层内粒子的纵向平均速度，而空心点则是在平坦下垫面上相应位置的对照点。其中粒子平均纵向速度up可由下式得到：
(3)
其中upi是单个粒子的纵向速度值，m/s.
结果表明，在粗糙子层粒子的平均纵向速度基本上不随高度变化而变化。该结果与文献所报道的冠层内粒子纵向速度随高度增加而增大的结果有所不同[12]。另一方面，从图11(a)-(c)中还可以看出，风速u=1 m/s及u=2 m/s时，在同一高度下，典型地形下的粒子纵向速度在整体上均小于平坦地形的粒子速度。这说明丘陵地形的存在可减缓粒子沿下风向的运动。从局部来看，在背风坡(x=800 mm)的位置，典型地形下粒子的纵向平均速度有增大的趋势，而平坦地形的粒子速度基本没有变化。这可能是在背风侧存在的涡旋对纵向粒子产生了加速影响。当风速u=2 m/s时，上述影响更加明显。如图11(d)-(f)所示，在峰2背风侧(x=4 300 mm)，其粒子纵向速度值甚至超过了平坦地形的相应粒子速度值。这可能是由于峰2背风侧存在的涡旋更大，粒子所获的动量更大造成的。
图11 粒径0.1～2 μm的气溶胶粒子平均纵向速度随高度分布
Fig.11 Variation of mean longitudinal particle (0.1～2 μm)velocity (up)as a function of height
气溶胶粒子的垂向平均速度包括下沉速度(正值)和上升速度(负值)，如果直接对测量垂向速度数据求算术平均值是不合理的。因此我们采用公式(4)、(5)分别计算了粒子下沉平均速度vpd和上升平均速度vpa：
(4)
(5)
其中：vpdi是单个下沉粒子的速度，m/s；vpai是单个上升粒子的速度，m/s.二者随高度的变化规律如图12所示。从图12(a)-(c)可以看出，随着高度的降低，粒子平均下沉速度减小，意味着粒子越向下，其下沉运动越缓慢。这说明，虽然粒子的垂向速度是向下的，但越靠近地形，其越受到空气扰动向上的作用力影响。从图12(d)-(e)可以看出，峰2上方粗糙子层粒子的下沉速度和上升速度随高度均没有显著变化。这可能是由于，峰2粗糙子层粒子的垂向运动状态是在如空气中的粒子间碰撞、空气湍流产生的升力、Magnus效应等作用力下达到了平衡[19-20]。
图12 粒径0.1～2 μm的气溶胶粒子平均垂向速度随高度分布
Fig.12 Variation of mean vertical particle (0.1～2 μm)velocity (vp)as a function of height
3 结束语
1)在典型丘陵地形的粗糙子层高度范围，粒子的湍流特性仍受到风速和复杂地形的影响。
2)粗糙子层中气溶胶粒子的纵向、垂向速度的概率分布受粒子粒径、风速以及地形的影响。当粒子粒径从0.1～2 μm增大到4～6 μm时，其纵向概率分布曲线的拟合系数R2降低。当风速从1 m/s增大到2 m/s，纵向速度范围从0.5～1.3 m/s增大到1～2.5 m/s，频率从0.12～0.17减小到0.05～0.09.粒子垂向速度概率分布结果表明，随风速增大、峰高增加，在粗糙子层高度范围向上运动的粒子比例增大。
3)在峰1、峰2地形下，粗糙子层内粒子的纵向平均速度随高度均无明显变化。对于垂向速度，在峰1处的粒子下沉平均速度随高度的降低而减小；而峰2处粒子的上升及下沉速度基本不随高度变化而变化，粒子呈平稳的状态。
[1] PANAGIOTIS G K,ELENI V,KONSTANTINOS K,et al.Primary and secondary organic aerosol in an urban/industrial site:Sources,health implications and the role of plastic enriched waste burning[J].Journal of Environmental Sciences,2021,99(1):222-238.
[2] SHAO M,WANG S,RUSSELL A G.Air pollution complex:understanding the sources,formation processes and health effects[J].Frontiers of Environmental Science &Engineering,2016,10(5):20.
[3] LIN Y,HUANG R,YAO X.Air pollution and environmental information disclosure:an empirical study based on heavy polluting industries[J].Journal of Cleaner Production,2021,278:124313.
[4] 刘玥晨,吴志军,胡敏.二次有机气溶胶相态的研究现状与展望[J].中国环境科学,2017,37(5):1637-1645.
LIU Y C,WU Z J,HU M.Advances in the phase state of secondary organic aerosol[J].China Environmental Science,2017,37(5):1637-1645.
[5] 苏玲玲.深圳市大气颗粒物滞留时间与沉积速率多核素示踪研究[D].深圳:深圳大学,2018.
[6] 李云鹏,张俊芳,李若洁,等.中性和稳定层结条件下建筑物周围流场与污染物扩散的风洞模拟研究[J].太原理工大学学报,2021,52(3):404-410.
LI Y P,ZHANG J F,LI R J,et al.Wind tunnel simulation study of flow and diffusion around a building under neutral and stable thermal stratification[J].Journal of Taiyuan University of technology,2021,52(3):404-410.
[7] 郭栋鹏,王冉,李云鹏,等.不同温度层结下建筑物周围流场的数值模拟研究[J].太原理工大学学报,2020,51(4):572-579.
GUO D P,WANG R,LI Y P,et al.Numerical study of the effect of thermal stratification on the flow field around a building[J].Journal of Taiyuan University of Technology,2020,51(4):572-579.
[8] 滕雪冬.环境模拟试验风洞的研制[D].哈尔滨:哈尔滨工业大学,2017.
[9] PING Y,DONG Z,QIAN G,et al.Height profile of the mean velocity of an aeolian saltating cloud:Wind tunnel measurements by Particle Image Velocimetry[J].Geomorphology,2007,89(3/4):320-334.
[10] BO T L,ZHENG X J,DUAN S Z,et al.Influence of sand grain diameter and wind velocity on lift-off velocities of sand particles[J].European Physical Journal E,2013,36(5):50.
[11] LI X,DONG Z,ZHENG X.Experimental analysis of sand particles’lift-off and incident velocities in wind-blown sand flux[J].Acta Mechanica Sinica,2005,21(6):564-573.
[12] JIANG C,DONG Z,WANG X.An improved particle tracking velocimetry (PTV)technique to evaluate the velocity field of saltating particles[J].Journal of Arid Land,2017,9(5):727-742.
[13] ZHANG H,ZHANG X,LI Q,et al.Research progress on estimation of the atmospheric boundary layer height[J].Journal of Meteorological Research,2020,34(3):482-498.
[14] ZOU J,SUN J,LIU G,et al.Vertical variation of the effects of atmospheric stability on turbulence statistics within the roughness sublayer over real urban canopy[J].Journal of Geophysical Research:Atmospheres,2018.
[15] 王国羽,孙鉴泞.城市粗糙子层湍流能谱特征分析[J].南京大学学报(自然科学版),2014,50(6):820-828.
WANG G Y,SUN J N.Characteristics of turbulence spectra in the urban roughness layer[J].Journal of Nanjing University(Natural Sciences),2014,50(6):820-828.
[16] XIAO F,DONG Z,GUO L,et al.Sand particle lift-off velocity measurements and numerical simulation of mass flux distributions in a wind tunnel[J].Journal of Arid Land,2017,9(3):331-344.
[17] BAUER B O,HOUSER C A,NICKLING W G.Analysis of velocity profile measurements from wind-tunnel experiments with saltation[J].Geomorphology,2004,59(1/4):81-98.
[18] DONG Z,LIU X,WANG X,et al.Experimental investigation of the velocity of a sand cloud blowing over a sandy surface[J].Earth Surface Processes &Landforms,2010,29(3):343-358.
[19] LIU L,HU F,CHENG X L.Extreme fluctuations of vertical velocity in the unstable atmospheric surface layer[J].Nonlinear Processes in Geophysics,2014,21(2):463-475.
[20] RICE M A.An experimental study of multiple grain-size ejecta produced by collisions of saltating grains with a flat bed[J].Sedimentology,2010,42:695-706.