王冰清，付海明

(東華大學 a. 環(huán)境科學與工程學院；b. 國家環(huán)境保護紡織工業(yè)污染防治工程技術(shù)中心，上海 201620)

?

樹冠微觀尺度流動阻力特性分析

王冰清a, b，付海明a， b

(東華大學 a. 環(huán)境科學與工程學院；b. 國家環(huán)境保護紡織工業(yè)污染防治工程技術(shù)中心，上海 201620)

摘要：采用可視化Basic程序設計語言(VBA)建立虛擬的三維樹冠模型，并使用計算流體動力學(CFD)模擬的方式，研究了不同孔隙率和風速下，樹冠內(nèi)部的流場特征和流動阻力. 用光學孔隙率及體積孔隙率來表征樹冠復雜的形態(tài)結(jié)構(gòu)，并討論了樹冠形態(tài)結(jié)構(gòu)參數(shù)對樹冠流動阻力的影響規(guī)律. 結(jié)果表明，樹冠的阻力系數(shù)與光學孔隙率滿足近似二次多項式關(guān)系，通過對數(shù)值模擬計算結(jié)果回歸分析，給出樹冠流動阻力系數(shù)與光學孔隙率的關(guān)系式，在真實條件下樹冠阻力系數(shù)的試驗數(shù)據(jù)對比的基礎上，提供了一個對真實樹冠阻力系數(shù)進行預測的依據(jù).

關(guān)鍵詞：樹冠； 阻力系數(shù)； 孔隙率； 微觀模型

森林是陸地生態(tài)系統(tǒng)中重要的因素，在全球的生態(tài)平衡中有著重要的作用.樹木具有保護水土、防風固沙、涵養(yǎng)水源、調(diào)節(jié)氣候、減少污染物等功能.自然或人為作用會使土壤流失、土壤結(jié)構(gòu)發(fā)生改變.防風林的栽種有助于防止土壤流失、減少風蝕，最終使風速降低.充分樹冠密度下樹木的種植或自然植被覆蓋的維護會降低風蝕的可能性.文獻[1]研究表明，在農(nóng)田中防風林使農(nóng)作物免受風災，其創(chuàng)造了一個提高農(nóng)作物生產(chǎn)率的有利微環(huán)境.若使防風林的效益最大化，還需要深入理解風和障礙物相互之間的物理作用.

盡管研究者們在防風林和大范圍孤立障礙物背風處的氣流測量和特征描述上做了很多努力，但是對單個植被內(nèi)氣流繞流相互作用而產(chǎn)生的阻力系數(shù)方面關(guān)注較少.正是由于此相互作用，在微觀層面上，觀察背風處流域是可靠的.通常，研究分為兩方面：(1)剪切應力分區(qū)，這里植被的作用是保護固定表面的土壤[2-10]；(2)防風林對流場和微氣候的影響[11-14].在第二方面中，大部分的工作用來描繪二維籬笆流場，因此這提高了對文獻[15]中二維人造籬笆的遮擋效應和空氣動力學的理解.文獻[16]總結(jié)了二維籬笆孔隙率與阻力系數(shù)的曲線關(guān)系的研究.實心籬笆的阻力系數(shù)最高，并隨著孔隙率增大而持續(xù)減小.

已有的研究[1]對于防風林、孤立樹木和灌木等三維多孔介質(zhì)繞流流動的理解很不完整.由于模型中代替數(shù)據(jù)的使用造成了這部分結(jié)論的缺失，文獻[7-8]的研究使用了文獻[16]中的實心粗糙度物體阻力系數(shù)來代表自然多孔的植物.風速減小、壓力擾動的原因正是由于滲透率和阻力，但是三維多孔介質(zhì)的原因是未知的[15].防風林阻力系數(shù)的測試大部分是通過障礙物周圍的風場動量損失的計算完成的[12-13，15，17].但也有其他嘗試，例如，通過縮放植被比例的風洞試驗進行阻力系數(shù)測量[1].風洞研究也測量了闊葉[18]、云杉樹枝[19]及整樹[20]的阻力.

