袁星月,方文杰,張懷清,安鋒,云挺,4*

(1. 南京林業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院,南京 210037;2. 中國林業(yè)科學(xué)研究院,北京 100091; 3. 中國熱帶農(nóng)業(yè)科學(xué)院橡膠研究所,海口 571101;4. 南京林業(yè)大學(xué)林草學(xué)院,南京 210037)

風(fēng)是一種普遍存在的自然現(xiàn)象,樹木與風(fēng)之間常常發(fā)生復(fù)雜的相互作用。樹木能夠降低風(fēng)速,防止水土流失,而強風(fēng)能夠?qū)淠驹斐刹豢珊鲆暤膿p害[1]。研究樹木與風(fēng)之間的相互作用,不僅能夠評估臺風(fēng)對觀賞樹種與經(jīng)濟樹種造成的損害,還有助于深入了解樹林的防風(fēng)效應(yīng),以優(yōu)化防風(fēng)林的結(jié)構(gòu)。

樹林與風(fēng)的研究根據(jù)數(shù)據(jù)來源不同可以分為3類:衛(wèi)星數(shù)據(jù)研究、實地實驗與電腦模擬實驗。通過對比分析臺風(fēng)襲擊前后的森林衛(wèi)星數(shù)據(jù),能夠?qū)︼Z風(fēng)造成的森林損害進(jìn)行全面綜合的評估分析[2]。實地實驗要求研究人員前往風(fēng)害發(fā)生頻率較高的森林,等待自然疾風(fēng)吹過或人工創(chuàng)造強風(fēng)吹襲森林以采集數(shù)據(jù)[3]。風(fēng)洞實驗可以被視為一種模擬實地環(huán)境的實驗方法,建造風(fēng)洞放入移植的樹木或人工樹木模型,研究不同樹林擋風(fēng)效果背后的變化規(guī)律[4]。電腦模擬實驗主要利用計算流體動力學(xué)(computational fluid dynamics,CFD),通過計算機模擬獲得流體經(jīng)過樹木后的信息[5]。在電腦模擬實驗中,往往將樹木冠層看作均勻的多孔介質(zhì),通過調(diào)整模擬參數(shù)達(dá)到不同的模擬效果[6]。不研究樹冠每片樹葉狀態(tài)的原因較多,一方面,當(dāng)涉及整個樹冠中幾千片樹葉時,僅僅是從樹冠中分離出每片樹葉這一步就較為復(fù)雜[7],而且每片樹葉的狀態(tài)也難以預(yù)測。單片樹葉在風(fēng)中可能呈現(xiàn)彎曲、扭轉(zhuǎn)甚至卷曲等多種狀態(tài)[8],防風(fēng)效果又受樹葉生長年限、葉柄剛度和葉柄長度等因素影響[9]。另一方面,一些風(fēng)洞實驗與多孔介質(zhì)電腦模擬的組合實驗已證實,多孔介質(zhì)能夠替代眾多樹葉進(jìn)行符合現(xiàn)實情況的預(yù)測[10-11],故電腦模擬實驗中常采用多孔介質(zhì)模擬樹冠。

總的來說,目前樹木與風(fēng)相互作用的研究已經(jīng)能滿足大部分環(huán)境下的需要,但依然存在值得深入研究的方向。一方面,依托自然環(huán)境進(jìn)行的實驗不可控、耗時長且有一定危險性;另一方面,現(xiàn)實中樹木形態(tài)各異,單株樹木內(nèi)部的形態(tài)結(jié)構(gòu)等表型特征難以使用數(shù)學(xué)方式描述和計算。

本研究首先仔細(xì)觀察現(xiàn)實中的闊葉樹與針葉樹,總結(jié)它們各自的特點,進(jìn)行了精細(xì)化建模;然后利用CFD技術(shù),分別使用k-ε雙方程模型和k-ω雙方程模型進(jìn)行多次模擬實驗,模擬了闊葉樹與針葉樹附近的流場變化情況,記錄總壓力與風(fēng)速數(shù)據(jù);同時,多個截面也呈穿插狀設(shè)立在闊葉模型與針葉模型中,記錄局部的壓力與速度數(shù)據(jù),為后續(xù)的分析提供數(shù)據(jù)支持。該方法可控性強,能夠形象展現(xiàn)不同樹木附近的流場,分析不同樹林在強風(fēng)下內(nèi)部氣流狀況,了解它們受氣流影響與降低風(fēng)速的程度,為研究不同樹木與風(fēng)之間的相互作用提供一種直觀有效的方式。

