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

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

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

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

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

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

1 數值模擬

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

(1)

(2)

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

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

(3)

(4)

(5)

(6)

式中，Pk為切應力相關的機械k的產生項，m2/s2；Pb為浮力相關的k的產生項，m2/s2；σk、σε、σh為湍流普朗特數；Cε1、Cε2、Cε3、Cμ為湍流模型經驗常數。表1給出了計算中模型常數。

表1 湍流模型常數Tab.1 Turbulence model constants

2 模擬描述

2.1 幾何描述

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

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

2.2 計算域和網格

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

圖3 網格無關性驗證Fig.3 The grid independence verification

2.3 邊界條件

(1)入流和出流邊界

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

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

穩(wěn)定條件：

(12)

(13)

不穩(wěn)定條件：

(14)

(15)

(16)

(2)對稱與壁面邊界

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

(17)

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

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

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

3 模型有效性分析

3.1 風洞實驗描述

Yassin[4]風洞實驗主要研究不同溫度層結建筑物對流場及其污染物擴散規(guī)律的影響，該風洞實驗段長16.0 m，寬1.2 m，高1.0 m。試驗以1∶300制作模型，模型長度(L)、寬度(W)、高度(H)均為100 mm的建筑物。為了與風洞實驗結果比較，數值模擬的計算區(qū)域設為16.0 m×1.2 m×1.0 m(長×寬×高)，風洞實驗模型見圖4，風洞實驗主要參數見表3，不同溫度層結下溫度層結下速度廓線、溫度廓線見圖5。

3.2 模擬設置

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

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

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

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

3.3 流場有效性分析

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

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

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

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

4 建筑物群對流場影響

4.1 模型的建立

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

4.2 結果分析

(1)縱向速度，u/UH

圖7顯示了歸一化流向速度u/uref在建筑物矩陣中心面y/H=0上的分布。不同穩(wěn)定度條件下，u/UH均出現了明顯的垂直梯度。街谷內部出現了大范圍的負向速度區(qū)域，且在街谷底部負向速度值較大，u/UH<-0.2的區(qū)域出現在街谷底部區(qū)域(z/H<0.5)。這說明不同穩(wěn)定度條件下，街谷底部出現了強烈的回流。不同穩(wěn)定度條件下，第一個街谷中u/UH<-0.2的區(qū)域明顯大于后續(xù)街谷，第二排建筑物之后各個街谷中u/UH的分布大體上是相似的，且相對于第一個街谷，不同穩(wěn)定度條件下u/UH<-0.2的區(qū)域均發(fā)生減小。這說明第一個街谷中渦旋的強度是非常強烈的，因為氣流通過第一排建筑物時產生的強烈的機械擾動對渦旋起主導作用。第二個及后續(xù)的街谷中，u/UH<-0.2的區(qū)域在不穩(wěn)定條件下最大，中性條件下次之，穩(wěn)定條件下最小，該區(qū)域在穩(wěn)定條件下減小最為明顯。這是因為第二排建筑物后，溫度層結對街谷渦旋的影響變得明顯，穩(wěn)定的大氣條件明顯減弱了渦旋強度。不同穩(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范圍內u/UH值均為負，這證明了街谷內1/2建筑高度以下范圍內以回流為主。z/H<0.5范圍內，水平位置x/H=-1過渡至1，|u/UH|值逐漸減小，表明回流速度從街谷內上游至下游部分發(fā)生減小。在x/H=-1和1處，不穩(wěn)定條件下的|u/UH|值略大于中性和穩(wěn)定條件的，這表明不穩(wěn)定條件下回流速度較大。街谷中，不同穩(wěn)定度條件下，u/UH隨高度增加逐漸趨于入流狀態(tài)，而穩(wěn)定條件下u/UH恢復至入流水平的高度略大于中性和不穩(wěn)定條件。這表明穩(wěn)定條件下建筑矩陣對較高處風速的衰減影響較大。

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

圖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)流場結構

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

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

5 結論

通過對大氣處于不同溫度層結時建筑物群對周圍流場結構影響的數值模擬研究，本文首先應用Yassin[4]風洞實驗結果對數值模擬有效性進行分析比較，其次研究不同溫度層結時理想建筑物群對流場結構的影響。結果表明：

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

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

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

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