溫度層結(jié)下建筑物群周圍流場影響的數(shù)值模擬研究

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

(1.太原科技大學(xué)，太原 030024；2.中國輻射防護(hù)研究院，太原 030006；3.山西省環(huán)境保護(hù)技術(shù)評估中心,太原 030024)

對于核設(shè)施廠址(如核電站等)，其存在一些陣列式的建筑物群(如圖1)，由于其復(fù)雜下墊面與大氣穩(wěn)定度的影響，使其周圍流場變得非常復(fù)雜，流線彎曲、速度劇烈地不連續(xù)及湍流分布不均勻。其中，大氣穩(wěn)定度對氣載放射性物質(zhì)擴(kuò)散的影響尤為重要。因此，為更加真實(shí)模擬氣載放射性物質(zhì)擴(kuò)散情況以及對環(huán)境的影響，需要研究不同溫度層結(jié)下復(fù)雜建筑物群對流動與污染物擴(kuò)散的影響。

圖1 某核設(shè)施廠址建筑物群Fig.1 Buildings at a nuclear facility site

對于這類特殊工程廠址對環(huán)境的影響起初使用風(fēng)洞模擬手段進(jìn)行[1-4]，隨著計(jì)算流體力學(xué)(CFD)的發(fā)展，CFD模擬技術(shù)已逐漸用于預(yù)測各種建筑物周圍流動特性以及建筑物周圍大氣污染物的擴(kuò)散規(guī)律[5-8]，但這些研究均是建立在中性層結(jié)的基礎(chǔ)上進(jìn)行的。關(guān)于溫度層結(jié)的數(shù)值模擬研究比較有限。Zhang等[9]使用k-ε模型在中性和穩(wěn)定層結(jié)條件下研究了標(biāo)準(zhǔn)體近場污染物的流動與擴(kuò)散規(guī)律，但其使用的是均一來流的入口條件，未考慮來流風(fēng)切變的影響。Santos等[10]使用k-ε模型，以立方體為研究對象，研究了不同溫度層結(jié)對污染物擴(kuò)散的影響，并與風(fēng)洞實(shí)驗(yàn)進(jìn)行比較。Ashrafi等[11]與Orkomi等[12]應(yīng)用CFD技術(shù)研究了不同溫度層結(jié)煙羽的抬升與擴(kuò)散規(guī)律。在我國，部分學(xué)者應(yīng)用CFD數(shù)值模擬技術(shù)在中性層結(jié)下對建筑物近場污染物流動與擴(kuò)散規(guī)律進(jìn)行研究[13-17]，而關(guān)于不同溫度層結(jié)下建筑物群近場的流動模擬研究尚不多。

本文采用STAR-CD提供的RNGk-ε模型對不同大氣穩(wěn)定度條件下建筑物對附近流場的影響進(jìn)行了模擬。首先應(yīng)用CFD技術(shù)對Yassin[4]風(fēng)洞實(shí)驗(yàn)進(jìn)行模擬，并將模擬結(jié)果與風(fēng)洞實(shí)驗(yàn)結(jié)果進(jìn)行對比分析，檢驗(yàn)其對溫度層結(jié)的模擬能力，建立溫度層結(jié)CFD數(shù)值模擬方法。并參照某核設(shè)施廠址，選取8×7建筑物群，對不同溫度層結(jié)(不穩(wěn)定、中性和穩(wěn)定)條件下對規(guī)則建筑物矩陣中流動的影響規(guī)律進(jìn)行初步研究，為評價(jià)核設(shè)施廠址復(fù)雜下墊面對流動影響及其近場氣載放射性污染物擴(kuò)散影響奠定基礎(chǔ)。

1 數(shù)值模擬

