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

(天津大學 水利工程仿真與安全國家重點實驗室，天津 300072)

1 引 言

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

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

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

圖1 防風網減風抑塵機理

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

2 數(shù)學模型建立

2.1 基于浸入邊界法的風場控制方程

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

(1)

(2)

(3)

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

(4)

(5)

(6)

(7)

2.2 防風網PM-IBM模型

防風網PM-IBM模型核心是根據(jù)網后壓降確定風場控制方程中的附加體積力大小。防風網兩側壓降與來流雷諾數(shù)Re相關，即

Re=ud/υ

(8)

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

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

(9)

(10)

(11)

綜上，當來流風垂直于防風網時，不考慮法向力源，則IB點流向力源及IB點附近網格的附加體積力計算形式為

(12)

(13)

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

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

(2) DEFINE_ADJUST。將歐拉網格流場信息插值到IB點并計算IB點附加體積力。

(3) DEFINE_SOURCE。利用式(13)計算網格附加體積力并嵌入動量方程中。

通過3個模塊，IBM方法能夠在流體方程中實現(xiàn)多孔介質，形成PM-IBM模型。利用該模型，不需要再引入另外的多孔介質模型及防風網邊界條件，該邊界條件由PM-IBM的散點(Immersed Boundary Points)代替。該方法同傳統(tǒng)多孔介質相比，無需繪制物面邊界，對于不同形狀邊界只需導入相應的邊界散點坐標，對復雜幾何邊界條件的適應能力強。

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

3 模型驗證

選用張寧[20]的風洞實驗進行模型驗證。試驗段長為6.75 m，寬為0.72 m，高為0.6 m，防風網高H0為0.03 m，開孔率為38.5%。實驗布置如 圖3 所示。

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

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

(14)

(15)

(16)

固壁邊界ux=uz=0

(17)

圖3 實驗布置模型[20]

圖4 計算域和網格劃分

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

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

(18)

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

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

(19)

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

(20)

式中z1和z2為高度，u1和u2為高度z1和z2處平均測量風速，κ為馮卡門常數(shù)，z0為空氣動力學地表粗糙度。數(shù)值模擬入口速度分布與實驗穩(wěn)定風速對比如圖5所示，數(shù)值模擬風速分布精度高，近底區(qū)誤差在1%以內，說明數(shù)值風洞能夠為防風網提供與物理風洞相近的來流邊界。

圖5 入口斷面風速分布驗證

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

4 計算結果及分析

根據(jù)模型驗證結果，基于PM-IBM模型采用Realizablek-ε模型研究開孔率及防風網高度影響。

4.1 開孔率的影響

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

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

圖7 不同開孔率風速分布

圖8 剩余風速沿程變化

4.2 防風網高度的影響

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

圖9 堆場布置及邊界條件

堆場附近風場結構如圖10所示，有網工況網與迎風面間形成低風速區(qū)，減風效果明顯；無網工況迎風面及堆頂風速較大，抑塵效益低。受堆場自身遮蔽效應影響，來流風被迫抬升，背風面均存在風速較低回流區(qū)。

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

圖10 近堆場風場結構(H=H0)

圖11 不同網高風速垂向分布

5 結 論

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

(1) PM-IBM模型能夠有效模擬防風網的阻風效應，其中Realizablek-ε模型模擬最為準確。

(2) 防風網最優(yōu)開孔率為30%～35%，網后最優(yōu)庇護區(qū)在高度<0.8H,沿程2H～8H之間。

(3) 當網高小于堆場高1.4倍時，堆頂風速隨網高增大近似線性降低；當網高大于堆場高1.4倍時，風速降低幅度趨于平緩。