許 棟， 張博曦， 及春寧， 楊海滔

(天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

1 引 言

防風(fēng)網(wǎng)是一種多孔墻式結(jié)構(gòu)，能夠降低網(wǎng)后局部風(fēng)速，形成對(duì)顆粒揚(yáng)塵的有效抑制，在港口堆場(chǎng)環(huán)境治理中廣泛應(yīng)用[1]。來(lái)流風(fēng)在防風(fēng)網(wǎng)上部形成繞流，與主流風(fēng)速匯集形成高速風(fēng)；下部由于阻風(fēng)效應(yīng)，網(wǎng)前高強(qiáng)度大尺度漩渦轉(zhuǎn)化為低強(qiáng)度小尺度漩渦，大幅降低網(wǎng)后起塵量[2-4]，減風(fēng)抑塵機(jī)理如圖1所示。防風(fēng)網(wǎng)風(fēng)場(chǎng)研究方法主要有風(fēng)洞實(shí)驗(yàn)及數(shù)值模擬，風(fēng)洞實(shí)驗(yàn)可信度高，但實(shí)驗(yàn)周期長(zhǎng)，設(shè)備制定復(fù)雜；隨計(jì)算機(jī)技術(shù)發(fā)展，基于計(jì)算流體力學(xué)方法建立的計(jì)算機(jī)數(shù)值風(fēng)洞得到廣泛應(yīng)用，王婷婷等[5]建立數(shù)值風(fēng)洞研究數(shù)值邊界精度，表明數(shù)值風(fēng)洞能夠有效模擬大氣邊界層；Scotta等[6]通過(guò)數(shù)值風(fēng)洞研究橋梁所受風(fēng)荷載，表明了數(shù)值風(fēng)洞的可靠性。

防風(fēng)網(wǎng)多孔結(jié)構(gòu)數(shù)值模擬主要有3D實(shí)體建模及多孔介質(zhì)PM模擬。3D實(shí)體建模方面，陳廷國(guó)等[7]針對(duì)防風(fēng)網(wǎng)截面形式，通過(guò)貼體網(wǎng)格構(gòu)建防風(fēng)網(wǎng)模型，表明蝶形網(wǎng)的抑塵性能優(yōu)于平面網(wǎng)；李建隆等[8]針對(duì)防風(fēng)網(wǎng)開(kāi)孔形狀，通過(guò)混合貼體網(wǎng)格建立防風(fēng)網(wǎng)模型，表明圓形開(kāi)孔網(wǎng)的抑塵效果最佳；3D實(shí)體建模對(duì)開(kāi)孔結(jié)構(gòu)細(xì)節(jié)模擬更準(zhǔn)確，但網(wǎng)格劃分復(fù)雜，網(wǎng)格質(zhì)量難以保證，建模耗時(shí)長(zhǎng)。多孔介質(zhì)方面，考慮防風(fēng)網(wǎng)阻流效應(yīng)實(shí)質(zhì)為壓降變化，引入多孔介質(zhì)力學(xué)參數(shù)(滲透系數(shù)及慣性系數(shù))實(shí)現(xiàn)壓降變化，力學(xué)參數(shù)由網(wǎng)板風(fēng)洞實(shí)驗(yàn)或刻畫(huà)精細(xì)的CFD模擬取得[9]。Maruyama[10]結(jié)合多孔介質(zhì)及大渦模擬對(duì)比風(fēng)洞實(shí)驗(yàn)網(wǎng)后流場(chǎng)分布；許棟等[11]通過(guò)多孔介質(zhì)模型探討數(shù)值邊界對(duì)計(jì)算精度的影響，表明多孔介質(zhì)能夠準(zhǔn)確模擬防風(fēng)網(wǎng)阻流效應(yīng)；Teitel[12]采用Forchheimer壓降方程模擬柔性開(kāi)孔網(wǎng)，表明多孔介質(zhì)力學(xué)參數(shù)受多因素影響。盡管多孔介質(zhì)建模簡(jiǎn)便，但受防風(fēng)網(wǎng)截面影響，前后壓降關(guān)系復(fù)雜，針對(duì)防風(fēng)網(wǎng)模擬的關(guān)鍵力學(xué)參數(shù)設(shè)定困難[11,12]；且多孔介質(zhì)依賴(lài)網(wǎng)格邊界，不同防風(fēng)網(wǎng)形式及工程布置需要重新劃分網(wǎng)格，建模效率低。

傳統(tǒng)浸入邊界法IBM將固壁邊界離散為拉格朗日坐標(biāo)下的IB點(diǎn)，在IB點(diǎn)處引入動(dòng)量方程附加體積力，實(shí)現(xiàn)不可穿透固壁邊界[9,13]。該方法能夠快速構(gòu)建物面邊界，但要求防風(fēng)網(wǎng)附近網(wǎng)格達(dá)孔徑級(jí)別，網(wǎng)格數(shù)目大，如許棟等[9]采用浸入邊界法對(duì)防風(fēng)網(wǎng)模擬的計(jì)算網(wǎng)格總數(shù)約為2000萬(wàn)。

圖1 防風(fēng)網(wǎng)減風(fēng)抑塵機(jī)理

本文結(jié)合多孔介質(zhì)及浸入邊界法，將壓降以源項(xiàng)形式加入流體動(dòng)量方程中，建立多孔介質(zhì)-浸入邊界模型(PM-IBM模型，下同)，實(shí)現(xiàn)邊界的部分穿透性，結(jié)合通用計(jì)算流體力學(xué)求解器，為防風(fēng)網(wǎng)截面選型和布置優(yōu)化等工程應(yīng)用提供高效建模工具；并探討不同湍流模型下PM-IBM模型的計(jì)算精度，分析網(wǎng)高及開(kāi)孔率對(duì)減風(fēng)抑塵效果的影響。

2 數(shù)學(xué)模型建立

2.1 基于浸入邊界法的風(fēng)場(chǎng)控制方程

近地層空氣受風(fēng)速和溫度等因素產(chǎn)生的密度相對(duì)改變率小于1%[14]，可視為不可壓縮粘性流體。浸入邊界法在不可壓縮粘性流體N-S方程中引入附加體積力，流體控制方程為

(1)

(2)

(3)

式中u為速度矢量,p為壓強(qiáng),t為時(shí)間,Re為雷諾數(shù),f為附加體積力矢量,為梯度算子,rhs(x,tn)為對(duì)流項(xiàng)、粘性項(xiàng)及壓力梯度在tn時(shí)刻之和,V(Xi,tn + 1)為可移動(dòng)固體邊界期望速度。本文采用四點(diǎn)正則插值函數(shù)，形式為

(4)

(5)

(6)

(7)

2.2 防風(fēng)網(wǎng)PM-IBM模型

防風(fēng)網(wǎng)PM-IBM模型核心是根據(jù)網(wǎng)后壓降確定風(fēng)場(chǎng)控制方程中的附加體積力大小。防風(fēng)網(wǎng)兩側(cè)壓降與來(lái)流雷諾數(shù)Re相關(guān)，即

Re=ud/υ

(8)

式中u為風(fēng)速，υ為空氣運(yùn)動(dòng)粘滯系數(shù)，d為柔性網(wǎng)線直徑或剛性網(wǎng)孔間距。

Re<10時(shí)，防風(fēng)網(wǎng)風(fēng)阻效應(yīng)可按多孔介質(zhì)滲流處理，適用達(dá)西定律，故網(wǎng)兩側(cè)壓降見(jiàn)式(9)。Re>10時(shí)，受孔隙慣性效應(yīng)影響，形成非達(dá)西流動(dòng)，壓降由Forchheimer方程[8]求解，見(jiàn)式(10)；當(dāng)Re>450時(shí)，K0趨近于0[18]；當(dāng)Re>250時(shí)，慣性系數(shù)β僅為開(kāi)孔率α的函數(shù)，A1(Re)為常數(shù)，約為 0.52[19]。故Re>450時(shí)，壓降由式(10)簡(jiǎn)化為式(11)。

(9)

(10)

(11)

綜上，當(dāng)來(lái)流風(fēng)垂直于防風(fēng)網(wǎng)時(shí)，不考慮法向力源，則IB點(diǎn)流向力源及IB點(diǎn)附近網(wǎng)格的附加體積力計(jì)算形式為

(12)

(13)

2.3 PM-IBM模型軟件實(shí)現(xiàn)

FLUENT是應(yīng)用最廣的通用CFD軟件之一。模型采用FLUENT用戶自定義函數(shù)(UDF)通過(guò)3個(gè)功能模塊實(shí)現(xiàn)PM-IBM模型，如圖2所示，各模塊功能如下。

(2) DEFINE_ADJUST。將歐拉網(wǎng)格流場(chǎng)信息插值到IB點(diǎn)并計(jì)算IB點(diǎn)附加體積力。

(3) DEFINE_SOURCE。利用式(13)計(jì)算網(wǎng)格附加體積力并嵌入動(dòng)量方程中。

通過(guò)3個(gè)模塊，IBM方法能夠在流體方程中實(shí)現(xiàn)多孔介質(zhì)，形成PM-IBM模型。利用該模型，不需要再引入另外的多孔介質(zhì)模型及防風(fēng)網(wǎng)邊界條件，該邊界條件由PM-IBM的散點(diǎn)(Immersed Boundary Points)代替。該方法同傳統(tǒng)多孔介質(zhì)相比，無(wú)需繪制物面邊界，對(duì)于不同形狀邊界只需導(dǎo)入相應(yīng)的邊界散點(diǎn)坐標(biāo)，對(duì)復(fù)雜幾何邊界條件的適應(yīng)能力強(qiáng)。

圖2 基于FLUENT的PM-IBM模型實(shí)現(xiàn)

3 模型驗(yàn)證

選用張寧[20]的風(fēng)洞實(shí)驗(yàn)進(jìn)行模型驗(yàn)證。試驗(yàn)段長(zhǎng)為6.75 m，寬為0.72 m，高為0.6 m，防風(fēng)網(wǎng)高H0為0.03 m，開(kāi)孔率為38.5%。實(shí)驗(yàn)布置如 圖3 所示。

數(shù)值風(fēng)洞涉及大量網(wǎng)格計(jì)算，覆蓋整個(gè)實(shí)驗(yàn)區(qū)域，計(jì)算量大，考慮區(qū)別于實(shí)驗(yàn)，數(shù)值風(fēng)洞能夠通過(guò)入流邊界和底部粗糙度快速形成穩(wěn)定大氣邊界層[5,11]，故選取網(wǎng)前10H0、網(wǎng)后30H0、高10H0的范圍作為計(jì)算域，能避免邊界對(duì)計(jì)算結(jié)果的影響[3,21]；防風(fēng)網(wǎng)采用PM-IBM模型，計(jì)算域左側(cè)為速度入口邊界(Velocity inlet)，右側(cè)為自由出流邊界(Outflow)，上部為對(duì)稱(chēng)邊界(Symmetry)，下部為不可滑移固壁邊界(Wall)，各邊界條件數(shù)學(xué)表達(dá)見(jiàn)式(14～17)。網(wǎng)格采用正交結(jié)構(gòu)網(wǎng)格，近網(wǎng)區(qū)局部加密，總網(wǎng)格數(shù)為69000，計(jì)算域及網(wǎng)格劃分如圖4所示。

入口邊界ux=f1(x,z)，uz=f2(x,z)

(14)

(15)

(16)

固壁邊界ux=uz=0

(17)

圖3 實(shí)驗(yàn)布置模型[20]

圖4 計(jì)算域和網(wǎng)格劃分

入流邊界速度分布采用對(duì)數(shù)律擬合為

uz=(u*/κ)ln(z/z0)

(18)

式中uz為高度z處的平均風(fēng)速，u*為地表摩擦風(fēng)速，κ=0.41為馮卡門(mén)常數(shù)，z0為空氣動(dòng)力學(xué)地表粗糙度。根據(jù)實(shí)驗(yàn)數(shù)據(jù)[20]，u*和z0由已知不同高度的風(fēng)速值求得，

z0=exp[(u1lnz2-u2lnz1)/(u1-u2)]

(19)

u*=u1κ/[ln(z1/z0)]

(20)

式中z1和z2為高度，u1和u2為高度z1和z2處平均測(cè)量風(fēng)速，κ為馮卡門(mén)常數(shù)，z0為空氣動(dòng)力學(xué)地表粗糙度。數(shù)值模擬入口速度分布與實(shí)驗(yàn)穩(wěn)定風(fēng)速對(duì)比如圖5所示，數(shù)值模擬風(fēng)速分布精度高，近底區(qū)誤差在1%以?xún)?nèi)，說(shuō)明數(shù)值風(fēng)洞能夠?yàn)榉里L(fēng)網(wǎng)提供與物理風(fēng)洞相近的來(lái)流邊界。