1 試驗方法

1.1 建立模型

根據(jù)闊葉樹和針葉樹形態(tài)結(jié)構(gòu),建立了兩種樹木模型。闊葉樹模型參考的樹木品種為樟樹,屬于南方常見的行道樹,實際照片見圖1a;針葉樹模型參考的樹木品種為黑松,也是廣泛種植的優(yōu)質(zhì)觀賞類針葉樹,實際照片見圖1d。最終構(gòu)建的兩種樹木模型形態(tài)不同,但大小接近,高度均為10 m,且冠幅均接近6.3 m,能夠較好表現(xiàn)兩種樹木的區(qū)別,展現(xiàn)不同樹種的形態(tài)與其防風(fēng)效應(yīng)的關(guān)系。闊葉樹模型由主干、分支及主干和分支對應(yīng)的不規(guī)則葉團簇構(gòu)成,針葉樹由主干、分支及主干上附著的圓錐狀葉團簇構(gòu)成。將葉團簇視為整體,模擬成多孔介質(zhì)[12]。

b和e分別為實地照片a和d對應(yīng)的單株闊葉樹與 單株針葉樹模型;c和f分別為闊葉樹林與附近流體域、 針葉樹林與附近流體域的具體建模圖。圖1 闊葉樹林與針葉樹林的實地照片與具體建模圖Fig. 1 Field photographs and concrete modeling images of broadleaf and coniferous forests

將9棵同種樹木模型以6.5 m為間距排列為3×3的方陣,以分析不同樹木組成的樹林特點。闊葉樹針葉樹的實際照片和對應(yīng)建模如圖1所示。

模擬時,風(fēng)從ZOX平面出發(fā),向Y軸正方向流動。記9棵樹中離ZOX平面最近的3棵樹木為第1排樹木,離ZOX平面較遠(yuǎn)和最遠(yuǎn)的兩排樹木分別為第2排和第3排樹木。因為直面強風(fēng),第1排樹木的狀態(tài)可以作為現(xiàn)實森林中迎風(fēng)面樹木的參考。類似的,中間的1棵樹(被其余8棵樹包圍)可以看成強風(fēng)作用下森林中心樹木的狀況,第3排樹木可以模擬森林背風(fēng)面樹木的狀況。

1.2 湍流計算公式

湍流是流體的一種流動狀態(tài),流體的流動可以根據(jù)雷諾數(shù)的大小分為層流和湍流兩種:層流指穩(wěn)定且有序的流動,而湍流則是無規(guī)則、復(fù)雜的流動。雷諾數(shù)描述了流體在慣性力和黏性力之間的相對性,當(dāng)流體速度較小、黏性較高時,雷諾數(shù)較小,流動趨向于層流,反之則接近湍流。

本試驗主要研究強風(fēng)狀態(tài)下風(fēng)樹的相互作用,設(shè)定的風(fēng)速高達(dá)30 m/s。高速流場中以湍流為主,幾乎不出現(xiàn)層流,所以只考慮湍流,選擇雷諾時均方程(Reynolds averaged navier stokes equations,RANS)模擬湍流相關(guān)參數(shù),即重點關(guān)注一段時間內(nèi)流體的平均流動狀態(tài)。雷諾時均方程是對納維-斯托克斯方程(navier stokes equations,NS)的改進(jìn),廣義上納維-斯托克斯方程由三大守恒控制方程組成,即質(zhì)量守恒方程、動量守恒方程和能量守恒方程。常溫空氣作為一種不可壓縮流體,流動時熱量變化很小,通常不考慮能量方程。雷諾時均方程中的質(zhì)量守恒方程和動量守恒方程由公式(1)和公式(2)描述[13]:

(1)

(2)

式中:ui(或uj)為xi(或xj)下的速度;i和j分別指空間直角坐標(biāo)系的任意兩個坐標(biāo)軸的方向(i,j=1,2,3);t為時間;P為壓力;ρ為空氣密度。公式(2)中用到了一種簡化數(shù)學(xué)表達(dá)的約定——愛因斯坦求和約定(Einstein summation convention),即使用一個項中的重復(fù)指標(biāo)代表對指標(biāo)的所有可能值求和[14],使得公式(2)實際包含了3個對應(yīng)不同方向的公式。公式(2)中的有效運動渦流黏度νeff和關(guān)于阻力系數(shù)的附加項Fd由公式(3)和公式(4)計算。

(3)

(4)

