植物葉片表面捕集氣溶膠粒子特性數(shù)值分析

賈孟霞，付海明

(東華大學(xué) 環(huán)境科學(xué)與工程學(xué)院，上海 201620)

隨著工業(yè)化、城市化進(jìn)程加快，工業(yè)廢氣、汽車尾氣、建筑粉塵的排放量顯著增加，大氣氣溶膠粒子污染日益加劇，給人類的生產(chǎn)、生活帶來巨大影響。在城市中合理種植綠色植物是降低粒子濃度較為經(jīng)濟(jì)有效的方式。宏觀上，植物表現(xiàn)出一定的“滯塵”效應(yīng)，即葉片可滯留并轉(zhuǎn)移粒子；微觀上，粒子通過布朗擴(kuò)散、攔截、慣性碰撞及重力沉降等機(jī)理沉積在葉片表面，進(jìn)而脫離大氣環(huán)境[1-2]。因此，合理地預(yù)測(cè)植物葉片的滯塵效果及其影響因素對(duì)城市規(guī)劃設(shè)計(jì)具有重要理論意義及實(shí)用價(jià)值。

針對(duì)植物“滯塵”現(xiàn)象的研究已有較多報(bào)道。一些學(xué)者主要通過現(xiàn)場(chǎng)測(cè)量的方式定量分析葉片捕集粒子的效果，例如：俞莉莉等[3]對(duì)揚(yáng)州市綠色植物滯塵效益進(jìn)行研究后發(fā)現(xiàn)，懸鈴木等葉面越粗糙的植物捕集粒子的能力越強(qiáng)；趙松婷等[4]對(duì)北京市29種園林植物進(jìn)行葉片采樣，通過電鏡分析評(píng)估粒子的滯留能力。另外一些學(xué)者重點(diǎn)關(guān)注了風(fēng)洞試驗(yàn)及植物捕集粒子效率的參數(shù)化表達(dá)，例如：Chamberlain[5]對(duì)人工黏草沉積粒子進(jìn)行風(fēng)洞試驗(yàn)，提出將捕集效率作為衡量植物滯塵能力強(qiáng)弱的依據(jù)；Slinn[6]基于文獻(xiàn)[5]的風(fēng)洞試驗(yàn)數(shù)據(jù)對(duì)植物冠層的捕集效率進(jìn)行擬合，得到了布朗擴(kuò)散、攔截、慣性碰撞作用下粒子捕集效率的參數(shù)化表達(dá)式；Zhang等[7]指出空氣濕度會(huì)影響粒子的沉積，并對(duì)不同品種植物的冠層捕集效率參數(shù)化表達(dá)式進(jìn)行優(yōu)化；Petroff等[8]對(duì)比了從森林、草地上測(cè)得的粒子沉積數(shù)據(jù)與現(xiàn)有的參數(shù)化模型模擬數(shù)據(jù)，結(jié)果表明參數(shù)化模型具有一定優(yōu)越性；Zhang等[9]基于顆粒動(dòng)力學(xué)分析方法進(jìn)行風(fēng)洞試驗(yàn)，測(cè)量了植物表面粒子的沉積速度，結(jié)果表明沉積效率受粒子粒徑和風(fēng)速的影響。

綜上所述，現(xiàn)有研究主要集中在宏觀層面下植物捕集粒子的效率預(yù)測(cè)及影響因素探究方面，而鮮有關(guān)于微觀尺度下單個(gè)植物葉片對(duì)粒子捕集行為影響的研究報(bào)道。本文在相似理論基礎(chǔ)上采用拉格朗日法追蹤粒子的運(yùn)動(dòng)軌跡，并在微觀尺度下研究植物單葉片表面對(duì)微細(xì)粒子的捕集效率，探究粒子粒徑、風(fēng)速和迎風(fēng)角對(duì)葉片捕集效率的影響，并嘗試給出葉片捕集粒子的效率表達(dá)式，以期為植物改善城市空氣質(zhì)量環(huán)境提供理論參考。

1 物理問題及研究方法

1.1 流場(chǎng)模型

