不同溫度層結下建筑物周圍流場的數(shù)值模擬研究

郭棟鵬，王 冉，李云鵬，趙 鵬，姚仁太

(1.太原科技大學 環(huán)境與安全學院，太原 030024；2.中國輻射防護研究院，太原 030006)

有潛在風險的污染物通過建筑物頂部釋放進入大氣環(huán)境中，其遷移擴散可能會對民眾的健康、舒適性造成潛在風險[1]。影響污染物遷移擴散的因素有很多，如建筑物周圍風速、排放高度、溫度層結等，其中，溫度層結會影響大氣邊界層(ABL)的厚度、結構以及邊界層內的大氣速度、溫度和湍流廓線，它們對于實際大氣中污染物擴散問題起主要作用。

關于此類問題，最初的探索是在風洞中實現(xiàn)的(MERONEY et al[2]，UEHARA et al[3]，YASSIN et al[4-5])，研究結果表明，溫度層結對建筑物周圍空氣流動與擴散產生影響。伴隨著計算流體力學(CFD)的發(fā)展，CFD開始用于建筑物周圍流體流動與擴散的模擬[6-10]，這些研究并未考慮溫度層結的影響。關于溫度層結的數(shù)值模擬研究比較有限。ZHANG et al[11]使用k-ε模型TEMPEST研究了在穩(wěn)定層結條件下立方建筑物周圍大氣的流場和擴散，但其在模擬過程中并未考慮來流風切變的變化。SANTOS et al[12]使用CFD研究了不同溫度層結下單個建筑物對污染物擴散的影響，并與風洞試驗結果進行對比分析。ASHRAFI et al[13]與ORKOMI et al[14]應用CFD技術研究了不同溫度層結下煙羽的抬升與擴散規(guī)律。在我國，部分研究者采用CFD技術研究了中性層結下建筑物對周圍流場與污染物擴散的影響[15-19]，而關于溫度層結對標準體周圍流動的影響研究尚不多。

本文采用STAR-CD提供的RNGk-ε模型對不同大氣穩(wěn)定度條件下標準建筑物對附近流場的影響進行了模擬。本文建立的物理模型與YASSIN[5]應用風洞試驗研究不同溫度層結下標準建筑物對周圍流場結構影響的模型相同，該風洞試驗段長16.0 m，寬1.2 m，高1.0 m.試驗以1∶300制作模型，模型為長度(L)、寬度(W)、高度(H)均為100 mm的建筑物。

本研究圍繞不同溫度層結條件下建筑物對流場影響的數(shù)值模擬展開，采用莫寧-奧撥霍夫長度(L)對入流湍流動能(k)及其耗散廓線(ε)進行修正，建立了溫度層結數(shù)值模擬方法，摸索出溫度層結模擬技術，研究了不同溫度層結對流場結構的影響，彌補了我國在溫度層結模擬方面的欠缺。

1 數(shù)值模擬

使用STAR-CD3.26作為計算平臺，為了與風洞試驗結果比較，數(shù)值模擬的計算區(qū)域設為16.0 m×1.2 m×1.0 m(長×寬×高)，模型尺寸與YASSIN[5]風洞試驗模型相同(見圖1).網格結構采用具有良好拓撲結構的六面體網格，計算區(qū)域網格總數(shù)約為200萬，區(qū)域內最大網格尺寸為30 mm，建筑物表面最小網格尺寸為5 mm，最小網格尺寸為0.1 mm.同時建立建筑物表面最小網格尺寸分別為2 mm(細網格)、10 mm(粗網格)的網格進行網格無關性驗證，如圖2所示，三套網格在建筑物頂部處的風速模擬結果均較為相近，因此最終采用最小網格尺寸5 mm.

圖1 建筑物及其測量位置圖Fig.1 Building and survey location map

圖2 三套網格無關性驗證Fig.2 Three sets of grid independence verification

1.1 基本方程

大氣環(huán)境中用來描寫流體運動的基本方程就是Navies-Stoke方程，各種數(shù)值模擬方法都是建立在該方程基礎上?？諝饬鲃铀俣纫话悴淮?，因此可將空氣當作不可壓縮流體，不可壓縮流體運動基本控制方程如下。