基于Reynolds averaged Navier-Stokes (RANS)方法的三維CFD模型STAR-CD V3.26被設(shè)計(jì)用于平坦和復(fù)雜地形中大氣流動和污染物擴(kuò)散模擬。其使用修正標(biāo)準(zhǔn)k-ε湍流模型求解穩(wěn)態(tài)流場，并通過有限體積法(FVM)離散流動和湍流的控制方程。本模擬使用Boussinesq近似簡化計(jì)算，即假設(shè)計(jì)算域中氣體為不可壓縮流體，且密度定常。連續(xù)性和動量守恒方程如下：

(1)

(2)

式中，xi(x1=x,x2=y,x3=z)分別為縱向、橫向和垂直方向；Ui為沿xi的平均速度分量，m/s；P為壓力，Pa；μt為湍流渦黏性；μ為運(yùn)動學(xué)黏性；ρ為密度，kg/m3；SB=gρβ(Tamb-T)，為浮力源項(xiàng)；T為溫度，℃；Tamb為環(huán)境溫度，℃；g為重力加速度，kg/s2；β=1/T，為熱膨脹系數(shù)，1/℃。

湍流動能(k)以及耗散率(ε)傳輸方程如下：

(3)

(4)

(5)

(6)

式中，Pk為切應(yīng)力相關(guān)的機(jī)械k的產(chǎn)生項(xiàng)，m2/s2；Pb為浮力相關(guān)的k的產(chǎn)生項(xiàng)，m2/s2；σk、σε、σh為湍流普朗特?cái)?shù)；Cε1、Cε2、Cε3、Cμ為湍流模型經(jīng)驗(yàn)常數(shù)。表1給出了計(jì)算中模型常數(shù)。

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

2 模擬描述

2.1 幾何描述

本文研究不同溫度層結(jié)條件下規(guī)則建筑物矩陣中對流動的影響。模擬使用的三維矩形障礙物高度(H)為10 m，其大小為H×3H×H(長×寬×高)。研究布局是由56個障礙物組成的矩陣，其規(guī)模為：8個障礙物分布于流向上，7個障礙物分布于橫向上。矩陣流向長度(Lx)和橫向?qū)挾?Ly)分別為36H和33H。建筑間隔在流向和橫向上分別為4H和2H。坐標(biāo)系原點(diǎn)(x/H=0和y/H=0)位于矩陣水平橫截面(z/H=0)的中心[見圖2(a)]。矩陣第一排建筑物迎風(fēng)面位于x/H=-18，最后一排建筑物背風(fēng)面位于x/H=18。建筑物矩陣幾何示意圖詳見圖2。

圖2 規(guī)則建筑群幾何布局Fig.2 Geometrical configuration of the regular buildings

2.2 計(jì)算域和網(wǎng)格

模擬中計(jì)算域在流向和橫向長度均為60H，計(jì)算域高度為25H，底部中心位于坐標(biāo)系原點(diǎn)。使用212×208×58個笛卡爾網(wǎng)格對計(jì)算域進(jìn)行離散化。網(wǎng)格在垂直方向、流向和橫向上分布是不均勻的。在矩陣中，水平方向上網(wǎng)格尺寸為0.2H，水平網(wǎng)格分辨率在矩陣上風(fēng)向和下風(fēng)向區(qū)域逐漸減小。在2H高度以下網(wǎng)格保持均勻最小的垂直間距，最小垂直網(wǎng)格尺寸為0.1H，垂直網(wǎng)格間距在2H高度以上逐漸增大。使用一套精細(xì)的網(wǎng)格280×268×62進(jìn)行網(wǎng)格獨(dú)立性分析，圖3比較了矩陣中心處使用精細(xì)網(wǎng)格和粗糙網(wǎng)格得到的流向速度垂直分布。散點(diǎn)圖中數(shù)據(jù)點(diǎn)大致分布在1∶1線上，網(wǎng)格獨(dú)立性分析結(jié)果表明兩套網(wǎng)格中流向速度差異是很小的。由于使用精細(xì)網(wǎng)格并不會造成研究結(jié)論的改變，所以本文選用粗糙網(wǎng)格進(jìn)行模擬。

圖3 網(wǎng)格無關(guān)性驗(yàn)證Fig.3 The grid independence verification