因此，研究樹冠層結(jié)構(gòu)參數(shù)與樹冠阻力系數(shù)的關(guān)系對于預測防沙護林的作用是至關(guān)重要的.本文采用可視化Basic程序設計語言 (VBA)建立虛擬的樹冠模型，并采用計算流體動力學(CFD)模擬計算得出虛擬樹冠阻力，通過對數(shù)值模擬計算結(jié)果回歸分析，給出樹冠流動阻力系數(shù)與光學孔隙率的關(guān)系式，并進行真實條件下樹冠阻力系數(shù)的試驗測試，將試驗結(jié)果與模擬結(jié)果進行比較，為真實樹冠的阻力系數(shù)預測提供計算依據(jù).

1物理模型

本文真實樹木與模擬模型中虛擬樹木的比例為10∶1.計算流域高為3h，長為5h，寬為防護林平均行距3.5 m，其中h為樹冠的平均高度.因此，如圖1所示為本文模擬的計算域，設樹冠高度為0.25 m，流域高為0.75 m(3h)，長為1.25 m(5h)，寬為0.35 m. 流域高度為3倍樹冠高度，故可忽略流域高度對流場求解結(jié)果的影響.為保證入口邊界和出口邊界處流動的充分發(fā)展，設定樹冠到流體入口邊界距離為0.625 m，即2.5倍樹冠高，樹冠到流體出口邊界距離為0.625 m，亦為2.5倍樹冠高.流域被劃分為8 385 633個非均勻四面體網(wǎng)格，并對樹冠內(nèi)部作局部網(wǎng)格加密處理.邊界條件：流域的底部設置為粗糙壁面，粗糙度取0.005 m,并假設壁面無明顯的障礙物和植被.除了入口和出口以外，計算流域的其他壁面包括樹枝和樹葉的表面均設為無滑移壁面邊界條件.入口設置為速度入口邊界條件，出口設置為相對壓力為0的壓力出口.事實上，樹枝的表面是粗糙的.但是之前的研究中沒有關(guān)于樹枝表面粗糙度的數(shù)據(jù)，文獻[21]研究得到的樹枝表面粗糙度為0.002 m，將其作為參考用于代表全尺寸的樹木表面粗糙度.流域中的三維微觀樹冠模型是由VBA生成的，其樹冠中的樹枝葉子是由計算機編程產(chǎn)生的，其外形及特征與真實樹冠十分相似，每一個虛擬樹冠都具有不同的隨機分布結(jié)構(gòu)，可近似認為本文的虛擬樹冠是“真實”的，如圖2所示.

圖1　樹冠模型模擬的計算域和邊界條件Fig.1　Computational domain and boundary conditions for simulation with model trees

圖2　微觀模型中的虛擬樹冠Fig.2　Picture of the virtual tree used in the microcosmic model

2風洞試驗

為了驗證三維微觀樹冠模型的可行性，選取槐樹樹枝進行風洞試驗以驗證數(shù)值模擬結(jié)果.風洞試驗裝置和測點位置(截面A)如圖3所示.風洞由有機玻璃材料制成.

1—風機；2—變頻裝置；3—軟接頭；4—孔板流量計；5—擴壓器；6—蜂窩柵極；7—測試段圖3　風動試驗裝置Fig.3　Diagram of the wind tunnel for the airflow experiment

將選取的樹枝置于測試段(圓形風筒)，每種樹枝設定5種不同光學孔隙率，共計進行20組試驗.光學孔隙率通過控制試驗樹枝的數(shù)量以及葉片數(shù)來實現(xiàn).光學孔隙率取值范圍為0～1.根據(jù)光學孔隙率定義[22]可知，在二維水平面內(nèi)，測試段7底面積與樹枝葉片面積重合部分之差與測試段7底面積之比即為光學孔隙率，也可稱其為疏透度，其定義：在垂直于某一水平方向上，樹冠邊緣垂直面上的透光孔隙的投影面積與該垂直面上的樹體投影總面積之比.試驗時，通過無極變頻裝置2(0～50 Hz內(nèi)無極調(diào)節(jié))和蜂窩柵極6均勻向測試段輸送速度1.0～ 15.4 m/s的連續(xù)風.并在測試段7的B、C兩點采用數(shù)字壓差計測量氣流通過樹枝前后的壓力損失，B、C兩點分別在測試段圓周上均勻取6個測試點以測出A、B兩點平均靜壓差.測試段7尾部通向大氣，每變換一次風速等待片刻，待風速均勻穩(wěn)定后再進行測試. 風速測量可以通過孔板流量計間接測得，也可以直接在測試段7使用熱線風速儀測量.為獲得更高的測量精度，本試驗采用后者.