A、質量守恒方程

(1)

式中：u，v，w分別為x，y，z方向風速。

B、動量守恒方程

(2)

式中：fx，fy，fz分別為x,y,z方向的運動阻力；p為環(huán)境大氣壓力；ρ為環(huán)境大氣密度。

1.2 湍流模型

近地層大氣的流動為復雜的湍流運動，目前采用雷諾(Reynolds)時均法對基本控制方程進行處理， Reynolds平均方程中引入了高階的二階脈動相關量，湍流模型就是把湍流的脈動值附加項與時均值聯(lián)系起來的一些特定關系式。而且空氣溫度變化不大，也即密度變化不大，因此認為空氣流動符合Bonssinesq假設。引入Boussinesq假定后，求解的關鍵就是如何求解湍流黏度，根據求解微分方程的不同，產生了不同的湍流模型，比如標準k-ε模型。而標準k-ε湍流模型不能準確地模擬繞流的復雜流動特征，RNGk-ε湍流模型能夠較好地模擬繞流問題[7]，所以本文選用RNGk-ε模型封閉N-S方程進行計算。RNGk-ε湍流模型方程如公式(3)-(8)，湍流模型常數(shù)值見表1.

表1 湍流模型常數(shù)Table 1 Turbulence model constants

(3)

(4)

(5)

(6)

(7)

(8)

式中：xi,xj分別代表縱向和垂直方向；ui,uj分別為xi,xj方向氣流的速度；k為湍流動能；ε為湍流耗散率；μ為空氣的運動學黏性；μt為湍流黏性；Sij為平均應變率張量；PT為湍流應力；PB為浮力應力。

1.3 邊界條件

本文入口邊界風廓線指數(shù)、不同高度處湍流強度與溫度的取值及其變化規(guī)律均與YASSIN[5]應用風洞試驗研究不同溫度層結下建筑物對周圍流場結構影響的模型相同，入口邊界條件風速沿高度方向的變化規(guī)律用指數(shù)方程描述(u=UH×(z/H)n，u為不同高度z處風速，UH為建筑物高度(H)處風速，風廓線指數(shù)n=0.25)，通過調整計算區(qū)域頂部與底部溫度實現(xiàn)不同穩(wěn)定度的模擬。計算時入口邊界條件風廓線、溫度廓線見圖3.

出口面邊界條件(outflow)為自由出口，出口邊界上流動已完全發(fā)展。計算域兩側設置為對稱邊界條件，來流風速為水平方向，可以認為速度沿切線方向的梯度為零，頂部采用滑移(Slip) 壁面條件。建筑物表面和地面設置一定的摩擦速度(u*)與粗糙度(z*)，采用無滑移(No Slip)的壁面條件并由標準壁面函數(shù)確定壁面附近流動。不同穩(wěn)定度條件下主要參數(shù)見表2.采用壓力-速度修正算法SIMPLE(Semi-Implicit Method for Pressure-Linked)求解方程，計算過程中對壓力、動量、湍流動能等作欠松弛處理。

1.4 不同溫度層結實現(xiàn)

本文通過莫寧-奧撥霍夫長度(L)描述不同溫度層結下大氣的運動，莫寧-奧撥霍夫長度主要描述近地面層湍流切應力和浮力對湍流耗散的影響，可以用下式表示：

(9)

圖3 不同溫度層結下速度廓線、溫度廓線Fig.3 Vertical distributions of mean velocity and Temperature in the simulated boundary layer under thermal stability

(10)

(11)

式中：u*為摩擦速度；T*為近地層溫度標；κ為卡門常數(shù)，0.4；qw為地表熱通量，W/m2；cp為比熱，J/(kg·K)；τw為地表剪應力。

表2 不同穩(wěn)定度條件下主要參數(shù)Table 2 Main parameters under different stability conditions

中性層結下地表熱通量為0，理查森數(shù)為0，根據莫寧-奧撥霍夫長度相似理論，中性層結下L為無窮大(∞)，Φm趨近于1(Φm是以z/L為變量的函數(shù))。在穩(wěn)定與不穩(wěn)定層結下，必須要考慮地表熱通量以及溫度梯度，根據莫寧-奧撥霍夫長度相似理論，穩(wěn)定與不穩(wěn)定層結下Φm不趨近于 1，不同高度處湍流動能廓線與湍流耗散廓線分別根據式(12)、(13)計算。