2.3 邊界條件

(1)入流和出流邊界

模擬中計(jì)算域側(cè)面邊界類型需要根據(jù)入流風(fēng)向確定，即上游邊界設(shè)為入流邊界，下游邊界設(shè)為出流邊界。本研究考慮垂直于建筑物矩陣的入流風(fēng)向(平行于x軸)并基于莫寧-奧布霍夫相似理論來描述表面層。入口處平均風(fēng)速和溫度垂直廓線由下列公式定義：

(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，絕熱直減率；κ為馮卡門常數(shù)，0.41；Qh為地面感熱通量，W/m2；Cp為定壓比熱，kJ/(kg·K)；Ψ,ψ分別為穩(wěn)定度修正函數(shù)，在中性條件下Ψ和ψ均等于0。

穩(wěn)定條件：

(12)

(13)

不穩(wěn)定條件：

(14)

(15)

(16)

(2)對稱與壁面邊界

計(jì)算域頂部及兩側(cè)面(平行于x軸)應(yīng)用對稱邊界條件，這些邊界上所有變量正向梯度均等于0。底部邊界定義為無滑移壁面條件并設(shè)置粗糙度z0等于0.02 m，該壁面上速度分量設(shè)為0并使用標(biāo)準(zhǔn)壁面函數(shù)計(jì)算其它變量。壁面動量由如下對數(shù)律法則描述：

(17)

式中，E為壁面粗糙度函數(shù)；u+=Up/u*，是無量綱速度，其中Up是平行于壁面的流體速度；y+=ρu*y/μ，是壁面至網(wǎng)格中心的無量綱距離；y為網(wǎng)格中心到壁面的距離。

不同模擬邊界層中參考高度處(zr=H)風(fēng)速(UH)和溫度(TH)分別設(shè)為2 m/s和298 K。表2給出了不同類型模擬邊界層的微氣象參數(shù)。

表2 模擬邊界層主要微氣象參數(shù)Tab.2 Micro-meteorology parameters of the simulated boundary layers

3 模型有效性分析

3.1 風(fēng)洞實(shí)驗(yàn)描述

Yassin[4]風(fēng)洞實(shí)驗(yàn)主要研究不同溫度層結(jié)建筑物對流場及其污染物擴(kuò)散規(guī)律的影響，該風(fēng)洞實(shí)驗(yàn)段長16.0 m，寬1.2 m，高1.0 m。試驗(yàn)以1∶300制作模型，模型長度(L)、寬度(W)、高度(H)均為100 mm的建筑物。為了與風(fēng)洞實(shí)驗(yàn)結(jié)果比較，數(shù)值模擬的計(jì)算區(qū)域設(shè)為16.0 m×1.2 m×1.0 m(長×寬×高)，風(fēng)洞實(shí)驗(yàn)?zāi)Ｐ鸵妶D4，風(fēng)洞實(shí)驗(yàn)主要參數(shù)見表3，不同溫度層結(jié)下溫度層結(jié)下速度廓線、溫度廓線見圖5。

3.2 模擬設(shè)置

為了與風(fēng)洞實(shí)驗(yàn)結(jié)果比較，數(shù)值模擬的計(jì)算區(qū)域與風(fēng)洞實(shí)驗(yàn)相同，數(shù)值模擬的計(jì)算區(qū)域設(shè)為16.0 m×1.2 m×1.0 m(長×寬×高)，模型尺寸與Yassin[4]風(fēng)洞實(shí)驗(yàn)?zāi)Ｐ拖嗤?見圖4)，計(jì)算域中網(wǎng)格總數(shù)為180萬。網(wǎng)格結(jié)構(gòu)采用具有良好拓?fù)浣Y(jié)構(gòu)的六面體網(wǎng)格，計(jì)算區(qū)域內(nèi)最大網(wǎng)格尺寸為30 mm，建筑物表面網(wǎng)格尺寸為5 mm。計(jì)算域中上風(fēng)向設(shè)為入流邊界，下風(fēng)邊界則為出流邊界。計(jì)算域頂部應(yīng)用對稱邊界條件，底部使用無滑移壁面條件。入流速度、溫度和湍流廓線的描述均與之前(2.3節(jié))保持一致，不同溫度層結(jié)下，擬合入流速度，溫度廓線見圖5，并依據(jù)表3給出的數(shù)據(jù)設(shè)置輸入?yún)?shù)。

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

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

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