施加于植被上的全部阻力與穿過樹冠內(nèi)部和周圍的氣流流動有關(guān)，大部分由孔隙率控制.選取5組樹冠模型進行二維和三維孔隙率的測量.二維光學孔隙率采用數(shù)字化掃描方法，首先取得清晰的樹木照片，用數(shù)字化掃描儀把照片轉(zhuǎn)化為數(shù)字化圖像文件，利用樹冠與背景之間的灰度差異，使用自編程計算軟件計算樹葉面積，求得透光孔隙面積，從而得到光學孔隙率或疏透度.三維體積孔隙率的測量使用簡單的排水體積法，假設樹是實心無孔隙的,然后樹被拆卸，每個枝干放入一裝滿水的量筒中，將每個枝干體積相加得到孔隙率.在計算模型中最重要的一個參數(shù)是葉面積密度(LAD).模型中虛擬樹木的LAD獲得方法與光學孔隙率相似.

物體表面的阻力系數(shù)(Cd)可由式(1)確定

(1)

式中：u為風速；ρ為流體密度；p為壓力.變量u、ρ可直接測量，而p必須直接測量或者通過Cd確定.為了簡化幾何形狀，可以通過確定的阻力曲線得到Cd.

3結(jié)果與討論

3.1樹冠流場特性

采用在動量和湍流輸運方程中附加源/匯項的方法模擬樹冠內(nèi)的繞流效應，但此方法要明確地解釋樹冠內(nèi)各向異性的空間結(jié)構(gòu)所導致的氣流復雜性是不可能的.事實上，樹冠內(nèi)繞流的相互作用會在樹冠枝干尾跡區(qū)域產(chǎn)生大量小范圍風速減小，使用先前的傳統(tǒng)方法無法實現(xiàn)這些小范圍的風速減小.迄今為止，文獻[21]通過在模型中引入有明確樹枝結(jié)構(gòu)，研究了無葉樹冠的流場，預測效果良好.本文的創(chuàng)新之處在于，三維微觀樹冠模型在樹枝周圍增加了“真實葉子”，能夠深入分析樹冠結(jié)構(gòu)對樹冠內(nèi)部及周圍區(qū)域繞流流動的影響.圖4(a)顯示了具有樹葉的樹冠層內(nèi)部流場速度分布圖，而圖4(b)顯示了樹冠層內(nèi)部流場分布的細節(jié)圖，這是先前的傳統(tǒng)方法無法實現(xiàn)的.由圖4可知，樹冠內(nèi)的風速整體是減小的，而在樹枝、樹葉周圍和背風處的風速卻是增加的.

(a) 整體速度云圖，h= 2.25 m (LAD值=25.89 m-1)

(b) A處速度云圖細節(jié)圖

3.2孔隙率

本文使用的三維微觀樹冠模型的孔隙率測試結(jié)果如表1所示. 由表1可知，光學孔隙率變化范圍為0.15～0.64，而體積孔隙率會大一些，其變化范圍為0.93～0.96.

表1　孔隙率測試結(jié)果

圖5給出了虛擬樹冠的光學孔隙率和體積孔隙率的關(guān)系.由圖5可知，光學孔隙率是體積孔隙率的冪函數(shù).在文獻[22]的研究基礎上，將表1孔隙率分析結(jié)果進行擬合，則光學孔隙率與體積孔隙率的關(guān)系如式(2)所示.

(2)

式中:Po為光學孔隙率；Pv為體積孔隙率.

圖5　微觀模型的體積孔隙率與光學孔隙率關(guān)系Fig.5　The relationship between volumetric porosity and optical porosity for the microcosmic model

