考慮羽流浮力影響的高層建筑周圍污染物擴散數(shù)值模擬

周晅毅，馬慧心，顧 明

（同濟大學土木工程防災(zāi)國家重點實驗室，上海 200092）

近年來，隨著城市化和工業(yè)化進程的不斷推進，環(huán)境污染問題引起了公眾的廣泛關(guān)注。建筑周圍氣體污染物的釋放與擴散給城市環(huán)境質(zhì)量及居民生命健康帶來了極大危害。對氣體污染物的擴散進行正確模擬，包括對時均濃度以及瞬時濃度的預(yù)測，可以為制定相應(yīng)的應(yīng)急管理措施提供重要參考信息。

由于羽流的運動會受到風和建筑結(jié)構(gòu)之間復(fù)雜相互作用的影響，建筑周圍污染物的擴散過程將會受到其自身的物理特性、大氣邊界層穩(wěn)定性、周邊建筑形態(tài)等因素的影響。在建筑后方尾流區(qū)等弱風區(qū)域，由環(huán)境空氣與污染物之間的密度差異所產(chǎn)生的羽流浮力效應(yīng)是驅(qū)動污染物擴散的重要因素。輕質(zhì)污染物由污染源處釋放后，由于其密度小于周圍空氣，所受到向上的浮力大于自身重力，從而產(chǎn)生了向上的合力，這里稱之為正浮力；而重質(zhì)污染物密度大于空氣，向上的浮力小于自身重力，產(chǎn)生了向下的合力，這里稱之為負浮力；正負浮力均會對污染物的濃度分布產(chǎn)生重要影響。因此，研究密度差異所引起的浮力效應(yīng)影響對污染氣體輸運的預(yù)測具有重要意義。

目前，建筑周圍污染物擴散的主要研究方法有以下三種：現(xiàn)場實測、風洞試驗和數(shù)值模擬［1］?，F(xiàn)場實測可以提供實際情況下氣體擴散的第一手資料，但受氣象條件等因素的影響，現(xiàn)場實測的結(jié)果通常難以得到重現(xiàn)。風洞試驗具有試驗邊界條件可控的優(yōu)勢，但對于考慮浮力效應(yīng)的湍流擴散過程，其在滿足密度弗勞德數(shù)的相似準則時需要較低的風速，這給滿足雷諾數(shù)獨立性的假定帶來了一定的困難［2］?；谟嬎懔黧w動力學理論的數(shù)值模擬是研究污染物擴散問題的另一有力工具，目前已被廣泛采用［3］。相對于現(xiàn)場實測和風洞試驗，數(shù)值模擬方法成本較低且研究周期短［4］；同時，可以克服風洞試驗受相似準則限制的弊端，并能夠提供細致的流場及濃度場數(shù)據(jù)，在現(xiàn)象分析及機理研究方面具有優(yōu)越性。常用的數(shù)值模擬方法主要包括雷諾平均方法（Reynolds-averaged Navier-Stokes simulation，RANS）和大渦模擬（large-eddy simulation，LES）。Olvera等［5］采用定常RANS方法對從立方體建筑后方釋放的中質(zhì)及輕質(zhì)羽流的擴散進行了數(shù)值模擬，發(fā)現(xiàn)正浮力會在尾流區(qū)引起較大的流動擾動，從而使建筑背風面的濃度較高，而建筑下游地面處的濃度較低。相對于單獨討論正負羽流浮力的作用，在相同幾何條件及模型設(shè)置下將二者同時考慮更有利于對比分析其對流場及濃度場的影響機理，進一步總結(jié)出相應(yīng)的規(guī)律。目前，同時考慮正負羽流浮力影響的研究僅有Tominaga和Stathopoulos［6］針對立方體建筑周圍的擴散問題所進行的模擬，該研究采用RANS方法進行，指出RANS方法的數(shù)值特性限制了建筑后方湍流及氣體污染物濃度分布的預(yù)測精度，但缺少對不同羽流浮力影響下濃度輸運通量及濃度瞬態(tài)變化特性的分析。