式中:ν為動力黏度;νt為湍流運動黏度;k為湍流動能;ε為耗散速率;Cμ=0.09為經(jīng)驗常數(shù);Cdf為阻力系數(shù);a為葉面積密度(leaf area density,LAD),單位體積的總植物葉面積,m2/m3。阻力系數(shù)Cdf與a有關(guān),使用收集到的N個樣本數(shù)據(jù)[2]進(jìn)行擬合可得到擬合函數(shù)。

前人針對阻力系數(shù)的研究表明,阻力系數(shù)隨葉面積密度上升呈現(xiàn)先上升后平緩的趨勢,單純使用多項式擬合結(jié)果較差,所以選取以自然常數(shù)為底數(shù)的指數(shù)函數(shù)為擬合函數(shù),該函數(shù)具有漸近線,符合阻力系數(shù)與葉面積密度的基本趨勢。設(shè)f(a)為擬合函數(shù),則有Cdf-f(a)為誤差函數(shù),使誤差平方和s最小化。

對s運用多元微分法[15]使其最小化,即要求s各參數(shù)的偏導(dǎo)等于0,得到方程組:

(5)

(6)

(7)

求解該方程組,得到該方程的最優(yōu)m0、m1、m2、m3、m4和m5分別為0.597 1,-1.235 0,1.712 0,3.943 0,-0.512 5,1.893 0。阻力系數(shù)Cdf與葉面積密度的最終擬合函數(shù)為:

(8)

阻力系數(shù)Cdf隨葉面積密度a的增大而增大,當(dāng)達(dá)到0.6后阻力系數(shù)基本保持不變(圖2)。不同樹種的葉面積密度不同,即使同為針葉樹,柏樹和松樹的葉面積密度也有不同[16],分別為2.2和0.4 m2/m3。研究表明,同種闊葉樹,楊樹在不同栽培策略下的葉面積密度會在6～10 m2/m3波動[3]。綜合考慮,本研究使用了1,5和10 m2/m33個不同葉面積密度,對應(yīng)的阻力系數(shù)Cdf分別為0.53,0.60 和0.60。

圖2 阻力系數(shù)與葉面積密度的擬合函數(shù)曲線Fig. 2 Fitted function plot of drag coefficient and leaf area density

雷諾時均方程中的湍流運動黏度,即式(3)中的νt無法直接計算,需要引入湍流模型進(jìn)行求解。湍流運動黏度可由k和ε共同推導(dǎo),共使用到k-ε雙方程模型和k-ω雙方程模型兩種湍流模型,都是針對樹木冠層改進(jìn)過的[5]。k-ε雙方程模型是目前應(yīng)用最廣泛的湍流模型,能較好地模擬計算湍流的發(fā)展?fàn)顩r;k-ω雙方程模型中,ω是比湍流耗散率,該湍流方程更加注重邊界層的細(xì)節(jié),同時需要更多計算資源。本研究主要使用k-ε雙方程模型,在使用到k-ω雙方程模型時會額外指出。

本研究使用的k方程和ε方程增加了與樹冠相關(guān)的附加項,具體的k方程和ε方程如公式(9)和公式(10)所示[13]:

(9)

(10)

式中:經(jīng)驗常數(shù)Cε1=1.42,Cε2=1.92;常數(shù)σk=1.0,σε=1.3;Sk和Sε分別是k方程和ε方程中與樹冠相關(guān)的附加項。流體穿過類似樹冠的多孔介質(zhì)時,發(fā)生了尾流產(chǎn)生和能量減少兩種現(xiàn)象,故k方程附加項Sk等于尾流產(chǎn)生項Pk和能量損失項Lk之和。

Sk=Pk+Lk

(11)

(12)

Lk=-2Cdfa|ui|k

(13)

耗散速率ε可以表示成湍流動能k的函數(shù)關(guān)系式:

(14)

式中:Cε=0.16,為經(jīng)驗常數(shù);l為混合長度。若l取為定值,則公式(14)兩邊對k求導(dǎo)變形可得到k和ε的具體數(shù)值關(guān)系,如公式(15)所示:

(15)

k方程附加項Sk,即公式(11～13),可與公式(15)的最終形式聯(lián)系起來,得到公式(10)中所示ε方程附加項Sε:

(16)

但l的大小會因為方程計算對象所處位置的變化而改變,公式(16)具有一定的局限性,故研究人員[13]通過研究使用經(jīng)驗系數(shù)修正了ε方程附加項Sε:

Sε=Pε+Lε

(17)

(18)