3.3 流場有效性分析

本文分別在沿建筑物中心線x/H=0.75、1.0、1.5和2.0(如圖4)四個不同位置進(jìn)行流場特征的有效性分析，在大氣環(huán)境模擬領(lǐng)域，通常采用歸一化速度消除不同模擬風(fēng)速引起的建筑物對流場結(jié)構(gòu)的影響差異。歸一化速度(u/UH)為局地縱向平均速度(u)與來流建筑物頂部縱向平均速度(UH)之比。下風(fēng)向不同距離處建筑物對周圍流場影響的歸一化速度比較結(jié)果見圖6。

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

由圖6可知，不同溫度層結(jié)下本文CFD數(shù)值模擬結(jié)果與Yassin風(fēng)洞實(shí)驗(yàn)結(jié)果吻合較好，變化趨勢基本一致。在建筑物背風(fēng)向近地面速度顯著減小，特別是當(dāng)大氣處于穩(wěn)定狀態(tài)時。隨著下風(fēng)距離的增大，氣流混合逐漸均勻，風(fēng)速逐漸恢復(fù)到來流狀態(tài)，但是大氣處于穩(wěn)定層結(jié)時，風(fēng)速恢復(fù)相對較慢，不穩(wěn)定層結(jié)風(fēng)速恢復(fù)相對較快。

當(dāng)大氣處于穩(wěn)定狀態(tài)時，在建筑物高度以下(z/H<1.0)近地面層速度減小略高于中性與不穩(wěn)定狀態(tài)，特別是在x/H=1.5、2.0處，而在x/H=0.75、1.0處，由于受建筑物機(jī)械擾動的影響，大氣的溫度層結(jié)對流場影響不明顯。

4 建筑物群對流場影響

4.1 模型的建立

通過標(biāo)準(zhǔn)體對模擬方法的驗(yàn)證，將本文建立的溫度層結(jié)數(shù)值模擬技術(shù)應(yīng)用于核設(shè)施廠址理想建筑物群流場結(jié)構(gòu)的模擬。本文使用56個高度(H)為10 m的矩形建筑物組成的矩陣來模擬核設(shè)施廠址理想建筑物群，每個障礙物大小 (x×y×z)為：H×3H×H。研究區(qū)域由8×7個建筑物組成，縱向(x軸)上每排分布有8個，橫向上每列分布有7個。建筑物矩陣縱向距離和橫向距離分別為36H和33H，建筑物之間縱向和橫向間隔分別為4H和2H。笛卡爾坐標(biāo)系原點(diǎn)位于障礙物矩陣中心，即x/H=0和y/H=0兩條線在z/H=0平面上的交匯處(見圖2a)。第一排建筑物迎風(fēng)面位于x/H=-18，最后一排建筑物背風(fēng)面位于x/H=18。該建筑物矩陣的幾何布局見圖2。

4.2 結(jié)果分析

(1)縱向速度，u/UH