隨著計算資源的不斷進步，近年來大渦模擬方法得到了廣泛關(guān)注［7］。LES方法直接求解大尺度漩渦，僅需要對各向同性的小尺度漩渦采用模型進行處理。相較于RANS方法而言，前人研究［8-10］發(fā)現(xiàn)在羽流浮力可以忽略的情況下LES方法對污染氣體擴散的預(yù)測精度更高。除此之外，LES 方法還可以提供羽流擴散的重要瞬態(tài)信息。但是，目前采用LES方法討論羽流浮力影響的文章相對較少。Wingstedt 等［11］對中質(zhì)和重質(zhì)氣體的擴散進行了大渦模擬，發(fā)現(xiàn)與中質(zhì)氣體相比重質(zhì)氣體的羽流外形會更寬更淺，且受下游立方體陣列高度的影響更大。Lin 等［12］以單一立方體建筑為研究對象，比較了RANS和LES方法對受正浮力影響的氣體濃度分布的預(yù)測精度。

如前所述，前人對于浮力效應(yīng)影響下污染物擴散的研究大多以立方體建筑為研究對象且側(cè)重于數(shù)值方法的對比驗證，目前很少有針對高層建筑周圍浮力氣體擴散的分析；同時，關(guān)于浮力對時均及脈動濃度場、污染物標量輸運機制影響的討論并不充分。因此，本文同時考慮正負浮力的影響，采用LES方法模擬了一個高層建筑周圍的流場及濃度場分布，總結(jié)了不同羽流浮力影響下高層建筑周圍污染物時均濃度場的分布特征，并對羽流浮力對濃度輸運機制、濃度脈動的影響進行了探討。

1 模擬案例概要

本文中的研究對象根據(jù)東京工藝大學的風洞試驗?zāi)Ｐ停?3］布置，如圖1所示。

圖1 結(jié)構(gòu)布置示意圖Fig.1 Schematic diagram of the building model

高層建筑模型置于中性大氣邊界層中，建筑模型沿x、y、z坐標方向的尺寸分別為0.08 m（B）×0.08 m（B）×0.16 m（H）。污染源直徑?為5 mm，其圓心距離建筑背風面的距離為0.25H。來流垂直于建筑的迎風面。為了研究正負羽流浮力對高層建筑周圍氣體污染物擴散的影響，設(shè)置了如表1所示的工況。這里，量綱一參數(shù)密度弗勞德數(shù)（Fr）按照公式（1）定義，用以表征污染物所受到的浮力與環(huán)境空氣的慣性力之間的相對大??；ρs和ρa分別為污染物密度和空氣密度，Δρ=ρs-ρa為污染物與空氣之間的密度差值，g為重力加速度。各工況的入口風剖面和湍動能剖面由試驗數(shù)據(jù)（EXP）擬合得到，如圖2中的實線所示。其中，工況1 在建筑高度處的順流向平均風速UH=1.40 m·s-1，與風洞試驗保持一致，用以驗證數(shù)值模擬的準確性。工況2～工況4中UH均為0.60 m·s-1，根據(jù)環(huán)境空氣與污染物之間的密度差異，有三個不同的密度弗勞德數(shù)，用以討論氣體污染物羽流浮力的影響。污染源處的氣體總釋放率為q=2.62×10-4UHH2，m3·s-1，對于工況1和工況2～工況4，污染源處所釋放的污染物體積分數(shù)F分別為5%和100%。

表1 數(shù)值模擬工況Tab.1 Cases of numerical simulation

圖2 數(shù)值模擬入口剖面Fig.2 Inflow profiles of numerical simulation

在本文中，風洞試驗及數(shù)值模擬的風速均采用參考風速UH進行量綱一處理；體積濃度c用來表征污染物的濃度分布，其與污染物質(zhì)量分數(shù)φ之間的關(guān)系如公式（2）所示

體積濃度數(shù)據(jù)均采用參考濃度C0進行量綱一處理，C0按照公式（3）定義：