Lε=-4Cpε2Cdfa|ui|ε

(19)

式中,Cpε1=1.8,Cpε2=0.6,均為經(jīng)驗常數(shù)。

關(guān)于k-ω雙方程模型,k方程與k-ε雙方程模型相同,ω方程見公式(20):

(20)

常數(shù)σω=0.5,ω方程中對應(yīng)的與樹冠相關(guān)的附加項Sω見公式(21):

(21)

樹冠相關(guān)的附加項Fd只在冠層內(nèi)部存在,在冠層之外的流場域中本節(jié)的各附加項都記為0:

Fd=Sk=Sε=Sω=0

(22)

不同研究往往在附加項上進(jìn)行研究與變動以完成針對樹冠的改進(jìn)與優(yōu)化,本研究提及的Fd、Sk、Sε和Sω4個附加項與樹冠參數(shù)密切相關(guān),屬于對原始雷諾時均方程的重要改進(jìn),在樹冠模擬實驗中有著重要作用。

1.3 網(wǎng)格化

網(wǎng)格化指將原來連續(xù)的求解區(qū)域劃分一系列不重復(fù)的有限單元[17],是對空間的離散化,便于下一步迭代求解。采用單元質(zhì)量用于表現(xiàn)網(wǎng)格劃分質(zhì)量,單元質(zhì)量越大越好,一般不能低于0.7。闊葉林單元總數(shù)為4 176 642,平均單元質(zhì)量為0.834 54;針葉林單元總數(shù)為14 071 152,平均單元質(zhì)量為0.838 66。具體網(wǎng)格如圖3所示。

圖3 闊葉樹林與針葉樹林的網(wǎng)格圖Fig. 3 Meshing images of broadleaf forest versus coniferous forest

1.4 龍格-庫塔法

龍格-庫塔法是求解微分方程的一種常見數(shù)值解法[18],在工程上應(yīng)用廣泛[19]。與解析解相對,數(shù)值解指通過使用數(shù)值計算方法對問題進(jìn)行近似求解得到結(jié)果[20]。解析解的運算耗時極長,還存在無法通過數(shù)學(xué)分析獲得解析解的情況,所以數(shù)值解使用廣泛,對湍流模擬有重要意義。

四階龍格庫塔法作為高階龍格-庫塔法,具有高精度、低誤差的優(yōu)點。以k-ε湍流模型為例,式(1)、(2)、(9)和(10)共同構(gòu)成了關(guān)于湍流求解的微分方程組,運用四階龍格-庫塔法,對微分方程組進(jìn)行求解。需要求解的微分方程和已知的初值簡化如公式(23)所示:

(23)

式中:A表示空間x、時間t中任意一個變量;B表示壓力P、速度u、動能k和耗散率ε中任意一個變量;c1和c2為對應(yīng)的初始值,整個微分方程中初始值是完全已知的。f為解關(guān)于湍流的微分方程組得到的關(guān)于B的一階導(dǎo)數(shù):

(24)

龍格-庫塔法屬于迭代法,以一次迭代為例,第n+1次迭代中,需要將第n次迭代中的已知數(shù)據(jù)Bn代入公式(25):

(25)

則可得到下一個迭代的Bn+1值,該值是由Bn加上步長和估算斜率的乘積得到的。K1、K2、K3和K4分別是用于計算估算斜率的開始斜率、兩個中間斜率和終點斜率。當(dāng)對時間偏導(dǎo)時,λ作為時間間隔,是人為選取的;對空間偏導(dǎo)時,λ由網(wǎng)格給出。公式(24)是由4個方程構(gòu)成的方程組,其中前3個方程分別是速度、動能和耗散率關(guān)于時間的偏導(dǎo)。它們代入公式(25)后,具體方程見公式(26)～(28)。

(26)

(27)

(28)

已知所有網(wǎng)格的壓力初始值,在求得某一個網(wǎng)格下一個迭代的新壓力數(shù)值后,可以通過推導(dǎo)計算出附近網(wǎng)格的壓力值,進(jìn)而得知所有網(wǎng)格的新壓力值,以此不斷迭代得到收斂解[21]。

對單個網(wǎng)格而言,求得速度后,求解壓力需要對公式(24)方程組最后一個公式求偏導(dǎo),可得到壓力泊松(Poisson)方程[22]見公式(29),推出下一迭代的壓力值:

(29)