圖7顯示了歸一化流向速度u/uref在建筑物矩陣中心面y/H=0上的分布。不同穩(wěn)定度條件下，u/UH均出現(xiàn)了明顯的垂直梯度。街谷內(nèi)部出現(xiàn)了大范圍的負(fù)向速度區(qū)域，且在街谷底部負(fù)向速度值較大，u/UH<-0.2的區(qū)域出現(xiàn)在街谷底部區(qū)域(z/H<0.5)。這說明不同穩(wěn)定度條件下，街谷底部出現(xiàn)了強(qiáng)烈的回流。不同穩(wěn)定度條件下，第一個街谷中u/UH<-0.2的區(qū)域明顯大于后續(xù)街谷，第二排建筑物之后各個街谷中u/UH的分布大體上是相似的，且相對于第一個街谷，不同穩(wěn)定度條件下u/UH<-0.2的區(qū)域均發(fā)生減小。這說明第一個街谷中渦旋的強(qiáng)度是非常強(qiáng)烈的，因?yàn)闅饬魍ㄟ^第一排建筑物時產(chǎn)生的強(qiáng)烈的機(jī)械擾動對渦旋起主導(dǎo)作用。第二個及后續(xù)的街谷中，u/UH<-0.2的區(qū)域在不穩(wěn)定條件下最大，中性條件下次之，穩(wěn)定條件下最小，該區(qū)域在穩(wěn)定條件下減小最為明顯。這是因?yàn)榈诙沤ㄖ锖?，溫度層結(jié)對街谷渦旋的影響變得明顯，穩(wěn)定的大氣條件明顯減弱了渦旋強(qiáng)度。不同穩(wěn)定度條件下街谷底部上游部分|u/UH|值大于下游部分。

圖7 不同穩(wěn)定度條件下，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范圍內(nèi)u/UH值均為負(fù)，這證明了街谷內(nèi)1/2建筑高度以下范圍內(nèi)以回流為主。z/H<0.5范圍內(nèi)，水平位置x/H=-1過渡至1，|u/UH|值逐漸減小，表明回流速度從街谷內(nèi)上游至下游部分發(fā)生減小。在x/H=-1和1處，不穩(wěn)定條件下的|u/UH|值略大于中性和穩(wěn)定條件的，這表明不穩(wěn)定條件下回流速度較大。街谷中，不同穩(wěn)定度條件下，u/UH隨高度增加逐漸趨于入流狀態(tài)，而穩(wěn)定條件下u/UH恢復(fù)至入流水平的高度略大于中性和不穩(wěn)定條件。這表明穩(wěn)定條件下建筑矩陣對較高處風(fēng)速的衰減影響較大。

圖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上的分布。在第一排建筑迎風(fēng)角附近，強(qiáng)烈的風(fēng)切變作用造成了w/UH強(qiáng)烈的垂直分布。街谷中，w/UH在中部和迎風(fēng)側(cè)附近為負(fù)值，在背風(fēng)側(cè)附近為正值。這表明氣流在街谷背風(fēng)側(cè)和迎風(fēng)側(cè)附近的垂直運(yùn)動分別是向上和向下的。街谷內(nèi)部，背風(fēng)側(cè)向上的運(yùn)動屬于上游渦旋的一部分，而迎風(fēng)側(cè)向下的運(yùn)動則屬于下游渦旋。從圖中可以看出，不同穩(wěn)定度條件下，背風(fēng)側(cè)附近|w/UH|值要略大于迎風(fēng)側(cè)，說明了背風(fēng)側(cè)向上的氣流運(yùn)動要強(qiáng)于迎風(fēng)側(cè)向上的運(yùn)動。這反映出街谷內(nèi)部上游渦旋流動要強(qiáng)于下游渦旋。街谷背風(fēng)側(cè)附近，w/UH值在不穩(wěn)定條件下要略大于中性和穩(wěn)定條件，且w/UH>0.2的區(qū)域在不穩(wěn)定條件下是最大的，在穩(wěn)定條件下w/UH均小于0.15。這是因?yàn)榉€(wěn)定的大氣條件抑制了垂直方向的氣流運(yùn)動，而在不穩(wěn)定條件下，這種運(yùn)動是更強(qiáng)烈的。

圖9 不同穩(wěn)定度條件下，y/H=0平面上歸一化速度w/UH的等值線云圖Fig.9 Contours of w/UH on the plane y/H=0 under different stability conditions

圖10 不同穩(wěn)定度條件下，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)流場結(jié)構(gòu)