式中：Qe=qF為污染源處污染物的體積流率。

2 數(shù)值方法與計算設(shè)置

本文采用LES方法對不同羽流浮力影響下的氣體擴散問題進行數(shù)值模擬，使用商業(yè)軟件ANSYS Fluent進行模擬。

2.1 大渦模擬方法

LES方法采用濾波函數(shù)對大尺度渦和小尺度渦進行區(qū)分，大尺度渦屬于可解尺度，直接進行求解；小尺度渦則需要采用模型處理。濾波后的質(zhì)量、動量守恒和組分輸運方程如式（4）—（6）所示：

式中：“·～”表征可解尺度變量；ui為速度分量；xi為空間坐標；t為時間；ρ為流場中氣體污染物和空氣所形成混合物的密度；μ為流體動力粘度；p為壓力；g為重力加速度；δi3為Kronecker delta 函數(shù)；φ為污染物的質(zhì)量分數(shù)；D為分子擴散系數(shù)。

亞格子尺度雷諾應(yīng)力τSGSij和亞格子尺度質(zhì)量通量JSGSi按照公式（7）—（8）求解：

式中：μSGS為亞格子湍流粘性，需要通過亞格子模型處理；ScSGS為亞格子湍流施密特數(shù)，在本文中設(shè)置為0.7。

本文所采用的亞格子湍流模型為動態(tài)Smagorinsky-Lilly 模型，文獻［14-16］已經(jīng)證明該模型是建筑環(huán)境周圍流動及污染物擴散模擬的合適選擇。各工況的計算時間步長如表1 所示，均滿足庫朗數(shù)（Courant-Friedrichs-Lewy numbers）小于1.0 的要求。各工況的量綱一采樣時間（t*=t·UH/H）均為840，經(jīng)驗證明該采樣時間足以獲得統(tǒng)計意義上收斂的變量時均值。

2.2 計算域及邊界條件設(shè)置

數(shù)值模擬的計算域如圖3a所示，其大小x·y·z為12.5H×7.5H×6.25H，計算域的坐標原點設(shè)置在建筑底面中心處。入口邊界設(shè)置為速度入口（Velocity inlet），平均風剖面及湍動能剖面如圖2 中的擬合曲線所示，湍流耗散率根據(jù)日本建筑學會（AIJ）行人風環(huán)境數(shù)值模擬指南［17］中的規(guī)定確定。入口處的速度脈動采用渦方法（Vortex method）［18］生成，渦的數(shù)目取為190［16］。出口邊界條件為自由出流（Outflow），地面及建筑表面設(shè)置為無滑移壁面邊界（No-slip wall），采用Kader［19］提出的混合函數(shù)進行近壁面處理；頂面及側(cè)面設(shè)置為對稱邊界（symmetry）；污染源設(shè)置為速度出口（Velocity inlet）。速度與壓力耦合采用SIMPLEC 算法，動量及組分輸運方程的離散采用有界中心差分格式（bounded central differencing scheme），壓力項的離散采用二階格式（second order），時間離散采用有界二階隱式格式（bounded second order implicit）。

圖3 計算域及網(wǎng)格示意圖Fig.3 Schematic diagram of computational domain and grid.

3 數(shù)值方法驗證

為了驗證所采用的數(shù)值方法對高層建筑周圍風場及污染物濃度場模擬的準確性，圖4將工況1中心平面（y/H=0）上建筑后方的時均順流向風速及時均濃度的模擬值與日本東京工藝大學的風洞試驗結(jié)果［13］進行了對比。由圖中可以看出，LES 方法對風場的模擬結(jié)果與實驗結(jié)果呈現(xiàn)出較好的一致性；對于平均濃度分布，數(shù)值模擬基本重現(xiàn)出了建筑后方污染物的分布特征，但在建筑下方（z/H＜0.5）LES的預(yù)測結(jié)果偏小。文獻［20］對該試驗的大渦模擬研究也得到了和本文相似的結(jié)果，在風洞試驗中，靠近地面位置的濃度測量存在一定的難度及不確定性，這可能是誤差產(chǎn)生的重要原因之一?？偟膩碚f，本文所構(gòu)建的數(shù)值模型可以較好地重現(xiàn)高層建筑周圍的流場及濃度場分布，可用于進一步討論羽流浮力效應(yīng)對污染物擴散的影響。