圖5 入口斷面風(fēng)速分布驗(yàn)證

基于PM-IBM模型，分別采用Standard模型、RNG模型和Realizable模型對(duì)防風(fēng)網(wǎng)進(jìn)行數(shù)值模擬，同傳統(tǒng)多孔介質(zhì)模型(流體采用Realizable模型)及實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，網(wǎng)后距離為2H0，4H0，6H0和8H0處的水平風(fēng)速如圖6所示。網(wǎng)后上部高速條帶區(qū)除Standardk-ε模型外均與風(fēng)洞實(shí)驗(yàn)較吻合，誤差在2%以?xún)?nèi)，Santiago等[22]研究表明，RNGk-ε模型和Realizablek-ε模型能夠準(zhǔn)確模擬網(wǎng)后上部高速條帶區(qū)，與本文結(jié)論一致。近網(wǎng)區(qū)底部數(shù)值模擬同實(shí)驗(yàn)有所差異，其中Realizablek-ε模型斷面風(fēng)速分布同實(shí)驗(yàn)誤差最小，偏差范圍2.6%～4.1%；RNGk-ε模型近網(wǎng)近底區(qū)存在回流。PM-IBM模型同傳統(tǒng)多孔介質(zhì)模型對(duì)比，近網(wǎng)區(qū)域流速基本一致，網(wǎng)后6H～8H范圍PM-IBM模型近底(0H～0.5H)流速略?xún)?yōu)于多孔介質(zhì)模型。PM-IBM模型總體精度能夠有效模擬防風(fēng)網(wǎng)阻風(fēng)效應(yīng)。