(12)

(13)

(14)

(15)

式中：Φε是以z/L為變量的函數(shù)。

2 結果分析與比較

本文主要研究不同溫度層結下標準立方體建筑物尾流區(qū)附近不同溫度層結對流動特征的影響，分別在沿建筑物中心線x/H=0.75、1.0、1.5和2.0(如圖1)四個不同位置進行流場特征的研究。

2.1 平均速度結果分析

本次數(shù)值模擬采用RNGk-ε湍流模型計算建筑物對其周圍流場的影響，在大氣環(huán)境模擬領域，通常采用歸一化速度消除不同模擬風速引起的建筑物對流場結構影響的差異。歸一化速度(u/UH)為局地縱向平均速度(u)與來流建筑物頂部縱向平均速度(UH)之比。下風向不同距離處建筑物對周圍流場影響的歸一化速度比較結果見圖4.不同層結下建筑物對其周圍流場影響的水平與垂直數(shù)值模擬結果見圖5和圖6.

由圖4可知，在回流區(qū)內速度減小，從建筑物頂開始的下風向，在近尾流區(qū)(包括回流區(qū))風速顯著減小，特別是當大氣處于穩(wěn)定狀態(tài)時，速度虧損略大于中性與不穩(wěn)定層結條件時，最大速度虧損為相應來流的30%左右，垂向速度在0.6

圖4 不同位置處不同模型歸一化速度(u/UH)隨高度的變化Fig.4 Vertical profiles of mean velocity component (u/UH) in the longitudinal direction

垂直速度廓線清楚地反映了流體的垂直方向上的流動。垂直速度在x/H=1.5、2.0處的變化較大，而在x/H=0.75、1.0處變化不明顯，這是由于受建筑物機械擾動的影響，大氣的溫度層結對流場影響不明顯，此時建筑物產生的機械湍流起主導作用。

由圖5和圖6可知，不同溫度層結下，數(shù)值模擬結果較好地模擬了建筑物附近尾流區(qū)及其回流區(qū)的信息，建筑物迎風面流場發(fā)生顯著變化，風向軸線上部分動量轉化成水平與垂直動量，建筑物后產生了回流并形成雙渦回流(馬蹄形渦)，部分流場變成垂直向下。中性層結下建筑物迎風側的靜駐點均在距地約2/3H處，空腔區(qū)的高度約為1.0H；穩(wěn)定層結下建筑物迎風側的靜駐點均在距地約1/2H處，空腔區(qū)的高度約為1.2H；不穩(wěn)定層結下建筑物迎風側的靜駐點均在距地約1/2H處，空腔區(qū)的高度約為1.1H.中性層結條件下建筑物后側再附著點約在2.3H處，穩(wěn)定層結條件下建筑物后側再附著點約在3.0H處，不穩(wěn)定層結條件下建筑物后側再附著點約在1.8H處，而不同溫度層結下數(shù)值模擬對建筑物頂部再附著點均沒有捕捉到。TOMINAGA[7]使用不同修正k-ε模型驗證立方體建筑物頂部與后側再附著點的位置分別在0.08H～0.4H與2.0H～3.0H處，本文研究結果與之相同。ZHANGet al[11]分別用數(shù)值模擬與水槽試驗進行了不同溫度層結下建筑物對流場與擴散的影響，發(fā)現(xiàn)穩(wěn)定層結條件下建筑物背風側回流區(qū)的長度大于中性層結條件下，主要由于穩(wěn)定層結條件下增加了流場的向下運動。本文研究結果與之相同。

圖5 水平剖面流場結構模擬結果Fig.5 Simulation results of flow field structure in horizontal section

圖6 垂直剖面流場結構模擬結果Fig.6 Simulation results of flow field structure in vertical section