每次迭代中各變量會進(jìn)行計算得到更新。重復(fù)迭代,當(dāng)某一循環(huán)下同一網(wǎng)格的壓力變化與上一循環(huán)相比小于0.1%或迭代循環(huán)次數(shù)達(dá)到300時,即退出循環(huán)。可得到任意時空下的壓力P、速度u、動能k和耗散率ε數(shù)值。

2 結(jié)果與分析

根據(jù)GB/T 28591—2012《風(fēng)力等級》,12級風(fēng)被定義為臺風(fēng),又稱一級颶風(fēng),風(fēng)速區(qū)間為32.7～37.0 m/s。本研究將入口風(fēng)速設(shè)置為30 m/s,兩種樹木模型與3種葉面積密度組合進(jìn)行6次實驗。

2.1 壓力分析

闊葉林與針葉林不同葉面積密度下的總壓力剖面云圖見圖4。圖4中每個分圖都包含10個壓力剖面云圖,這些半透明剖面以3.25 m為間距排列,貫穿整個樹林。

a～c)分別是葉面積密度為1,5和10 m2/m3的闊葉林結(jié)果; d～f)分別是葉面積密度為1,5 和10 m2/m3的針葉林結(jié)果。圖4 闊葉林與針葉林不同葉面積密度下的總壓力剖面云圖Fig. 4 Total pressure profile cloud at different leaf area densities of broadleaf forest and coniferous forest

總體來看,不同葉面積密度下,闊葉林與針葉林附近的氣流壓力變化趨勢類似。壓力極大值區(qū)域的出現(xiàn)位置都為迎風(fēng)面第一排樹的前方,大范圍低壓區(qū)都出現(xiàn)在了樹林后方。遇到物體阻礙后,氣流在物體前方產(chǎn)生高壓區(qū),在物體后方產(chǎn)生低壓區(qū),整體來看對物體施加與流動方向同向的力。

隨著葉面積密度的上升,最前排樹木承受的壓力不斷增大,樹林內(nèi)部的氣壓差也不斷上升。面臨同等強風(fēng),枝葉稀疏的樹木比枝葉茂密的樹木承受的壓力更小,不容易受損。強風(fēng)吹落大量樹葉后,樹木整體狀態(tài)會趨于穩(wěn)定。

為了更好地展示不同位置的壓力情況,從圖4每個分圖中單獨取出了4個剖面,分別是迎風(fēng)面第1排樹木前的剖面1(距第1排樹主干3.25 m),第1排樹木和第2排樹木之間的剖面2(距第1排和第2排樹主干均為3.25 m),第2排樹木和第3排樹木之間的剖面3以及第3排樹木后面的剖面4(距第3排樹主干3.25 m)。以葉密度最高的闊葉樹林為例,4個剖面的位置示意圖見圖5a,不同樹木模型中四剖面的位置相同。從下風(fēng)側(cè)觀測圖5a所示四剖面的壓力云圖,不同實驗的壓力剖面圖如圖6所示。

圖5 剖面位置和測量風(fēng)速位置的示意圖Fig. 5 Schematic diagram of the sections and the measured wind speed position

a～c)葉面積密度為1,5和10 m2/m3的闊葉林結(jié)果;d～f)葉面積密度為1,5和10 m2/m3的針葉林結(jié)果。 1、2、3和4分別對應(yīng)剖面1、剖面2、剖面3和剖面4。圖6 闊葉林與針葉林不同葉面積密度下的4個壓力剖面Fig. 6 Four pressure profiles at different leaf area densities in broadleaf and coniferous forests

從左到右觀察圖6,壓力剖面1中間壓力大周圍壓力小,后3個壓力剖面中間壓力小周圍壓力大。剖面4展現(xiàn)的高壓力區(qū)域邊界與樹木樹冠本身的形態(tài)邊界類似,闊葉樹為不規(guī)則圓狀而針葉樹是上小下大的箭頭狀,說明壓力的產(chǎn)生區(qū)域與承受流體的多孔介質(zhì)形態(tài)有一定的關(guān)聯(lián)。在形狀相似的單樹構(gòu)成的樹林中,樹林邊緣的樹木承受的壓力低于樹林內(nèi)部的樹木,更易受損。這一特點在圖6e3中表現(xiàn)最為明顯,樹林邊緣的樹木附近壓力僅為1 210 Pa,樹林內(nèi)部的樹木附近壓力達(dá)到1 234 Pa。

為了進(jìn)一步分析,選擇最高葉面積密度的闊葉樹林中的壓力剖面3,即圖6c3對應(yīng)的情況,進(jìn)行k-ω湍流模型模擬,闊葉樹林兩種湍流模型的剖面3壓力云圖對比結(jié)果見圖7。