圖4 數(shù)值方法驗證Fig.4 Validation of numerical method

4 結(jié)果與討論

4.1 平均流場

圖5給出了各工況中心平面（y/H=0）處的時均垂直速度分量及流線分布。圖中三角形表示污染源所在的位置，在后文的所有圖片中均有相同的含義。流線圖直觀展示了建筑后方回流區(qū)的大小，對于本文所討論的高層建筑，其后方回流區(qū)的長度為0.6H。整體上，羽流浮力對流場的影響較?。坏诳拷ㄖ筹L側(cè)及污染源處，輕質(zhì)氣體所受到的正浮力使垂直速度分量明顯增大，較大的上升氣流使回流區(qū)的流線有明顯的上洗現(xiàn)象；而重質(zhì)氣體工況由于向下負羽流浮力的存在，氣流的爬升受到抑制，產(chǎn)生了相反的效果。

4.2 平均濃度場

圖6為各工況在中心平面（y/H=0.1）以及z/H=0.1 平面處的時均濃度分布云圖。由于污染源位于建筑后方尾流區(qū)內(nèi)（x/H=0.5，y/H=0），在尾流渦的作用下，污染物自污染源處釋放后朝著建筑背風墻面輸運，于是在建筑背風墻面和污染源之間有較高的濃度。對于輕質(zhì)氣體工況，污染物羽流中心線與地面呈一個明顯夾角，正浮力和尾流渦的共同作用使建筑背風側(cè)產(chǎn)生了較強上洗氣流（圖5 a），因此污染物沿建筑背風墻面爬升并產(chǎn)生了較高的濃度（圖6 a）。相較于無浮力及負浮力工況，在正浮力的作用下更多的污染氣體突破了尾流區(qū)的限制，對建筑頂面產(chǎn)生了較大影響。對于重質(zhì)氣體工況，受向下的負羽流浮力的影響，污染物沿建筑背風墻面的爬升被明顯削弱，在建筑下方（z/H＜0.5）積聚形成了高濃度區(qū)域并沿著順風向輸運。同時，在z/H=0.1平面處，重質(zhì)氣體沿橫風向的擴散較輕質(zhì)氣體有明顯的增強。

圖5 中心平面（y/H=0）處量綱一垂直速度分量（W/UH）云圖及流線分布Fig.5 Distribution of time-averaged streamlines at vertical plane y/H=0 contoured by the time-averaged dimensionless vertical velocity W/UH.

圖6 量綱一時均濃度（C/C0）分布云圖Fig.6 Contours of dimensionless time-averaged concentration (C/C0)

為了進一步評價羽流浮力對高層建筑后方污染物濃度累積情況的影響，圖7給出了三個工況在建筑后方沿順流向若干截面的平均濃度分布情況，各個截面的尺寸均與建筑背風墻面相等。在建筑后方回流區(qū)內(nèi)（0.25＜x/H＜0.85），受不同羽流浮力的影響，三個工況的面平均濃度差異顯著。在靠近建筑背風面處，輕質(zhì)氣體工況的截面平均濃度最高，重質(zhì)氣體工況的平均濃度最低，這主要是由于在正浮力的作用下較多的污染物沿建筑背風面爬升所導(dǎo)致的。然而，隨著截面位置沿順流向推移，這種現(xiàn)象發(fā)生了反轉(zhuǎn)，重質(zhì)氣體面平均濃度明顯高于其他兩種氣體，尤其是在污染源所處的平面上（x/H=0.5）。在回流區(qū)下游（x/H＞0.85），三工況之間的面平均濃度差異沿順流向逐漸減小，但輕質(zhì)氣體始終保持最低的濃度積累水平?？梢?，負浮力和正浮力分別抑制和增強了建筑后方，尤其是回流區(qū)內(nèi)的氣體擴散。

