城市大气微环境大涡模拟研究

中国科学G辑:物理学力学天文学2008年第38卷第6期:1~12城市大气微环境大涡模拟研究崔桂香①*,史瑞丰①,王志石②,许春晓①,张兆顺①①清华大学工程力学系,北京100084;②澳门大学科技学院,澳门999078;*E-mail:[e-mail]收稿日期:2007-10-30;接受日期:2008-03-11国家自然科学基金(批准号:10572073)(NSFC-10572073)和澳门科技发展基金(批准号:0221/Z006/A)(FDCT022/2006/A)资助项目摘要良好的空气质量是居民生活环境的重要条件之一,大气微环境研究与居民区的空气质量密切相关.大气环境由多尺度的空气运动主宰,而城市居民区空气流动属于微尺度大气流动.介绍城市居民区微尺度大气环境的现代数值方法——大涡数值模拟方法.针对居民小区复杂流动的特点,讨论如何正确应用大涡数值模拟方法,比较不同数值方法对预测城市微环境流动的影响,如:下垫面阻力元法与浸没边界法,亚格子模式的影响等.以模型小区和实际小区为例计算了风场(包括湍流统计特性),污染物浓度等物理量分布,数值计算结果和模型小区及实际小区的实验测量数据符合良好.关键词微环境流动大涡模拟湍流扩散良好的空气质量是优美人居环境的主要条件之一.随着城市化和人民生活水平的提高,城市居民区的空气质量愈来愈受到重视.随着城市规划科学性和人居环境质量要求的不断提高,传统的经验预测或工程估计方法难以满足要求.在计算机资源迅速改善的今天,用计算流体力学的新方法研究环境流动已成为经济有效的手段.然而,实现大气微环境的准确大涡模拟预测,仍需要在流动的物理机制和湍流模拟方法两方面进行深入研究.城市居民小区的大气环境空间范围较小,属于微尺度气象问题,也是城市大气环境的子系统[1],具有复杂的湍流流动特征,需要较准确地分辨建筑物等几何边界,才能准确预测风场和污染物浓度分布.另一方面,居民小区的大气环境是大气边界层内的流动现象,数值预测与中尺度大气预测有不同的特点,随机扰动与小尺度流动结构都可能对微环境流动产生较大影响,必须尽可能精细地考察区域内流动与物理量变化.因此,不能把下垫面简化为粗糙表面,而需要将建筑物等组成的下垫面作为求解域边界,准确满足边界条件.居民小区大气环境处在大气边界层内,而大气边界层的流动是非定常的湍流运动,因此,崔桂香等:城市大气微环境大涡模拟研究2居民小区大气环境的数值预测属于非定常复杂湍流的数值计算.湍流数值模拟是流动数值模拟的前沿课题,由于湍流含有不同尺度的脉动,通常称之为各种尺度的“湍涡”.分辨所有尺度湍涡的数值模拟称为直接数值模拟(DirectNumericalSimulation,DNS),但是目前计算机水平不足以分辨复杂湍流中所有尺度的湍涡.工程计算中常用时间平均法,将所有湍流脉动对平均运动的作用归结为雷诺应力和平均湍流扩散,进而对雷诺应力和平均湍流扩散作模式假定,这种方法称作雷诺统计平均法(ReynoldsAverageNavier-Stokes,RANS).显然,RANS只能提供时间平均量,不能提供脉动量,难以准确地提供大气运动的演化过程.大气科学家[2]最早提出一种湍流大涡数值模拟方法(LargeEddySimulation,LES),这种方法的基本思想是将湍流运动中的小尺度脉动过滤掉,将小尺度脉动对大尺度(大涡)的作用建立模型.近年来,这种方法的研究和应用取得很大发展.本文先简要介绍大气环境数值模拟方法的基本原理,然后结合具有实验结果的模型小区和实际居民小区的算例,说明微环境流动的大涡数值模拟方法和相应的结果.1数值方法1.1大涡模拟的控制方程大气运动属于低速空气运动,它的控制方程满足Navier-Stokes方程和质量扩散方程,大涡模拟的控制方程则通过过滤运算得到,本文采用以下过滤函数()()()d,iiuuGV=−∫xyxy(1)()()1231,G=ΔΔΔη2,iiηΔ≤()0,G=η|ηi|Δi/2,式中Δi是过滤长度,将(1)式用于Navier-Stokes方程和质量扩散方程,便可得到如下大涡模拟的控制方程组:0,iiux∂=∂(2)2,ijiiijijijjjuuupuftxxxxxτνρ∂∂∂∂∂+=−+++∂∂∂∂∂∂(3)2.jjcjjjjTcccuDStxxxx∂∂∂∂+=++∂∂∂∂∂(4)方程(2)是质量守恒方程.方程(3)是动量方程,其中ν和D分别为分子黏性系数和分子扩散系数,if表示外力,例如,包括浮力,方程中除了分子黏性应力外(右边第2项),还有过滤掉的小尺度湍涡对大涡的作用τij,称作亚格子应力;方程(4)是污染物浓度的守恒方程,右边是分子扩散项、亚格子质量通量项和污染物的源项cS,例如机动车或电厂的排放.在大涡数值模拟中亚格子应力和亚格子质量通量都是未知量,需要用湍流模型予以封闭.在湍流理论中有各种亚格子模式[3],复杂流动中常用的模式是涡黏和涡扩散模式.中国科学G辑:物理学力学天文学2008年第38卷第6期31.2亚格子模式方程(3)和(4)中τij与Tj分别是亚格子应力和亚格子质量通量,定义为,ijijijuuuuτ=−(5).jjjTcucu=−(6)本文采用两个适于复杂边界湍流计算的涡黏和涡扩散模式,一个是Smagorinsky模式,表达式可写作:212,,3ijtijkkijtSCSτντδν=+=Δ(7),jtjTxθκ∂=∂2tDSκ=Δor,tttPrνκ=(8)式中νt和κt分别称为亚格子涡黏系数和涡扩散系数,()122,ijijSSS=()/2,ijijjiSuxux=∂∂+∂∂方程(7)和(8)中C和D为是模式常数,Prt为湍流普朗特数,可由模式常数计算Prt=C/D,最简单的亚格子模型是Smagorinsky模式,常数2,sCC=0.18,sC=此时,湍流普朗特数可取作Prt=0.7~0.8.大量大涡模拟计算证实Smagorinsky模式耗散过大,因此计算中常采用减小模式系数以克服耗散过大的缺点,通常可取0.08~0.1.sC=为了提高模拟精度,准确地模拟亚格子应力对可解尺度物理量的影响,采用包含局部湍流特性的涡黏模型,动态确定模型常数C[4],这是本文采用的另一个亚格子模式,即Lagrange动力模式:,ijijijijCMLMM=(9)(9)式中k,ijijijLuuuu=−k22212,ijijijMSSSS⎡⎤=Δ−Δ⎢⎥⎣⎦Δ1和Δ2是两次过滤的尺度,且有Δ1Δ2,本文计算中取Δ2=2Δ1,上标‘−’和‘~’分别表示由Δ1和Δ2两次过滤的量.尖括号表示系综平均,可在均匀空间体积、均匀面或线上平均.然而,复杂湍流很难找到均匀面进行平均,因此,本文采用Meveneau[5]提出的在流体质点轨迹上统计平均的方法,具体公式如下:,LMMMCII=(10)其中ILM和IMM分别是ijL和ijM在迹线上积分并插值到网格点上的值,最终表达式为:()()()()()111,nnnnnnLMLMijijLMIIttLMIT+++−−ΔΔ=−⎡⎤⎡⎤⎢⎥⎣⎦⎣⎦xxuxx(11)()()()()()111,nnnnnnMMMMijijMMIIttMMIT+++−−ΔΔ=−⎡⎤⎡⎤⎢⎥⎣⎦⎣⎦xxuxx(12)式中()181.5LMMMTII−=Δ是根据实验的近似计算公式,上标n和n+1表示时间步.同样方法可导出(8)式中的涡扩散系数D.1.3数值格式与方法有限体积法是流体数值计算领域熟知的离散化方法[6],具有很好地守恒性,但经典的有限崔桂香等:城市大气微环境大涡模拟研究4体离散格式精度只有二阶,为了提高精度,加大时间步长,本文采用非交错网格、Padé4阶紧致插值的高精度有限体积法,时间方向用4阶龙格库塔积分,详细计算方法可参见徐岚的博士论文[7].计算域边界上合理的物理条件和正确的数值处理方法是获得准确模拟结果的基本保证,可以根据数学原理给定方程有解的边界条件(数学上称为适定条件),但是要满足实际情况的边界条件并不容易,如:下垫面需要准确满足无滑移边界条件,由于城市居民区有各种各样几何形状的建筑、街道等,要满足无滑移条件并非易事.边界条件数值方法的难点在于,为了准确满足复杂几何边界的无滑移条件,需要生成贴体网格,其难度和计算量都非常大,对于运动边界,每个时间步都要重新生成网格,计算量就更大,考虑到环境问题中下垫面的复杂程度,在目前的计算条件下采用贴体网格几乎是不可能实现的.因此,本文应用新的数值技术——浸没边界方法(ImmersedBoundaryMethod,IBM),基本思想是在动量方程中增加强迫力项,使之满足无滑移边界条件:2,ijijiiiijijjjuuuupfgtxxxxxτνρ∂∂∂∂∂=−−++++∂∂∂∂∂∂(13)式中gi只有在包含边界点的网格上不等于零,它的作用是使固壁上满足无滑移条件.gi的给法有不同形式,本文采用直接强迫力方法如下:(),nnnibiigRHSVut=−+−Δ上标n表示当前时间步.RHS包含对流项、压力项等方程(13)的右端项,Vbi是物面的速度,于是边界附近的网格点xi上n+1时刻速度便可由下式计算:()1,nnnniiiuutRHSg+−=Δ+(14)IBM的优点在于可采用规则的直角坐标系网格实现复杂边界的无滑移条件,可采用与求解域内点一致的解法求解,既满足无滑移条件,又使网格简单方便.为了考察IBM方法的优劣,本文还在同一算例中采用环境流动计算的平均阻力元方法进行比较.本文数值求解是在MPI并行系统上,采用分区并行的方法,求解离散化的方程组,分区界面采用重叠网格迭代的方法保持解的连续性和准确性.2模型小区大气环境的大涡数值模拟2.1计算条件本算例求解域为一模型小区建筑群,由39个立方体交错排列组成,其分布如图1所示,每个立方体长宽高均为0.12m,立方体间距为0.24m.本文计算条件与Davidson[8]实验条件一致,入流风速U0=4.0m/s,以入流风速为特征速度,建筑物高H为特征长度的雷诺数Re=24300.实验中在第1排建筑物前1.2m处释放二氧化碳气体作为污染源,模型实验在风洞中完成,实验过程中测量建筑群间的气流速度和二氧化碳浓度.本文算例以主流风向为流向(计算域长度),水平面上垂于主流风的方向为横向(计算域宽度),与水平面垂直的方向为垂向(计算域高度),位置坐标分别记作x1,x2,x3.计算域的大小由中国科学G辑:物理学力学天文学2008年第38卷第6期5小区面积决定,计算域高度取决于边界层的厚度,一般可根据地面障碍物高度试算,以计算域顶部流动接近自由流为准,即顶部流动不受地面障碍物的影响,本文计算域高度取立方体高度的10倍,它大于风洞实验的边界层厚度0.1m;宽度和长度分别取立方体高度的30倍和70倍.采用的网格数为162×171×30,在地面和建筑物附近进行了局部加密,如图2所示.该网格已通过网格无关性检验.图1建筑群模型(a)俯视,(b)侧视图2计算网格(a)俯视,(b)侧视动量方程的边界条件为:侧向边界提充分发展的零梯度条件;出流边界采用无反射条件,如出流边界在x1方向,可写作1//0iiiutuux∂∂+⋅∂∂=;顶部给镜像条件;地面和建筑物表面等底部边界上无滑移条件,用浸没边界法实现数值计算;值得注意的是主流风的入流边界上,即,计算域迎风边界上的湍流脉动是很重要的参数,因为小区风场处于大气边界层的内层,气流脉动比较剧烈,能否给出接近实际问题的入流湍流脉动,直接关系到对湍流输运模拟的准确性,即会影响温度或质量的湍流扩散的模拟结果.本文根据文献[8]的物理实验给定平均风剖面与随机脉动之和,表达式如下:平均剖面为:()*0lninUUzzκ=(*U为表面摩擦速度)00,0,0.0025,0.4,ininVWzHκ====(15)0inUU=(0U为是顶部自由流风速),0,0,8,ininVWzH==≥(16)脉动量