圖7 闊葉樹林在不同湍流模型下剖面3的壓力云圖對比Fig. 7 Comparison of pressure cloud image in Section 3 of broad-leaved forest under different turbulence models

圖7中兩張分圖的壓力上下限及壓力分布類似,圖7b捕捉到了更多細(xì)節(jié)。圖7b的低壓區(qū)面積比圖7a更大,整個剖面的壓力上下限也更大,證實了k-ω湍流模型對邊界的細(xì)節(jié)處理更優(yōu)秀。

綜合統(tǒng)計模型總體數(shù)據(jù)與4個剖面數(shù)據(jù),可以得到表1所示的各區(qū)域壓力數(shù)據(jù)。

不同葉面積指數(shù)的針葉林平均總模型壓力極大值(2 112.70 Pa)大于闊葉林(1 898.73 Pa),在針葉模型樹冠邊角處易產(chǎn)生局部高壓區(qū)。闊葉林總模型整體壓力差大于針葉林,闊葉林內(nèi)部氣流比針葉林更混亂,抗風(fēng)能力較弱,這可能跟闊葉樹樹冠體積較大有關(guān);分析表1數(shù)據(jù)可以發(fā)現(xiàn),各剖面的壓力極值沒有負(fù)壓,而總模型存在負(fù)壓。壓力為負(fù)值表明壓力低于參考壓力,該區(qū)域產(chǎn)生了吸力,附近的枝干與葉團簇更容易受損?？偰Ｐ拓?fù)壓區(qū)集中在樹干后方極小區(qū)域,該區(qū)域壓力急劇下降,和現(xiàn)實中樹木在強力風(fēng)荷載下枝條易彎折的情況是吻合的。最后觀察各剖面之間的聯(lián)系:相比其他剖面,剖面1有最高的平均壓力極大值(1 690.91 Pa)、最高的平均壓力極小值(1 283.96 Pa)和最大的平均壓力差(406.95 Pa)。后幾個剖面的壓力極值和壓力差與剖面1相比均有下降且后幾個剖面的數(shù)值較為接近,說明當(dāng)狂風(fēng)來襲時樹林迎風(fēng)面第一排樹木承受了最大的壓力。

表1 CFD模擬中的各區(qū)域壓力數(shù)據(jù)Table 1 Regional pressure data in CFD simulation

2.2 樹林內(nèi)部風(fēng)速分析

闊葉林與針葉林風(fēng)速矢量圖見圖8。圖8最后一排樹木的枝干處和樹林后方均出現(xiàn)細(xì)小的藍(lán)色箭頭,表現(xiàn)了風(fēng)速的降低。黃色箭頭聚集在整體樹林的四周,即樹林底部、側(cè)面與樹林樹冠上方的風(fēng)速都有提升,這是一種基礎(chǔ)的流動情況[23]。當(dāng)流體繞過物體流動時,因為存在物體的阻礙,流體的流通面積變小了,所以在物體側(cè)面,流體需要以更快的速度通過,才能保證和上游流通能力相同。

圖8 闊葉林與針葉林速度矢量圖Fig. 8 Broadleaf and coniferous forest velocity vector plots

矢量圖也展現(xiàn)了闊葉樹和針葉樹的區(qū)別。由于樹冠結(jié)構(gòu)不同,氣流經(jīng)過針葉樹樹冠頂部尖端受到的影響更為突出。遠(yuǎn)離樹冠的地方風(fēng)速高,越靠近樹冠風(fēng)速越低。闊葉樹樹冠由多個葉團簇組合而成,且不同葉團簇之間有空隙,所以在葉團簇之間,氣流的影響非常明顯,迎風(fēng)面第1排闊葉樹的葉團簇之間有大量黃綠色較高速的矢量箭頭,表明闊葉林內(nèi)部風(fēng)速分布相對不均勻,波動大。

在圖5a所示四剖面的位置再次設(shè)置4個風(fēng)速剖面用以采集風(fēng)速數(shù)據(jù),不同葉面積密度下不同模型不同位置的對應(yīng)風(fēng)速數(shù)據(jù)見表2。

表2 CFD模擬中的各區(qū)域最高風(fēng)速數(shù)據(jù)Table 2 Maximum speed data for each region in CFD simulation