式(2)函數(shù)與數(shù)據(jù)的相關(guān)系數(shù)為0.805，滿足實心物體模型孔隙率為0的客觀要求，若無物體，光學孔隙率和體積孔隙率均為1.由式(2)可知，文獻[22]總結(jié)的光學孔隙率與體積孔隙率的關(guān)系曲線的斜率明顯小于本文得到的兩者關(guān)系曲線的斜率，這主要是由于研究的植物結(jié)構(gòu)特性及葉子的長寬比不同.文獻[22]使用的是松樹，其體積增長速度快于面積，因此體積孔隙率的增長比例大于光學孔隙率.而本文使用的是虛擬槐樹，由于槐樹葉面積較大，其面積增長速度快于體積，因此體積孔隙率的增長比例小于光學孔隙率.不僅由于葉面結(jié)構(gòu)的不同，同種植被擁有不同深度和寬度比時(本文深寬比為1∶1)，冪指數(shù)都可能會改變.

3.3阻力研究

本文使用的風速范圍對應的雷諾數(shù)(Re)范圍為0.4×105～1.7×105，Re與Cd的關(guān)系如圖6所示. 文獻[16]測得實心物體的表面阻力系數(shù)為0.6.由圖6可知，在本文所考慮的所有Re范圍內(nèi)，有孔隙的介質(zhì)均比實心物體具有更大的阻力系數(shù).隨著風速的增大，Re值增大，但相同孔隙率下的每種樹冠模型的阻力系數(shù)均減小.

(a) u=3 m/s

(b) u=10 m/s

different wind speed

在低Re值和高Re值下，阻力系數(shù)Cd與光學孔隙率Po的關(guān)系曲線如圖7所示.

圖7　光學孔隙率與阻力系數(shù)關(guān)系圖Fig.7　The relationship between Po and Cd

由圖7可知，隨著Po從0增加到1，Cd值先增大，達到最大值后減小,直到孔隙率為1時減到0.由圖7還可以明顯發(fā)現(xiàn)，三維阻力系數(shù)隨Po的變化明顯不同于文獻[16]總結(jié)的二維籬笆孔隙率與阻力系數(shù)的曲線關(guān)系，其二維風阻的阻力系數(shù)隨著光學孔隙率增大而逐漸減小.將圖7中的數(shù)據(jù)進行擬合,可得到如下關(guān)系式

1.087×10-6Re+0.843，R2=0.95

(3)

由式(3)可知，Re對Cd的影響較小，若忽略Re對關(guān)系式的影響，擬合得到

Cd=-4.46Po2+3.72Po+0.843

(4)

將不同Re值下阻力系數(shù)隨光學孔隙率的關(guān)系曲線與其他文獻進行對比,結(jié)果如圖8所示.由圖8可知，隨著光學孔隙率從0增加到1，Cd均是先增大，達到最大值后減小,直到孔隙率為1時減到0.由圖8中還可以明顯發(fā)現(xiàn)，文獻[22]總結(jié)的阻力系數(shù)隨光學孔隙率的關(guān)系曲線中，其Cd最大值對應的光學孔隙率小于本文得到的光學孔隙率，這與圖5中光學孔隙率和體積孔隙率的關(guān)系曲線相對應，其主要原因是研究的植物結(jié)構(gòu)特性及葉子的長寬比不同.隨著孔隙率的增加，低Re值和高Re值下的阻力系數(shù)曲線開始匯聚，這表明隨著孔隙率的增加，阻力系數(shù)Cd開始獨立于雷諾數(shù)Re.

圖8　光學孔隙率與阻力系數(shù)關(guān)系圖與其他文獻對比Fig.8　The relationship between Po and Cd compared with other literature

根據(jù)文獻[23]提出的樹冠動量源項Su：

Su=-Cd×ρAfu2

(5)

式中：Cd為阻力系數(shù)；ρ為空氣密度；Af為葉面積指數(shù)；u為空氣平均流速.

樹冠阻力Δp的表達式為

Δp=Cd×ρu2/2

(6)

將式(4)代入式(6)可得：

3.72Po+0.843)×ρu2/2

(7)

本文采用槐樹進行試驗，分別取5種光學孔隙率，則可得到5條不同光學孔隙率的壓力損失-速度曲線.分別將這5條曲線所對應的光學孔隙率以及速度代入式(7)得到5條模擬曲線.但模擬公式得到的壓差偏大，故需乘以修正系數(shù)k，即：