從整個流場的特征來看，穩(wěn)定層結條件下建筑物兩側與頂部分離區(qū)域比中性層結與不穩(wěn)定層結條件下小，建筑物對水平方向的影響高于垂直方向，主要由于穩(wěn)定層結下，垂直方向大氣的運動抑制了迎風面垂直速度。穩(wěn)定層結條件下建筑物背風側回流區(qū)的長度大于中性與不穩(wěn)定層結條件下，主要由于穩(wěn)定層結條件增加了流場的向下運動，這會引起氣流產生切變，從而增加氣流的湍流動能，但是這種趨勢被浮力產生的湍流運動抑制。當處于弱穩(wěn)定條件下，浮力引起的湍流的抑制作用不足以抵消由剪切力增加引起的湍流，會造成回流區(qū)長度的增長?？傊?，數(shù)值模擬能較好地模擬受建筑物影響時氣流的位移區(qū)、尾流區(qū)、空腔區(qū)以及建筑物迎風面氣流的分離與近地面的停滯回流、建筑物頂部氣流的回流、建筑物兩側的逆于來流方向的氣流和建筑物下風方向空腔區(qū)內的馬蹄形渦，本次數(shù)值模擬在以上方面均有較好的體現(xiàn)。

2.2 湍流動能結果分析

不同溫度層結建筑物對其周圍湍流動能影響的垂直剖面數(shù)值模擬結果見圖7.

圖7 建筑物對湍流動能影響垂直剖面Fig.7 Vertical profile of the influence of buildings on turbulent kinetic energy

由圖7可知，湍流動能最大值出現(xiàn)在建筑物的迎風面，主要是由建筑物的擾動引起較高的速度梯度造成，建筑物頂部湍流動能也比較高，回流區(qū)湍流動能相對較低，但是仍然比來流方向高。建筑物頂部分離區(qū)域產生較強的湍流動能并且一直延伸到下游，在近地面湍流動能相對較小。大氣處于穩(wěn)定層結時，湍流動能的變化顯著低于中性層結與不穩(wěn)定層結時，并且在建筑物頂部，湍流動能低于中性層結與不穩(wěn)定層結時。主要原因是由于大氣處于穩(wěn)定層結時，垂直方向大氣的溫度梯度處于逆溫狀態(tài)，大氣在垂直方向上運動相對較小，從而湍流動能的變化較小，而湍流動能的變化對流場與濃度場的影響較大。當大氣處于不穩(wěn)定層結時，由于自由流和建筑物的相互作用，降低了建筑物下游靠近地面的溫度梯度，從而減少了該區(qū)域的湍流動能產生。在遠離建筑物的區(qū)域，溫度曲線逐漸恢復，接近地面的湍流動能再次增加。在不穩(wěn)定條件下，湍流動能較大，氣流在建筑物頂和側壁上的流動分離遠小于在中性和穩(wěn)定條件下的流動分離，使得建筑物后的再循環(huán)區(qū)域的長度也較小。

3 結論

通過對大氣處于不同溫度層結時標準建筑物對周圍流場結構影響的數(shù)值模擬研究，結果表明：

1) 當大氣處于穩(wěn)定層結時，在建筑物尾流區(qū)范圍內近地層風速較低，建筑物尾流區(qū)范圍內最大速度虧損大于中性與不穩(wěn)定層結時，尾流區(qū)中回流區(qū)的長度更長。中性條件下模擬建筑物后側再附著點約在2.0H處，穩(wěn)定條件下模擬建筑物后側再附著點約在3.0H處，不穩(wěn)定條件下模擬建筑物后側再附著點約在1.8H處。

2) 隨著下風向距離的增加，熱穩(wěn)定性對湍流動能的影響逐漸顯現(xiàn)。在建筑物附近，由于建筑物擾動與熱穩(wěn)定性的共同影響，湍流動能相對較大。

本文圍繞不同溫度層結條件下建筑物對流場結構影響的CFD數(shù)值模擬展開，采用L對入口湍流動能與耗散進行修正，建立了溫度層結CFD數(shù)值模擬方法，摸索出溫度層結模擬技術，研究了不同溫度層結對流場結構的影響。結果表明，CFD數(shù)值模擬技術對其周圍流場研究是一種較為有效的途徑，特別是當建筑物結構比較復雜時，數(shù)值模擬技術更加能顯現(xiàn)出其對流場和污染物擴散模擬的優(yōu)點。