植物葉片的形狀多種多樣，本文選擇較為常見的扁平形葉片，將其置于長方體流場(chǎng)中，假設(shè)風(fēng)向垂直于葉片。鑒于三維流場(chǎng)的運(yùn)算較為復(fù)雜，截取植物葉片所在的xy平面將流場(chǎng)簡化為二維。植物葉片捕集粒子的流場(chǎng)及計(jì)算域如圖1所示，其中，葉片厚度為0.5 cm，長度l為20 cm。基于l對(duì)葉片中心與計(jì)算域邊界間的距離進(jìn)行無量綱處理，得到葉片中心與速度入口的距離為4l，與壓力出口的距離為8.5l，與上、下邊界的距離均為3.75l。迎風(fēng)角θ為葉片中線與水平軸x軸的夾角。

圖1 植物葉片捕集粒子的流場(chǎng)及計(jì)算域Fig.1 Flow field and computational field of plant leaf capture particles

根據(jù)植物葉片的長度l、風(fēng)速v以及運(yùn)動(dòng)黏度μ計(jì)算雷諾數(shù)Re可知，流場(chǎng)中空氣的運(yùn)動(dòng)處于紊流狀態(tài)。為簡化計(jì)算，假設(shè)其為二維定常不可壓縮流動(dòng)，滿足的連續(xù)性方程和Navier-Stokes方程如式(1)～(3)[10]所示。

(1)

(2)

(3)

式中：vx、vy分別為x、y軸方向的流體速度分量，m/s；μ為空氣動(dòng)力學(xué)黏度，Pa·s；ρ為空氣密度，kg/m3；p為空氣壓強(qiáng)，Pa。

由于葉片、流場(chǎng)的尺寸與粒子粒徑的數(shù)量級(jí)相差較大，在Fluent軟件中求解流場(chǎng)時(shí)，若網(wǎng)格劃分尺寸過小則計(jì)算機(jī)無法求解，故將葉片及計(jì)算流域原型尺寸縮小到1/1 000，即ln/lm=1 000，下標(biāo)n、m分別代表原型和模型，此時(shí)稱原型和模型具有相似性，即在對(duì)應(yīng)的時(shí)空點(diǎn)，標(biāo)量的特征尺寸大小成比例，矢量速度大小成比例、方向相同[11-12]，并且滿足Stkn=Stkm，斯托克斯數(shù)計(jì)算式如(4)所示。

(4)

式中：Stk為斯托克斯數(shù)，其描述粒子慣性作用與擴(kuò)散作用相對(duì)大小；ρp為粒子密度，kg/m3；dp為粒子直徑，m；μ為空氣動(dòng)力學(xué)黏度，Pa·s；v為風(fēng)速，m/s；l為葉片長度，m。由式(4)可以得到如式(5)所示的規(guī)律。

(5)

為得到更為精確的流場(chǎng)和粒子軌跡，在ICEM CFD軟件中對(duì)模型流域進(jìn)行網(wǎng)格劃分，同時(shí)在樹葉周圍進(jìn)行局部加密處理，并進(jìn)行網(wǎng)格無關(guān)化驗(yàn)證。在Fluent軟件中應(yīng)用式(1)～(3)求解流體運(yùn)動(dòng)方程，邊界條件設(shè)置如下：入口采用速度入口，上、下邊界設(shè)置為對(duì)稱邊界條件，出口設(shè)置為壓力出口，壁面采用無滑移邊界條件[13]，壓力速度耦合計(jì)算采用Simple算法。

1.2 粒子運(yùn)動(dòng)軌跡方程及求解

描述流體運(yùn)動(dòng)的方法有兩種：歐拉法和拉格朗日法。拉格朗日法是指以某一個(gè)流體質(zhì)點(diǎn)的運(yùn)動(dòng)為研究對(duì)象，觀察這一質(zhì)點(diǎn)的運(yùn)動(dòng)參數(shù)(位置、速度、加速度等)隨時(shí)間的變化規(guī)律。本文描述粒子的運(yùn)動(dòng)采用拉格朗日法，將粒子視為離散相，觀察其運(yùn)動(dòng)軌跡。當(dāng)粒徑dp<0.1 μm時(shí)，粒子懸浮在空氣中易受到空氣分子的無規(guī)律撞擊，從而產(chǎn)生隨機(jī)不定向運(yùn)動(dòng)即布朗運(yùn)動(dòng)，此時(shí)主要受布朗隨機(jī)力的影響。隨著粒徑的增大，布朗擴(kuò)散作用逐漸減弱，而慣性作用加強(qiáng)，此時(shí)顆粒物受到流體的斯托克斯阻力?？諝?/p>

