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

王冰清，付海明

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

?

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

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

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

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

關鍵詞：樹冠； 阻力系數； 孔隙率； 微觀模型

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

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

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

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

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倍樹冠高度，故可忽略流域高度對流場求解結果的影響.為保證入口邊界和出口邊界處流動的充分發展，設定樹冠到流體入口邊界距離為0.625 m，即2.5倍樹冠高，樹冠到流體出口邊界距離為0.625 m，亦為2.5倍樹冠高.流域被劃分為8 385 633個非均勻四面體網格，并對樹冠內部作局部網格加密處理.邊界條件：流域的底部設置為粗糙壁面，粗糙度取0.005 m,并假設壁面無明顯的障礙物和植被.除了入口和出口以外，計算流域的其他壁面包括樹枝和樹葉的表面均設為無滑移壁面邊界條件.入口設置為速度入口邊界條件，出口設置為相對壓力為0的壓力出口.事實上，樹枝的表面是粗糙的.但是之前的研究中沒有關于樹枝表面粗糙度的數據，文獻[21]研究得到的樹枝表面粗糙度為0.002 m，將其作為參考用于代表全尺寸的樹木表面粗糙度.流域中的三維微觀樹冠模型是由VBA生成的，其樹冠中的樹枝葉子是由計算機編程產生的，其外形及特征與真實樹冠十分相似，每一個虛擬樹冠都具有不同的隨機分布結構，可近似認為本文的虛擬樹冠是“真實”的，如圖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風洞試驗

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

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

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

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

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

(1)

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

3結果與討論

3.1樹冠流場特性

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

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

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

3.2孔隙率

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

表1　孔隙率測試結果

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

(2)

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

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

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

3.3阻力研究

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

(a) u=3 m/s

(b) u=10 m/s

different wind speed

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

圖7　光學孔隙率與阻力系數關系圖Fig.7　The relationship between Po and Cd

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

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

(3)

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

Cd=-4.46Po2+3.72Po+0.843

(4)

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

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

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

Su=-Cd×ρAfu2

(5)

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

樹冠阻力Δp的表達式為

Δp=Cd×ρu2/2

(6)

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

3.72Po+0.843)×ρu2/2

(7)

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

3.72Po+0.843)×ρu2/2

(8)

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

圖9　試驗數據與經過修正的擬合結果對比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為試驗修正系數，校正試驗中所產生的誤差，其數值為0.62～1.25.光學孔隙率小的修正系數選取較小數值，光學孔隙率大的修正系數選取較大的數值.

4結語

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

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

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

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

參考文獻

[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—)，女，遼寧撫順人，碩士研究生，研究方向為空氣過濾凈化與空氣品質. E-mail：wbq1827@163.com 付海明(聯系人), 男，高級工程師，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