圖7 建筑后方面平均濃度（C/C0）分布Fig.7 Plane-averaged dimensionless time-averaged concentration(C/C0)behind the building

4.3 濃度輸運

氣體污染物濃度的標量輸運受到平均流動及湍流的共同影響，二者對濃度輸運所做的貢獻通常用對流通量（UiC）和湍流通量（）表示。濃度通量的分布反映了濃度遷移的本質(zhì)，因此分析不同羽流浮力工況的對流通量和湍流通量，對進一步解釋羽流浮力對污染物擴散機理的影響是至關(guān)重要的［8］。圖8給出了三個工況在中心平面（y/H=0）內(nèi)的量綱一對流及湍流通量分布。由圖中可以看出，各個工況的對流通量和湍流通量在量級上基本一致，也就是說平均流動和湍流對污染物在豎直方向上的擴散貢獻相當，這種現(xiàn)象也曾被Gousseau等［10］在針對不受浮力影響的污染物擴散問題研究中觀察到。圖8反映出，在正浮力的影響下靠近建筑背風側(cè)的對流及湍流通量均有明顯的增強，使得污染物更多的向上方輸運；而負浮力則產(chǎn)生了相反的效應(yīng)。圖9定量給出了三個工況在建筑背風面后方建筑高度處的對流及湍流通量分布，在相對遠離建筑背風側(cè)的x/H＞0.375位置處，三個工況的對流通量（圖9 a）均出現(xiàn)了負值，這主要是由于尾流渦的存在（圖5）而導(dǎo)致的。在同一位置處，三個工況的湍流通量始終為正值，即向上輸運；且輕質(zhì)氣體工況的湍流輸運強度明顯大于中質(zhì)及重質(zhì)工況，因此，輕質(zhì)氣體在湍流作用下更多的從尾流區(qū)擴散出去，而重質(zhì)氣體則在尾流區(qū)內(nèi)出現(xiàn)了堆積，如圖6所示。

圖9 建筑高度處垂直方向污染物濃度通量對比Fig.9 Comparison of vertical concentration fluxes of pollutant at the building height

4.4 濃度脈動

對于會對人們生命健康產(chǎn)生不利影響的有害氣體，除平均濃度外，瞬時濃度也是需要關(guān)注的重點。根據(jù)Li 和Meroney［21］的定義，圖10 給出了不同羽流浮力影響下建筑后方的局部濃度脈動強度ICL局部濃度脈動強度是將濃度的均方根值使用當?shù)氐钠骄鶟舛冗M行歸一化處理，在羽流浮力的影響下，三個工況呈現(xiàn)出明顯不同的分布特征。對于輕質(zhì)工況，污染源正上方的ICL明顯大于中質(zhì)工況，也就是說由于正浮力的作用導(dǎo)致該處出現(xiàn)了較平均濃度而言很大的濃度脈動。重質(zhì)工況的濃度脈動明顯比中質(zhì)工況小很多，且ICL的較大值分布在污染源順流向的后方。值得一提的是，輕質(zhì)工況的ICL峰值超過了4.0，而重質(zhì)工況的ICL峰值僅在2.5左右?？偟膩碚f，輕質(zhì)氣體所受到的正羽流浮力使?jié)舛让}動進一步增強，而重質(zhì)氣體受到的負羽流 浮力則抑制了濃度脈動。