中粒子濃度較低時(shí)，可忽略粒子間的相互作用力，則單個(gè)粒子的無規(guī)則Langevin運(yùn)動(dòng)軌跡方程可描述為式(6)[14-16]。

(6)

粒子所受斯托克斯阻力[17]可表示為

(7)

式中：ex、ey分別為x和y軸方向的單位矢量；C為Cunningham修正系數(shù)，C=1+Kn[1.257+0.4exp(-1.1/Kn)]，Kn為粒子的克努森數(shù)。

粒子所受布朗隨機(jī)力[18]可表示為

(8)

式中：ζx、ζy為單位變量獨(dú)立高斯隨機(jī)數(shù)；Δt為時(shí)間步長，取粒子的弛豫時(shí)間，s；T為熱力學(xué)溫度，K；κB為玻爾茲曼常數(shù)，κB=1.38×10-23J/K；D為粒子擴(kuò)散系數(shù)。

1.3 粒子的捕集效率

采用統(tǒng)計(jì)學(xué)方法計(jì)算葉片捕集粒子的效率[19]，如式(9)所示。

(9)

式中：E為模擬捕集效率均值；NG為粒子發(fā)射數(shù)；NC為粒子捕集數(shù)；HG為粒子釋放面高度，m；rl為樹葉的垂直方向投影長度，rl=lsinθ，m。

2 結(jié)果與討論

利用Fluent 16.0軟件對(duì)植物捕集粒子模型進(jìn)行求解。為討論不同機(jī)理下粒子的捕集效率，模擬中氣溶膠粒子需包含布朗擴(kuò)散起主導(dǎo)作用的超細(xì)微粒(dp<0.1 μm)、攔截能力較強(qiáng)的亞微米級(jí)粒子(0.1～1.0 μm)、慣性碰撞起主導(dǎo)作用的微米級(jí)粒子(1～50 μm)、重力沉降起主導(dǎo)作用的超大粒子(>50 μm)[20]。為討論風(fēng)速v和迎風(fēng)角θ對(duì)捕集效率及粒子分布的影響，設(shè)計(jì)參數(shù)如表1所示。

表1 模型主要設(shè)計(jì)參數(shù)Table 1 The main design parameters of the model

2.1 粒子繞葉片表面運(yùn)動(dòng)的軌跡特征

為便于觀察，令葉片位于尺寸為0.010 m×0.006 m平面的中心，如圖2中4條虛線匯聚所成的虛線框所示?？紤]到我國北方地區(qū)近50年來平均風(fēng)速為2.71 m/s[21]，選取vn=3 m/s；同時(shí)為觀察布朗擴(kuò)散、慣性碰撞機(jī)理下粒子的運(yùn)動(dòng)，分別選擇超細(xì)微粒(dp=0.005 μm)、微米級(jí)粒子(dp=5.000 μm)，繪制此時(shí)不同迎風(fēng)角(θ=15°、30°、45°、60°)下的粒子運(yùn)動(dòng)軌跡圖如圖3所示，其中灰色線條為葉片表面繞流流線。

圖2 粒子在葉片表面繞流的研究區(qū)域Fig.2 The main region of particles around the blade surface

由圖3可知，空氣在不同迎風(fēng)角的葉片表面發(fā)生繞流，導(dǎo)致流線發(fā)生偏轉(zhuǎn)，粒子的運(yùn)動(dòng)軌跡也相應(yīng)發(fā)生變化。當(dāng)粒徑為0.005 μm時(shí)，粒子受到的各方向空氣分子撞擊作用不平衡，表現(xiàn)出無規(guī)則的運(yùn)動(dòng)軌跡，即布朗運(yùn)動(dòng)，此時(shí)擴(kuò)散機(jī)理占據(jù)主導(dǎo)地位；隨著粒徑的增大(dp=5.000 μm)，粒子運(yùn)動(dòng)變得較為規(guī)律，當(dāng)氣流經(jīng)過葉片表面發(fā)生偏轉(zhuǎn)時(shí)，粒子由于慣性作用的增強(qiáng)，不能完全跟隨空氣分子及時(shí)避開葉片，導(dǎo)致發(fā)生碰撞的粒子被葉片表面所捕集[22]，遵循慣性碰撞捕集機(jī)理。