3.72Po+0.843)×ρu2/2

(8)

由此得出的模擬數(shù)據(jù)與試驗數(shù)據(jù)進行對比，結(jié)果如圖9所示.表2列出了5種光學孔隙率下試驗數(shù)據(jù)與經(jīng)過修正的模擬結(jié)果的平均相對誤差.

圖9　試驗數(shù)據(jù)與經(jīng)過修正的擬合結(jié)果對比Fig.9　Comparison of the modified simulation results and measured data

Po0.190.270.440.550.70平均相對誤差/%3.902.126.043.082.78

3.72Po+0.843)

(9)

式中：k為試驗修正系數(shù)，校正試驗中所產(chǎn)生的誤差，其數(shù)值為0.62～1.25.光學孔隙率小的修正系數(shù)選取較小數(shù)值，光學孔隙率大的修正系數(shù)選取較大的數(shù)值.

4結(jié)語

為深入分析三維微觀樹冠內(nèi)部結(jié)構(gòu)參數(shù)對樹冠繞流流動阻力系數(shù)的影響，采用VBA程序構(gòu)建復雜的三維虛擬樹冠模型.采用數(shù)值模擬和試驗測試相結(jié)合的方法研究不同孔隙率和風速下樹冠的內(nèi)部流場特征和流動阻力.用樹冠葉面積密度(LAD)、光學孔隙率(Po)以及體積孔隙率(Pv)來表征樹冠復雜的形態(tài)結(jié)構(gòu)，提出樹冠結(jié)構(gòu)參數(shù)間的關(guān)聯(lián)關(guān)系，并討論樹冠形態(tài)結(jié)構(gòu)參數(shù)對樹冠流動阻力的影響規(guī)律,并得到下述結(jié)論.

(1) 通過在流場中放置一個由VBA生成的三維虛擬樹冠模型來預測樹木枝干及“真實葉子”對樹冠內(nèi)部及周圍區(qū)域繞流流動的影響，能夠捕捉到空氣繞流流動特性的細節(jié)，深入分析樹冠結(jié)構(gòu)對樹冠內(nèi)部及周圍區(qū)域繞流流動的影響.

(2) 孔隙率的分析結(jié)果表明了光學空隙率和體積孔隙率之間的關(guān)系，盡管經(jīng)過試驗檢驗這兩者間的關(guān)系可能是真實的，但不同樹冠結(jié)構(gòu)的光學孔隙率和體積孔隙率之間的關(guān)系有待進一步深入研究.

(3) 樹冠的阻力系數(shù)與光學孔隙率滿足近似二次多項式關(guān)系，故采用光學孔隙率分析樹冠阻力系數(shù)，并且該多項式可以直接應用到防風林的單一樹冠風阻作用的預測中.

參考文獻

[1] HEISLER G M, DEWALLE D R. Effects of windbreak structure on wind flow [J]. Agriculture, Ecosystems and Environment, 1988,22(23): 41-69.

[2] MARSHALL J K. Drag measurements in roughness arrays of varying density and distribution [J].Agricultural Meteorology, 1971, 8(3): 269-292.

[3] GILLETTE D A, STOCKTON P H. The effect of nonerodible particles on wind erosion of erodible surfaces[J]. Journal of Geophysical Research, 1989, 94(10): 12885-12893.

[4] MUSICK H B, GILLETTE D A. Field evaluation of relationships between a vegetation structural parameter and sheltering against wind erosion[J]. Land Degradation & Rehabilitation, 1990, 2(2): 87-94.

[5] STOCKTON P H, GILLETTE D A. Field measurement of the sheltering effect of vegetation on erodible land surfaces[J]. Land Degradation and Rehabilitation, 1990, 2(2): 77-85.

[6] IVERSEN J D, WANG W P, RASMUSSEN K R,et al. Roughness element effect on local and universal saltation transport[J]. Acta Mechanica, 1991, 2: 65-75.

[7] RAUPACH M R. Drag and drag partition on rough surfaces [J]. Boundary-Layer Meteorology, 1992, 60(3): 375-395.

[8] RAUPACH R A, GILLETTE D A, LEYS J F. The effect of roughness elements on wind erosion thresholds[J]. Journal of Geophysical Research, 1993, 98(2): 3023-3029.