4 計(jì)算結(jié)果及分析

根據(jù)模型驗(yàn)證結(jié)果，基于PM-IBM模型采用Realizablek-ε模型研究開(kāi)孔率及防風(fēng)網(wǎng)高度影響。

4.1 開(kāi)孔率的影響

開(kāi)孔率α是影響防風(fēng)網(wǎng)遮蔽效應(yīng)的關(guān)鍵因素，本文對(duì)20%～45%開(kāi)孔率防風(fēng)網(wǎng)空氣動(dòng)力學(xué)特性進(jìn)行數(shù)值模擬，計(jì)算域如圖4所示，網(wǎng)后流速分布如圖7所示，可以看出，網(wǎng)后存在明顯風(fēng)速較低區(qū)域；20%～30%開(kāi)孔率網(wǎng)后5～7倍網(wǎng)高存在回流現(xiàn)象，易引起顆粒起揚(yáng)；開(kāi)孔率>35%時(shí)，回流區(qū)基本消失。此外，隨開(kāi)孔率增大，防風(fēng)網(wǎng)阻風(fēng)效應(yīng)減弱，風(fēng)速增大。

定義無(wú)量綱剩余風(fēng)速u(mài)z/Uz為網(wǎng)后風(fēng)速與入流斷面等高處的風(fēng)速比，不同高度剩余風(fēng)速沿程分布見(jiàn)圖8?？梢钥闯觯叨葄≥0.6H0范圍，剩余風(fēng)速沿程呈先減后增趨勢(shì)；高度z≤0.4H0范圍，受局部回流影響，20%～25%開(kāi)孔率剩余風(fēng)速沿程于5倍網(wǎng)高附近存在局部增加現(xiàn)象，抑塵效益降低。最小剩余風(fēng)速隨開(kāi)孔率增大而增大，20%開(kāi)孔率減風(fēng)效果最優(yōu)，該結(jié)果同Lee等[23]以雙頓粒子測(cè)速技術(shù)所得20%開(kāi)孔率防風(fēng)網(wǎng)對(duì)網(wǎng)后順流向速度折減最大的結(jié)論一致。綜合剩余風(fēng)速及網(wǎng)后回流因素，建議最優(yōu)防風(fēng)網(wǎng)開(kāi)孔率為30%～35%，最優(yōu)庇護(hù)區(qū)在z≤0.8H0和x= 2H0～8H0范圍內(nèi)，此區(qū)域剩余風(fēng)速不超過(guò)0.3，抑塵效益高。