圖3 植物葉片繞流流場(chǎng)分布及粒子軌跡(vn=3 m/s)Fig.3 Orbital flow field distribution and particle trajectory of plant leaves (vn=3 m/s)

2.2 風(fēng)速對(duì)葉片捕集粒子的影響

不同迎風(fēng)角下，捕集效率隨風(fēng)速表現(xiàn)出一致的規(guī)律性，選擇迎風(fēng)角θ=60°，繪制不同風(fēng)速下葉片的捕集效率，如圖4所示。由圖4可知，捕集效率隨粒徑的增大呈先減小再增大的趨勢(shì)。對(duì)于粒徑較小(dp<0.1 μm)的粒子，捕集效率隨風(fēng)速的增大而減小，下降趨勢(shì)由急速變?yōu)槠骄?；?duì)于粒徑為0.1～1.0 μm的粒子，捕集效率無明顯變化趨勢(shì)；對(duì)于粒徑dp>1.0 μm的粒子，捕集效率隨風(fēng)速的增大而增大。由圖4(b)可知，隨著風(fēng)速的增大，捕集效率由下降轉(zhuǎn)變?yōu)樯仙鶎?duì)應(yīng)的粒徑逐漸減小。

圖4 捕集效率與粒徑和風(fēng)速的關(guān)系Fig.4 Relationship between capture efficiency and particle sizes/wind speeds

2.3 迎風(fēng)角對(duì)捕集效率的影響

現(xiàn)實(shí)生活中植物葉片通常呈多種不規(guī)律角度，迎風(fēng)角對(duì)捕集效率的影響尚待討論，通過微觀尺度定性分析其對(duì)捕集效率的影響，可為日后研究整株植物不同迎風(fēng)角下的捕集效率奠定基礎(chǔ)。選取風(fēng)速vn=5 m/s，將不同迎風(fēng)角下葉片的捕集效率繪制成圖5。由圖5可知，不同迎風(fēng)角時(shí)的捕集效率之間無顯著差異。在迎風(fēng)角低于90°條件下，當(dāng)粒徑dp<1.0 μm時(shí)，布朗擴(kuò)散捕集效率隨迎風(fēng)角的增大而減小，隨粒徑的進(jìn)一步增大，捕集效率也隨之增大。

圖5 捕集效率與粒徑和迎風(fēng)角的關(guān)系Fig.5 Relationship between capture efficiency and particle sizes/windward angles

2.4 捕集效率公式擬合

根據(jù)粒子沉積機(jī)理，分段闡述捕集效率，即布朗擴(kuò)散效率EB和慣性碰撞效率Eimp。布朗運(yùn)動(dòng)中采用擴(kuò)散系數(shù)刻畫粒子的運(yùn)動(dòng)快慢[23]，引入施密特?cái)?shù)Sc描述粒子的動(dòng)量擴(kuò)散與質(zhì)量擴(kuò)散，如式(10)[23]所示。隨著粒徑的增大，慣性作用加強(qiáng)，采用斯托克斯數(shù)Stk描述粒子的運(yùn)動(dòng)，如式(4)所示。Stk越大，顆粒與氣體流線的重合率越小，即顆粒抵抗跟隨氣體流線流動(dòng)的能力增強(qiáng)[24]。

(10)

式中：a為20 ℃時(shí)的空氣運(yùn)動(dòng)黏度，m2/s；D為布朗擴(kuò)散系數(shù)，m2/s。

2.4.1 布朗擴(kuò)散效率的擬合

采用Origin 9.0軟件將不同風(fēng)速、粒徑、迎風(fēng)角下布朗擴(kuò)散效率EB的數(shù)據(jù)擬合成施密特?cái)?shù)Sc、迎風(fēng)角θ及風(fēng)速v的公式，如式(11)所示。

EB=1.219sinθ-0.139Sc-0.547v-0.467

(11)

式中：θ為弧度制迎風(fēng)角。