總模型中最大風(fēng)速出現(xiàn)在第1排樹的樹干附近。第1排樹的樹干直接暴露在強風(fēng)中,因此在它們附近出現(xiàn)了高風(fēng)速區(qū)。而由于樹木的阻礙作用,各剖面的最大風(fēng)速都出現(xiàn)在樹林側(cè)面和樹冠上方。

2.3 下風(fēng)側(cè)風(fēng)速垂直分布圖分析

樹木模型見圖9a,在入口風(fēng)速為30 m/s,葉面積密度為1,5和10 m2/m3時,兩種樹林模型在下風(fēng)側(cè)兩個距離樹林10 m(1倍樹高)和20 m(2倍樹高)的風(fēng)速數(shù)據(jù)見圖9b、c和d?？v坐標(biāo)表示測量高度,橫坐標(biāo)表示測量高度對應(yīng)的風(fēng)速。兩個測量風(fēng)速的具體位置參見圖5b中所對應(yīng)的綠色和藍(lán)色線位置。

圖9 下風(fēng)側(cè)風(fēng)速垂直分布圖Fig. 9 Vertical wind speed distribution on the leeward side

整體而言,圖9中b、c和d 3個分圖表明葉面積密度越大,樹林擋風(fēng)的效果越好,二者是非線性關(guān)系:高葉面積密度數(shù)值再次提升時,樹林的擋風(fēng)效果提升較小,趨于穩(wěn)定。從具體數(shù)值來看,葉面積密度由1 m2/m3增加到5 m2/m3時,闊葉林針葉林下風(fēng)側(cè)10 m處最低風(fēng)速平均下降4.43 m/s;而葉面積密度由5 m2/m3增加到10 m2/m3時,下風(fēng)側(cè)10 m處最低風(fēng)速平均下降3.73 m/s,變化較小。

所有數(shù)據(jù)趨勢相似:風(fēng)在經(jīng)過樹木模型樹冠層(高度2～10 m)附近時,速度明顯降低。樹冠層上方,風(fēng)速逐漸上升又回到入口風(fēng)速的30 m/s左右,13 m高度以上幾乎看不到樹木對風(fēng)的影響。2～10 m是樹木冠層的主要分布區(qū),風(fēng)速顯著降低且離樹林越近降低越多。12 m高度處風(fēng)速略高于入口面的30 m/s風(fēng)速,這是阻攔引起的加速現(xiàn)象。

風(fēng)速在20 m測量點處的變化幅度小于10 m測量點處,并且樹木模型的種類也影響風(fēng)速的垂直分布。闊葉樹整體冠層趨近于球狀,而針葉樹的樹冠形態(tài)近似圓錐,冠層重心偏下。所以葉面積密度較低時(圖9b),通過兩種樹木的風(fēng)分別在約6和4 m高度處達(dá)到最低速度。但隨著葉面積密度的上升,風(fēng)速的極小值出現(xiàn)的高度都向5 m高度處移動,樹木冠層的形態(tài)對樹后氣流的影響逐漸減弱。

3 結(jié)論與討論

本研究建立了闊葉樹和針葉樹的樹木模型,并保留了兩種樹木的形態(tài)特征,改變?nèi)~面積密度進(jìn)行多次實驗,定量分析比較不同條件下壓力速度的分布情況。

本研究結(jié)論在兩方面與現(xiàn)有實驗方法的結(jié)果相似。首先在損傷性評價方面,模擬結(jié)果顯示葉面積密度數(shù)值高的樹木承受的壓力更大,在風(fēng)中的受損概率也比葉面積密度低的樹木大,更易發(fā)生倒伏。說明闊葉林內(nèi)部氣流比針葉林內(nèi)部氣流混亂,闊葉林內(nèi)部的局部最大風(fēng)速大于針葉樹內(nèi)部,更易受損?？紤]到同等樹高和冠幅下,闊葉樹的樹冠體積大于針葉樹,也就是說樹木的形態(tài)特征可能是判斷風(fēng)損害樹木程度的重要因素,樹冠體積小的樹木抗風(fēng)性能優(yōu)于樹冠形態(tài)不規(guī)則且體積大的樹木,更易發(fā)生枝干彎折與葉片脫落,這與華南地區(qū)及安徽農(nóng)業(yè)大學(xué)風(fēng)災(zāi)后實地統(tǒng)計數(shù)據(jù)相符[24-25]。