圖7 不同開(kāi)孔率風(fēng)速分布

圖8 剩余風(fēng)速沿程變化

4.2 防風(fēng)網(wǎng)高度的影響

網(wǎng)高是影響庇護(hù)區(qū)的重要結(jié)構(gòu)因素之一，本文設(shè)置網(wǎng)高H=H0，1.2H0，1.4H0，1.6H0和1.8H0和未設(shè)防風(fēng)網(wǎng)情況，其中初始堆場(chǎng)高度H0=0.03 m，開(kāi)孔率為35%。堆場(chǎng)前端距防風(fēng)網(wǎng)2H0，計(jì)算域及邊界條件設(shè)置如圖9所示，堆場(chǎng)采用不可穿透浸入邊界施加，防風(fēng)網(wǎng)采用PM-IBM模型。

圖9 堆場(chǎng)布置及邊界條件

堆場(chǎng)附近風(fēng)場(chǎng)結(jié)構(gòu)如圖10所示，有網(wǎng)工況網(wǎng)與迎風(fēng)面間形成低風(fēng)速區(qū)，減風(fēng)效果明顯；無(wú)網(wǎng)工況迎風(fēng)面及堆頂風(fēng)速較大，抑塵效益低。受堆場(chǎng)自身遮蔽效應(yīng)影響，來(lái)流風(fēng)被迫抬升，背風(fēng)面均存在風(fēng)速較低回流區(qū)。