將CFD模擬得到的捕集效率數(shù)據(jù)與擬合得到的布朗擴(kuò)散效率EB繪制成圖6，圖中散點(diǎn)為模擬所得數(shù)據(jù)，曲線為擬合結(jié)果。由圖6可知：6條彼此不重合的曲線代表不同風(fēng)速下的捕集效率，曲線斜率隨風(fēng)速的增大而減??；在相同風(fēng)速下，不同迎風(fēng)角的捕集效率曲線出現(xiàn)一定程度的重合，進(jìn)一步表明迎風(fēng)角對(duì)布朗擴(kuò)散捕集效率的影響甚微。此外，當(dāng)Sc<105即粒徑dp<0.3 μm時(shí)，擬合結(jié)果與模擬結(jié)果吻合度較高，相關(guān)系數(shù)r=0.99，對(duì)于小粒徑的粒子，捕集效率隨Sc的增大呈規(guī)律性遞減，隨風(fēng)速的增大而減小。

圖6 布朗擴(kuò)散捕集效率擬合與模擬數(shù)據(jù)驗(yàn)證Fig.6 Brownian diffusion capture efficiency fitting and simulation data verification

2.4.2 慣性碰撞效率的擬合

微米級(jí)粒子在被葉片捕集時(shí)，慣性碰撞起主導(dǎo)作用，故選擇粒徑尺寸dp>1 μm時(shí)模擬得到的捕集效率，與斯托克斯數(shù)、風(fēng)速、迎風(fēng)角進(jìn)行擬合，得到慣性捕集效率Eimp的表達(dá)式，如式(12)所示。

(12)

將模擬所得數(shù)據(jù)與擬合公式所得計(jì)算結(jié)果繪制成圖7。由圖7可知，當(dāng)Stk>10-5時(shí)，即對(duì)應(yīng)粒徑(dp>1 μm)擬合所得結(jié)果與模擬數(shù)據(jù)的相關(guān)性較好，r=0.97。由圖7可知：在銳角范圍內(nèi)，迎風(fēng)角的增大促使捕集效率Eimp增大，這可能是因?yàn)樵龃笥L(fēng)角度將使得葉片作為障礙物的捕集面積增大；對(duì)于微米級(jí)粒子，Eimp隨著Stk的增加而增大。

圖7 慣性碰撞捕集效率擬合與模擬數(shù)據(jù)驗(yàn)證Fig.7 Inertia collision capture efficiency fitting and simulation data verification

2.4.3 總效率的擬合結(jié)果

依據(jù)經(jīng)典過濾理論[25]，定義植物葉片捕集粒子的總效率E=EB+Eimp，如式(13)所示。

E=1.219sinθ-0.139Sc-0.547v-0.467+

(13)

將擬合結(jié)果與模擬所得結(jié)果繪制成圖8。由圖8可知，粒徑為0.1～1.0 μm時(shí)擬合結(jié)果與模擬結(jié)果差異較大，這可能是因?yàn)閿M合時(shí)忽略了攔截作用所造成的粒子捕集?？傮w而言，擬合結(jié)果與模擬結(jié)果的吻合度較高。

圖8 總效率擬合與模擬數(shù)據(jù)驗(yàn)證Fig.8 Verification between the total efficiency of the fit and the simulated data

2.5 植物單葉片捕集總效率預(yù)測(cè)公式的驗(yàn)證

文獻(xiàn)中布朗擴(kuò)散、慣性碰撞及總效率的參數(shù)化表達(dá)式如式(14)～(16)所示。

(14)

(15)

(16)

式中：Ea為文獻(xiàn)總效率；EBa為文獻(xiàn)[6]的布朗擴(kuò)散效率；Eimpa為文獻(xiàn)[7]的慣性碰撞效率；cv/cd為黏性阻力比，取0.3；γ為經(jīng)驗(yàn)常數(shù)，植物表面取值為0.67；α為經(jīng)驗(yàn)常數(shù)，與植被冠層種類有關(guān)，取混合闊葉樹和針葉樹所對(duì)應(yīng)的參數(shù)值0.8。