其次在擋風(fēng)效果方面,下風(fēng)側(cè)風(fēng)速垂直分布數(shù)據(jù)給出了定量的評價,不同葉面積密度下,各樹林下風(fēng)側(cè)1倍樹高的地方,最低風(fēng)速與入口風(fēng)速相比下降了60.13%～98.51%。Poh等[6]使用的實地測量數(shù)據(jù)顯示,針葉樹防風(fēng)林下風(fēng)側(cè)1倍樹高的地方最低風(fēng)速下降了約60%,該實測數(shù)據(jù)在本研究模擬數(shù)據(jù)范圍內(nèi)。類似的風(fēng)洞實驗中,鄭波等[26]使用定制的闊葉林模型與通用樹林模型多次排列成不同樹林測得各樹林下風(fēng)側(cè)一倍樹高的地方,最低風(fēng)速與入口風(fēng)速相比下降了68.5%～98.9%,與本研究模擬的上下限誤差分別為8.37%和0.39%,這可能是不同實驗環(huán)境中葉面積密度與樹林行數(shù)不同導(dǎo)致的正常波動。在改變?nèi)~面積密度的風(fēng)洞實驗中,Cao等[4]對同一樹木逐步摘除葉片后進(jìn)行多組風(fēng)洞實驗,得到的數(shù)據(jù)顯示:葉面積密度越大,樹木擋風(fēng)效果越好,但樹木擋風(fēng)效果增加到一定程度后,再提升葉面積密度,擋風(fēng)效果改善變緩,與本研究實驗結(jié)論呈現(xiàn)一致的趨勢。

本研究方法具有如下優(yōu)點:首先,電腦模擬實驗繞過了測速器材的限制,能夠得到樹林中每個角落的詳細(xì)數(shù)據(jù);其次,入口風(fēng)速靈活可控,能夠為實驗創(chuàng)造出適合的環(huán)境,降低時間和資金投入。最后,樹木模型體現(xiàn)了闊葉樹的分支結(jié)構(gòu)和針葉樹的樹冠形態(tài),能夠較好地展現(xiàn)樹林內(nèi)部的氣流狀態(tài)。本次實驗驗證了k-ε湍流模型和k-ω湍流模型的差異性,產(chǎn)生了不同葉面積密度下的多組壓力速度數(shù)據(jù),為樹木在風(fēng)力作用下的行為提供了數(shù)據(jù)支持。這些結(jié)果能夠幫助進(jìn)一步深入研究樹木風(fēng)力響應(yīng)的機制,推測樹木的結(jié)構(gòu)變化和振動特性,對于優(yōu)化數(shù)值模擬方法以及應(yīng)用于實際工程和環(huán)境中也具有重要意義。

本方法也存在一些可進(jìn)一步優(yōu)化的方向:樹木模型還有細(xì)化的空間,可以嘗試構(gòu)建得更精細(xì);研究對象都是由同種樹木構(gòu)成的純林,可以擴展到針葉混交林或更為復(fù)雜的森林模型,將灌木等低矮樹種納入模型;通過引用更貼合實例的葉面積密度數(shù)據(jù)開展更有針對性的研究。

研究樹木在強風(fēng)作用下的各類參數(shù)在防護林、綠化帶等設(shè)施的構(gòu)建中有指導(dǎo)意義。本研究基于雷諾時均方程、改進(jìn)的k-ε湍流模型和k-ω湍流模型,采用CFD模擬,研究了不同葉面積密度下不同樹木模型的抗風(fēng)表現(xiàn)。結(jié)果表明,葉面積密度對樹林阻風(fēng)效果的影響大于樹木模型對樹林阻風(fēng)效果的影響。隨著葉面積密度的增大,樹林阻風(fēng)效果不斷增強。葉面積密度由1 m2/m3增加到5 m2/m3時,下風(fēng)側(cè)10 m處闊葉樹風(fēng)速的最大下降百分比由60.13% 提高到77.70%,針葉樹由79.50%提高到91.36%。針葉樹的擋風(fēng)能力強于闊葉樹,同在下風(fēng)側(cè)10 m處且葉面積密度同為1 m2/m3時,針葉林的最低風(fēng)速為闊葉林最低風(fēng)速的51.43%。本研究使用的方法可以通過調(diào)整參數(shù)與建模,擴展細(xì)分到不同種類的闊葉林與針葉林;結(jié)果可用于分析確定森林在大風(fēng)條件下的各類參數(shù),預(yù)測強風(fēng)襲擊后樹木的損毀情況,對不同森林種植策略下的風(fēng)損害進(jìn)行量化分析。同時可以幫助建設(shè)更合理的防風(fēng)林,有利于更好地節(jié)省資金,在強風(fēng)下保護設(shè)施。