一種計(jì)算機(jī)模擬的針闊林木抗風(fēng)性能分析

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

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

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

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

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

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

1 試驗(yàn)方法

1.1 建立模型

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

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

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

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

1.2 湍流計(jì)算公式

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

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

(1)

(2)

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

(3)

(4)

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

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

對(duì)s運(yùn)用多元微分法[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)。不同樹(shù)種的葉面積密度不同,即使同為針葉樹(shù),柏樹(shù)和松樹(shù)的葉面積密度也有不同[16],分別為2.2和0.4 m2/m3。研究表明,同種闊葉樹(shù),楊樹(shù)在不同栽培策略下的葉面積密度會(huì)在6～10 m2/m3波動(dòng)[3]。綜合考慮,本研究使用了1,5和10 m2/m33個(gè)不同葉面積密度,對(duì)應(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

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

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

(9)

(10)

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

Sk=Pk+Lk

(11)

(12)

Lk=-2Cdfa|ui|k

(13)

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

(14)

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

(15)

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

(16)

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

Sε=Pε+Lε

(17)

(18)

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

(19)

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

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

(20)

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

(21)

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

Fd=Sk=Sε=Sω=0

(22)

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

1.3 網(wǎng)格化

網(wǎng)格化指將原來(lái)連續(xù)的求解區(qū)域劃分一系列不重復(fù)的有限單元[17],是對(duì)空間的離散化,便于下一步迭代求解。采用單元質(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 闊葉樹(shù)林與針葉樹(shù)林的網(wǎng)格圖Fig. 3 Meshing images of broadleaf forest versus coniferous forest

1.4 龍格-庫(kù)塔法

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

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

(23)

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

(24)

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

(25)

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

(26)

(27)

(28)

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

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

(29)

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

2 結(jié)果與分析

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

2.1 壓力分析

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

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

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

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

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

圖5 剖面位置和測(cè)量風(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分別對(duì)應(yīng)剖面1、剖面2、剖面3和剖面4。圖6 闊葉林與針葉林不同葉面積密度下的4個(gè)壓力剖面Fig. 6 Four pressure profiles at different leaf area densities in broadleaf and coniferous forests

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

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

圖7 闊葉樹(shù)林在不同湍流模型下剖面3的壓力云圖對(duì)比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更大,整個(gè)剖面的壓力上下限也更大,證實(shí)了k-ω湍流模型對(duì)邊界的細(xì)節(jié)處理更優(yōu)秀。

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

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

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

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

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

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

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

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

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

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

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

樹(shù)木模型見(jiàn)圖9a,在入口風(fēng)速為30 m/s,葉面積密度為1,5和10 m2/m3時(shí),兩種樹(shù)林模型在下風(fēng)側(cè)兩個(gè)距離樹(shù)林10 m(1倍樹(shù)高)和20 m(2倍樹(shù)高)的風(fēng)速數(shù)據(jù)見(jiàn)圖9b、c和d?？v坐標(biāo)表示測(cè)量高度,橫坐標(biāo)表示測(cè)量高度對(duì)應(yīng)的風(fēng)速。兩個(gè)測(cè)量風(fēng)速的具體位置參見(jiàn)圖5b中所對(duì)應(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個(gè)分圖表明葉面積密度越大,樹(shù)林擋風(fēng)的效果越好,二者是非線性關(guān)系:高葉面積密度數(shù)值再次提升時(shí),樹(shù)林的擋風(fēng)效果提升較小,趨于穩(wěn)定。從具體數(shù)值來(lái)看,葉面積密度由1 m2/m3增加到5 m2/m3時(shí),闊葉林針葉林下風(fēng)側(cè)10 m處最低風(fēng)速平均下降4.43 m/s;而葉面積密度由5 m2/m3增加到10 m2/m3時(shí),下風(fēng)側(cè)10 m處最低風(fēng)速平均下降3.73 m/s,變化較小。

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

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

3 結(jié)論與討論

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

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

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

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

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

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