文獻(xiàn)[6-7]中研究條件均是風(fēng)向垂直于葉片，故式(11)～(13)中迎風(fēng)角度取90°，其余參數(shù)取值如表1所示，將計(jì)算得到的布朗擴(kuò)散、慣性碰撞以及總捕集效率分別與文獻(xiàn)結(jié)果進(jìn)行對(duì)比，如圖9～11所示。由圖9可知，對(duì)于小粒徑(dp<1.0 μm對(duì)應(yīng)Sc<105)粒子而言，擬合結(jié)果與文獻(xiàn)結(jié)果的變化趨勢(shì)基本一致，本文EB的擬合結(jié)果為多條斜率的曲線，而文獻(xiàn)[6]僅有1條曲線。這是由于本文將風(fēng)速這一影響因素考慮進(jìn)去，由擬合公式可知，斜率越大風(fēng)速越小。由圖10可知，對(duì)于大粒徑(dp>1.0 μm對(duì)應(yīng)Stk>10-5)粒子而言，Eimp擬合結(jié)果與文獻(xiàn)值的整體變化趨勢(shì)一致且吻合度較高，說明本文提出的方法具有可行性。由圖11可知，總效率E的擬合與文獻(xiàn)結(jié)果變化趨勢(shì)一致，但在粒徑處于0.1～1.0 μm時(shí)，擬合結(jié)果的值存在過高估計(jì)的現(xiàn)象，這可能時(shí)由于本文沒有考慮粒子的攔截作用。

圖9 布朗擴(kuò)散效率的擬合結(jié)果與文獻(xiàn)結(jié)果Fig.9 Fitting results and literature results of Brownian diffusion efficiency

圖10 慣性碰撞效率的擬合結(jié)果與文獻(xiàn)結(jié)果Fig.10 Fitting results and literature results of inertia collision efficiency

圖11 總效率的擬合結(jié)果與文獻(xiàn)結(jié)果Fig.11 Fitting results and literature results of total efficiency

3 結(jié) 語

本文針對(duì)植物葉片捕集粒子的物理現(xiàn)象，應(yīng)用物理模型相似理論及數(shù)值分析法求解植物葉片表面附近的流場(chǎng)，同時(shí)考慮粒子布朗擴(kuò)散和慣性碰撞的聯(lián)合作用，應(yīng)用拉格朗日法追蹤粒子運(yùn)動(dòng)軌跡并求解植物葉片表面粒子捕集行為。分析計(jì)算得到捕集效率，討論了粒徑dp、迎風(fēng)角度θ、風(fēng)速v對(duì)捕集效率的影響，同時(shí)對(duì)效率與施密特?cái)?shù)Sc、斯托克斯數(shù)Stk進(jìn)行擬合得到效率的預(yù)測(cè)式，通過與文獻(xiàn)進(jìn)行對(duì)比，得到以下結(jié)論：

(1) 隨著粒徑的增大，捕集效率呈先減小后增大的趨勢(shì)，這與葉片捕集粒子的機(jī)理相關(guān)。粒徑較小時(shí)粒子的布朗擴(kuò)散起主導(dǎo)作用，隨著粒徑的增大，慣性作用增強(qiáng)，慣性碰撞起主導(dǎo)作用。

(2) 對(duì)于dp<1.0 μm的粒子，迎風(fēng)角度對(duì)捕集效率的影響很小，隨迎風(fēng)角的增大，捕集效率呈下降趨勢(shì)；對(duì)于dp>1.0 μm的粒子，迎風(fēng)角度θ的影響顯而易見，隨著迎風(fēng)角的增大，葉片捕集的粒子數(shù)也呈現(xiàn)增加趨勢(shì)。

(3) 當(dāng)粒徑dp<0.1 μm時(shí)，捕集效率隨風(fēng)速的增大而降低；當(dāng)粒徑在0.1～1.0 μm時(shí)，捕集效率隨風(fēng)速變化不明顯；當(dāng)dp>1.0 μm時(shí)，捕集效率與風(fēng)速呈正相關(guān)。

(4) 根據(jù)不同的捕集機(jī)理分別擬合出布朗擴(kuò)散、慣性碰撞作用下捕集效率公式，應(yīng)用疊加原理得到總效率預(yù)測(cè)公式。本文中捕集效率隨粒徑的整體變化趨勢(shì)與文獻(xiàn)研究結(jié)果較吻合，但在粒徑為0.1～1.0 μm時(shí)，本文應(yīng)用預(yù)測(cè)公式計(jì)算所得的總效率偏大。