圖11 給出了各個工況在監(jiān)測點位置處的脈動濃度概率密度分布（probability density distribution，PDF）圖。這里共選取了三個監(jiān)測點，其中A 點（0.05，0，0.128）和B點（0.05，0，0.03 2）靠近受羽流浮力影響較大的建筑背風側(cè)分布，C 點（0.2，0，0.032）布置在尾流渦之外，各個監(jiān)測點的位置如圖11所示。由圖中可以看出，在三個監(jiān)測點位置處，不同羽流浮力影響下各個工況的脈動濃度均不服從正態(tài)分布。表2中對各個監(jiān)測點處污染物濃度的統(tǒng)計特性進行了總結(jié)，σc為污染物瞬時濃度的均方根值，這里以污染物濃度時程數(shù)據(jù)的95%分位值作為峰值濃度（CP）。輕質(zhì)氣體工況在靠近建筑背風側(cè)的A、B兩點處相對于其他兩個工況而言具有較高的平均濃度，同時也具有較高的峰值濃度，但濃度正值波動的出現(xiàn)概率相對較低。而遠離建筑背風側(cè)的C點則情況相反，在該位置處重質(zhì)氣體工況具有較高的平均濃度和峰值濃度，以及相對較小的濃度正值波動概率?？偟膩碚f，各個監(jiān)測點處的峰值濃度可達平均濃度的2～3 倍，對于有毒有害氣體的防控應(yīng)該引起特別的注意。

圖11 監(jiān)測點處脈動濃度概率密度分布圖Fig.11 Probability density distribution of fluctuating concentration at three key points.

表2 污染物濃度統(tǒng)計特性Tab.2 Statistical characteristic of pollutant concentration.

5 結(jié)論

本文對高層建筑周圍的污染物擴散問題進行了LES模擬研究，設(shè)置三個工況同時考慮了由污染物和周圍空氣的密度差異所產(chǎn)生的正負羽流浮力的影響。通過將中質(zhì)氣體工況下的模擬結(jié)果與風洞試驗數(shù)據(jù)對比驗證了數(shù)值模擬方法的準確性，并在此基礎(chǔ)上進一步分析了不同羽流浮力對時均濃度、濃度輸運以及濃度脈動的影響，得出的主要結(jié)論如下：

（1）輕質(zhì)氣體所受正浮力使污染物沿建筑背風面的爬升增強，從而導(dǎo)致靠近建筑背風側(cè)濃度較高；而重質(zhì)氣體受負浮力的影響在建筑下方堆積，沿順流向及橫風向發(fā)生明顯擴散?？偟膩碚f，正負浮力分別增強和抑制了建筑后方尾流區(qū)內(nèi)對污染物的稀釋作用。

（2）輕質(zhì)氣體工況靠近建筑背風側(cè)的豎直方向?qū)α骷巴牧魍肯噍^于其他兩個工況均有明顯的增強，因此更多的污染物突破了尾流渦的限制，從上方擴散出去；而重質(zhì)工況的對流及湍流通量均較小，導(dǎo)致其在尾流區(qū)內(nèi)發(fā)生了堆積。

（3）輕質(zhì)氣體所受正浮力增強了靠近污染源及建筑背風側(cè)的濃度脈動，而重質(zhì)氣體所受的負浮力則對濃度脈動有一定的抑制作用。對于不同浮力工況而言，靠近建筑背風側(cè)及尾流區(qū)外的監(jiān)測點處峰值濃度均可達到平均濃度的2～3倍。

本文僅對弗勞德數(shù)為±3.5的兩個浮力氣體擴散工況進行了探討，但在實際的氣體污染物擴散過程中根據(jù)羽流浮力強弱的不同可能涉及到更大的弗勞德數(shù)范圍，在后續(xù)的研究工作中將設(shè)計補充工況，從而總結(jié)得出適用范圍較大的一般性結(jié)論。本文針對高寬比為1：2的單體高層建筑進行了研究，而對于高寬比較小的超高層建筑、建筑群中受到環(huán)境干擾的建筑，由于受到復(fù)雜流動結(jié)構(gòu)的影響，污染物的擴散特性也會有所不同，這需要在后續(xù)的工作中做進一步的討論。除此之外，風洞試驗數(shù)據(jù)的缺乏是制約數(shù)值模擬精度驗證的一個問題，后續(xù)研究工作中計劃針對高層建筑周圍的浮力氣體擴散進行風洞試驗研究，獲得可靠的脈動濃度數(shù)據(jù)，進一步驗證數(shù)值模擬的精度。