圖12顯示了前3個街谷中y/H=0平面上速度場流線圖。不同穩(wěn)定度條件下，所有街谷內(nèi)部可以觀察到兩個渦旋的出現(xiàn)，這流動模式即OKE[18]提出的尾跡干擾流動。不同穩(wěn)定度條件下，街谷1中的渦旋結(jié)構(gòu)略不同于后續(xù)街谷(街谷2和3)，特別是在穩(wěn)定條件下，這種差異更為明顯。穩(wěn)定條件下街谷1中上游渦旋規(guī)模較小且并不完整。在街谷2和3中流動趨于穩(wěn)定，不同穩(wěn)定度條件下街谷中上游渦旋的規(guī)模明顯大于下游渦旋。相較于街谷1，街谷2和3中上游渦旋規(guī)模略微減小，上游渦旋中心距街谷背風(fēng)側(cè)的水平距離略也較短，而下游渦旋規(guī)模出現(xiàn)略微增大。不穩(wěn)定條件下，上游渦旋的中心高度和尺度均大于中性和穩(wěn)定條件。下游渦旋在穩(wěn)定條件下也略小于不穩(wěn)定和中性條件。

圖12 前3個街谷中y/H=0垂直剖面的流場結(jié)構(gòu)Fig.12 Streamlines of the velocity field on the plane y/H=0

5 結(jié)論

通過對大氣處于不同溫度層結(jié)時建筑物群對周圍流場結(jié)構(gòu)影響的數(shù)值模擬研究，本文首先應(yīng)用Yassin[4]風(fēng)洞實(shí)驗(yàn)結(jié)果對數(shù)值模擬有效性進(jìn)行分析比較，其次研究不同溫度層結(jié)時理想建筑物群對流場結(jié)構(gòu)的影響。結(jié)果表明：

(1) CFD 模擬結(jié)果與Yassin[4]風(fēng)洞實(shí)驗(yàn)結(jié)果能較好地吻合，并且研究發(fā)現(xiàn)穩(wěn)定層結(jié)建筑尾流區(qū)范圍內(nèi)速度和湍流動能減小。

(2) 不同溫度層結(jié)下，街谷內(nèi)流向速度受建筑物影響均出現(xiàn)很大程度的減弱。在街谷底部均有強(qiáng)烈的回流出現(xiàn)，并且不穩(wěn)定條件下街谷底部回流速度大于中性和穩(wěn)定條件。在第2街谷之后，不同溫度層結(jié)下街谷內(nèi)u/UH分布開始達(dá)到平衡。不穩(wěn)定條件下，建筑物對縱向通道內(nèi)u/UH的影響較小，縱向速度分布在與入流狀態(tài)基本一致。而中性和穩(wěn)定條件下，通道內(nèi)u/UH隨距離的增大逐漸出現(xiàn)明顯衰減，并且z/H<2.0范圍內(nèi)u/UH分布逐漸趨于線性。穩(wěn)定條件下建筑物矩陣對縱向速度的最大影響高度大于不穩(wěn)定和中性條件。

(3) 不同溫度層結(jié)下街谷內(nèi)上游和下游部分均有渦旋出現(xiàn)。并且溫度層結(jié)對渦旋大小和中心位置有著顯著的影響。不穩(wěn)定條件下，街谷渦旋強(qiáng)度大于中性和穩(wěn)定條件。穩(wěn)定條件下由于大氣層結(jié)的抑制，街谷渦旋強(qiáng)度是較弱的。

不同溫度層結(jié)下，由于建筑物群內(nèi)部與環(huán)境大氣溫度存在溫度差，從而影響了流場結(jié)構(gòu)與湍流特征，本文通過與風(fēng)洞實(shí)驗(yàn)結(jié)果的對比，建立了溫度層結(jié)數(shù)值模擬技術(shù)，為下一步應(yīng)用于評價(jià)不同溫度層結(jié)下核設(shè)施建筑物群近場氣載放射性污染物的流動與擴(kuò)散奠定了基礎(chǔ)。