[9] NICKLING W G, MCKENNA N C. Development of deflation lag surfaces[J]. Sedimentology, 1995, 42(3): 403-414.

[10] WOLFE S A, NICKLING W G. Shear stress partitioning in sparsley vegetated desert canopies[J]. Earth Surface and Landforms, 1996, 21(7): 607-619.

[11] HAGEN L J, SKIDMORE E L. Turbulent velocity fluctuations and vertical flow as affected by windbreak porosity[J]. Transactions of the ASAE, 1971, 14(4): 634-637.

[12] SEGINER I, SAGI R. Drag on a windbreak in two-dimensional flow[J]. Agricultural Meteorology, 1971/1972, 9: 323-333.

[13] SEGINER I. Windbreak drag calculated from the horizontal velocity profile[J]. Boundary-Layer Meteorology, 1972, 3(1): 87-97.

[14] WILSON J D. On the choice of a windbreak porosity profile[J]. Boundary Layer Meteorology, 1987, 38(1): 37-49.

[15] WANG H, TAKLE E S. On three-dimensionality of shelterbelt structure and its influences on shelter effects[J]. Boundary-Layer Meteorology, 1996, 79(1): 83-105.

[16] TAYLOR P A. Flow and transport in the natural environment: Advances and applications [M]. New York : Springer-Verlag, 1988: 270-292.

[17] WANG S, TAKLE E S. A numerical simulation of boundary-layer flows near shelterbelts[J]. Boundary-Layer Meteorology, 1995, 75(1): 141-173.

[18] VOGEL S. Drag and reconfiguration of broad leaves in high winds[J]. Journal of Experimental Botany, 1989, 40(8): 941-948.

[19] GRANT R. The influence of the physical attributes of a spruce shoot on momentum transfer[J]. Agricultural and Forest Meteorology, 1985, 24: 7-18.

[20] MAYHEAD G J. Some drag coefficients for British trees dericed from wind tunnel studies[J]. Agricultural Meteorology, 1973, 12: 123-130.

[21] WIERINGA J. Updating the davenport roughness classification[J]. Wind Engineering and Industrial Aerodynamics, 1992, 41(1), 357-368.

[22] GRANT P F, NICKLING W G. Direct field measurement of wind drag on vegetation for application to windbreak design and modeling[J]. Land Degradation & Development, 1998, 9(1):57-66.[23] SANZ C. A note onκ-εmodelling of vegetation canopy air-flows [J]. Boundary-Layer Meteorology, 2003, 108(1): 191-197.

文章編號：1671-0444(2016)03-0426-06

收稿日期：2015-05-12

基金項目：國家自然科學基金資助項目(41371445)

作者簡介：王冰清(1990—)，女，遼寧撫順人，碩士研究生，研究方向為空氣過濾凈化與空氣品質(zhì). E-mail：wbq1827@163.com 付海明(聯(lián)系人), 男，高級工程師，E-mail: fhm@dhu.edu.cn

中圖分類號：X 513; S 731.2

文獻標志碼：A

Flow Resistance Characteristics Analysis of Micro-scale Canopy

WANGBing-qinga, b,FUHai-minga, b

(a. School of Environmental Science and Engineering; b. State Environmental Protection Engineering Center for Pollution Treatment and Control in Textile Industry, Donghua University, Shanghai 201620, China)

Abstract:A three-dimensional(3D) architecture of canopy was produced by VBA (visual basic application) and CFD (computational fluid dynamics) was used to investigate the flow characteristics and flow resistance through 3D microcosmic vegetation canopies with different porosities and various velocities. Optical porosity and volumetric porosity were used to analyze the complicated morphology of vegetation canopy and influences of canopy structure parameters on flow resistance. The results show that, the relationship between drag coefficient and optical porosity can be approximately predicted as a quadratic polynomial, and a model expression is presented through simulation results. The model expression is validated using experimental results in a wind tunnel with tree branches. Based on the comparison of drag coefficient data obtained through the field experiments in real environment, this method provides a basis for predicting the real drag coefficient.

Key words:canopy; drag coefficient; porosity; microcosmic model