不同網(wǎng)高堆場(chǎng)表面速度垂向分布如圖11所示，堆場(chǎng)表面高度2H范圍內(nèi)，有網(wǎng)風(fēng)速明顯小于無(wú)網(wǎng)工況，迎風(fēng)面風(fēng)速最大減小約65%，堆頂及背風(fēng)面最大減小約55%。對(duì)堆頂，H<1.4H0時(shí)，堆場(chǎng)表面風(fēng)速隨網(wǎng)高增加線性減?。籋≥ 1.4H0時(shí)，網(wǎng)高與風(fēng)速變化的線性關(guān)系趨于平緩；對(duì)背風(fēng)面，網(wǎng)高增加對(duì)表面風(fēng)速影響趨于0；對(duì)迎風(fēng)面，表面風(fēng)速隨網(wǎng)高增加逐漸減小?？紤]防風(fēng)網(wǎng)建設(shè)成本，建議防風(fēng)網(wǎng)高度為堆場(chǎng)高度的1.4倍。

圖10 近堆場(chǎng)風(fēng)場(chǎng)結(jié)構(gòu)(H=H0)

圖11 不同網(wǎng)高風(fēng)速垂向分布

5 結(jié) 論

本文首次將多孔介質(zhì)與浸入邊界法相結(jié)合，提出了能準(zhǔn)確、快速模擬防風(fēng)網(wǎng)邊界的PM-IBM模型，避免了對(duì)開(kāi)孔邊界的直接描述，使網(wǎng)格數(shù)量大幅降低，在保證計(jì)算精度和計(jì)算效率前提下，實(shí)現(xiàn)結(jié)構(gòu)化網(wǎng)格下防風(fēng)網(wǎng)的高效數(shù)值建模，結(jié)果如下。

(1) PM-IBM模型能夠有效模擬防風(fēng)網(wǎng)的阻風(fēng)效應(yīng)，其中Realizablek-ε模型模擬最為準(zhǔn)確。

(2) 防風(fēng)網(wǎng)最優(yōu)開(kāi)孔率為30%～35%，網(wǎng)后最優(yōu)庇護(hù)區(qū)在高度<0.8H,沿程2H～8H之間。

(3) 當(dāng)網(wǎng)高小于堆場(chǎng)高1.4倍時(shí)，堆頂風(fēng)速隨網(wǎng)高增大近似線性降低；當(dāng)網(wǎng)高大于堆場(chǎng)高1.4倍時(shí)，風(fēng)速降低幅度趨于平緩。