nvPM采樣輸運(yùn)系統(tǒng)損失因子計(jì)算與改進(jìn)研究

楊曉軍，靳曉宇，胡英琦，劉智剛

(中國民航大學(xué)航空工程學(xué)院，天津 300300)

0 引 言

隨著經(jīng)濟(jì)的高速發(fā)展，城市化、工業(yè)化規(guī)模的擴(kuò)大，空氣中各種環(huán)境污染物的排放量迅速增加，細(xì)顆粒物污染尤為嚴(yán)重[1]。航空發(fā)動機(jī)排放的非揮發(fā)性顆粒物絕大部分是由燃燒室產(chǎn)生的碳質(zhì)顆粒、發(fā)動機(jī)進(jìn)氣道的顆粒和金屬顆粒構(gòu)成，也是大氣中煙塵顆粒污染的重要來源之一。與其他大多數(shù)的燃燒源相比，其排放的顆粒物不僅直徑遠(yuǎn)小于10 μm（大多處于納米級別），而且由于飛機(jī)巡航狀態(tài)時(shí)排放的顆粒物，既會進(jìn)入大氣邊界層，也將注入對流層上部和平流層下部[1-3]，根據(jù)聯(lián)合國政府間氣候變化專門委員會（IPCC）[4]對其影響評估，可知航空顆粒物排放被認(rèn)為是云形變化和全球氣候變化的潛在強(qiáng)力影響因素。同時(shí)納米級的微粒還會進(jìn)入人類的呼吸系統(tǒng)，嚴(yán)重危害人體健康。

為此，能夠準(zhǔn)確地監(jiān)測發(fā)動機(jī)排放情況是十分必要的。國際民航組織[5]（ICAO）頒布了相關(guān)測量建議標(biāo)準(zhǔn)，為不受環(huán)境干擾，通過一根不超過35 m長的采樣管，把采集的尾氣輸運(yùn)到測量儀器處，但在采樣管中，粒子會發(fā)生大量損失使得測量結(jié)果小于真實(shí)值。由此美國汽車工程師協(xié)會（SAE）E31委員會頒布了相關(guān)測量校正標(biāo)準(zhǔn)（AIR6504和ARP6481[6-7]），允許由測量儀器在飛機(jī)發(fā)動機(jī)排放平面處所測得的nvPM質(zhì)量和數(shù)量濃度數(shù)據(jù)，估算nvPM質(zhì)量和數(shù)量損失校正因子，并提供了詳細(xì)的計(jì)算校正因子流程。

在SAE提供的計(jì)算流程中，假設(shè)發(fā)動機(jī)排放平面處的nvPM數(shù)量濃度按對數(shù)正態(tài)分布函數(shù)分布，使用假設(shè)的平均粒徑和被設(shè)定為固定值的對數(shù)正態(tài)分布函數(shù)的標(biāo)準(zhǔn)偏差和顆粒密度進(jìn)行估算。Abegglen等[8]在研究粒徑和發(fā)動機(jī)推力對顆粒密度的影響時(shí)，發(fā)現(xiàn)密度隨推力和粒徑變化而變化，推力與粒徑呈正相關(guān)關(guān)系；Lobo等[9]在比較3種采樣測量參考系統(tǒng)時(shí)，測量得出了CFM56-7B26/3型發(fā)動機(jī)的質(zhì)量和數(shù)量濃度排放指數(shù)，與其對應(yīng)的分布函數(shù)，且發(fā)現(xiàn)平均粒徑與推力呈正相關(guān)趨勢；Delhaye等[10]進(jìn)行了MERMOSE項(xiàng)目，發(fā)現(xiàn)平均粒徑隨推力增大而增大，從17～55 nm不等，且在推力為85%時(shí)檢測到最大數(shù)量濃度排放；Lobo等[11]在不同的LTO工況下發(fā)現(xiàn)在低推力下測量得到的nvPM質(zhì)量濃度跳動較大，驗(yàn)證得出測量體系在低推力下測量質(zhì)量濃度不穩(wěn)定。大多數(shù)研究可以表明損失校正計(jì)算流程存在假設(shè)數(shù)據(jù)與實(shí)際不符和測量結(jié)果不穩(wěn)定等現(xiàn)象，但是對于改進(jìn)校正計(jì)算方法所涉及的相關(guān)研究較少。

由此本研究分析并改進(jìn)了計(jì)算校正因子的方法，使用新的計(jì)算校正工具，結(jié)合所掌握的 CFM56-7B26/3型發(fā)動機(jī)的相關(guān)數(shù)據(jù)，對新的校正工具進(jìn)行了驗(yàn)證與結(jié)果分析。隨后進(jìn)一步分析該型發(fā)動機(jī)數(shù)據(jù)，對計(jì)算工具進(jìn)行了再度優(yōu)化，使其能夠僅輸入發(fā)動機(jī)推力即可得到結(jié)果，令計(jì)算校正工具能更加方便使用。

1 數(shù)據(jù)和方法

1.1 采樣標(biāo)準(zhǔn)方法

在ICAO和SAE頒布的標(biāo)準(zhǔn)采樣測量和校正實(shí)驗(yàn)流程中，使用在發(fā)動機(jī)排放平面處呈十字釘耙狀布置的不銹鋼合金探針引出排放尾氣。為保證尾氣中的顆粒不受環(huán)境干擾，使其進(jìn)入一個(gè)不超過35 m長的采樣管，管中連接稀釋分離儀器以保證流量和濃度在管末布置的采樣儀器測量量程范圍內(nèi)，最后由采樣儀器測得結(jié)果。在整個(gè)過程中會造成顆粒大量損失，主要有輸運(yùn)、儀器、收集損失。

1.1.1 輸運(yùn)損失

輸運(yùn)損失是指尾氣中的顆粒在采樣管中發(fā)生的損失[6]，主要有擴(kuò)散、彎曲、慣性、靜電、熱泳5種損失機(jī)制，每種機(jī)制都經(jīng)過大量實(shí)驗(yàn)和理論分析驗(yàn)證，推導(dǎo)得出粒徑與穿透概率的關(guān)系式[7]。其中若顆粒100%穿過，則穿透概率為1。關(guān)系式主要變量為粒徑，也包含各種測量條件，如采樣管相關(guān)參數(shù)、載氣和環(huán)境條件參數(shù)。校正標(biāo)準(zhǔn)主要通過5種機(jī)制疊加相乘來估算粒子通過采樣管的穿透概率（如圖1中P Num和P Mass），也就是穿透率?？偞┩嘎孰S粒徑變化趨勢如圖1所示。

1.1.2 儀器損失

儀器損失主要是粒子穿過測量儀器時(shí)所發(fā)生的損失，標(biāo)準(zhǔn)中詳細(xì)描述了管路儀器布置，其中主要儀器有稀釋器（Diluter）、凝結(jié)粒子計(jì)數(shù)器（condensation particle counter，CPC）、旋風(fēng)分離器（cyclone separator，Cyc）、揮發(fā)粒子去除器（volatile particle remover，VPR）。標(biāo)準(zhǔn)內(nèi)詳細(xì)給出了儀器所對應(yīng)的粒徑與穿透率關(guān)系式，如圖1所示。

圖1 粒徑與穿透率關(guān)系圖

1.1.3 收集損失

收集損失主要指收集探針處發(fā)生的損失。因考慮因素過多，成因復(fù)雜，所以SAE標(biāo)準(zhǔn)中雖然提及該部分的熱泳損失穿透率方程，但并未在計(jì)算中考慮和應(yīng)用。如圖1中Thermo曲線所示。

1.2 損失計(jì)算方法

SAE測量校正標(biāo)準(zhǔn)解釋了除收集部件處損失之外，估算nvPM質(zhì)量和數(shù)量濃度損失校正因子的迭代計(jì)算方法。詳細(xì)迭代計(jì)算流程如圖2所示。

圖2 系統(tǒng)損失校正因子的迭代計(jì)算流程圖

估算校正因子之前，輸入所有采樣管段的條件參數(shù)和各個(gè)儀器調(diào)試校驗(yàn)后的穿透率結(jié)果，最后輸入測量儀器顯示的質(zhì)量與數(shù)量濃度。按上述方法開始運(yùn)行迭代計(jì)算，估算得出排放平面處的平均粒徑和質(zhì)量、數(shù)量濃度校正因子。

在迭代運(yùn)算中，假設(shè)發(fā)動機(jī)排放平面處的nvPM數(shù)量濃度按對數(shù)正態(tài)分布函數(shù)分布，且平均粒徑（Dmg）先假設(shè)為6 nm，標(biāo)準(zhǔn)偏差為固定值1.8，顆粒的密度為常數(shù)1 000 kg/m3。運(yùn)用上述假設(shè)和各粒徑在不同損失機(jī)制中的穿透率關(guān)系式，估算出各粒徑在發(fā)動機(jī)排放平面出口處的數(shù)量濃度，由密度推出各粒徑的質(zhì)量濃度，各粒徑疊加算出總質(zhì)量和數(shù)量濃度，認(rèn)定輸入儀器測量的質(zhì)量和數(shù)量濃度比率（以下簡稱質(zhì)數(shù)比）與出口平面一致，用估算的質(zhì)數(shù)比與測量的質(zhì)數(shù)比進(jìn)行相對誤差運(yùn)算，當(dāng)兩者的卡方誤差Δ不滿足小于 ε（ε=1×10-9）的條件時(shí)，改變上述假設(shè)輸入的平均粒徑，直到滿足卡方誤差小于ε，認(rèn)定最后求得的Dmg為實(shí)際排放平面對數(shù)分布的期望值μ （即幾何平均粒徑），用平均粒徑估算出發(fā)動機(jī)排放平面處的質(zhì)量和數(shù)量濃度為實(shí)際濃度，應(yīng)用損失穿透率計(jì)算損失后的各濃度，兩者相比即為質(zhì)量校正因子和數(shù)量校正因子。

1.3 航空發(fā)動機(jī)排放數(shù)據(jù)

根據(jù)目前所掌握的多種類型的發(fā)動機(jī)排放數(shù)據(jù)，大體上顯示出發(fā)動機(jī)的平均粒徑與標(biāo)準(zhǔn)差有著正相關(guān)關(guān)系（如圖3所示）。由此猜想，原有估算方式中假設(shè)標(biāo)準(zhǔn)偏差為常數(shù)與實(shí)際情況不相符。為改進(jìn)計(jì)算校正方法，選取CFM56-7B26/3型發(fā)動機(jī)數(shù)據(jù)進(jìn)行研究，主要包含12組該型發(fā)動機(jī)不同推力下的質(zhì)數(shù)比、平均粒徑等，按推力從小到大的順序進(jìn)行排列，同時(shí)在排列后發(fā)現(xiàn)恰好后7組的質(zhì)數(shù)比趨勢也是從小到大的。而后在討論中，經(jīng)過分析排放數(shù)據(jù)與推力相關(guān)聯(lián)系，將計(jì)算工具進(jìn)行了再度優(yōu)化，令其使用更加便捷高效。

圖3 多種發(fā)動機(jī)的標(biāo)準(zhǔn)偏差與平均粒徑關(guān)系圖

1.4 修改損失計(jì)算方法

如圖4所示，主要改進(jìn)為直接輸入測量得出的平均粒徑Dmg和測量質(zhì)量、數(shù)量濃度，再由迭代計(jì)算找出最合適的標(biāo)準(zhǔn)偏差。使用平均粒徑和估算出的標(biāo)準(zhǔn)偏差進(jìn)行排放平面處數(shù)量濃度分布函數(shù)估算，再求得質(zhì)量和數(shù)量校正因子。按照SAE所提供的采樣測量校正標(biāo)準(zhǔn)、損失穿透率方程和改進(jìn)的計(jì)算流程，運(yùn)用Matlab軟件，對改進(jìn)的流程進(jìn)行編程。

圖4 改進(jìn)的迭代運(yùn)算流程圖

2 結(jié)果與分析

2.1 工具驗(yàn)證

工具驗(yàn)證數(shù)據(jù)采用在原有損失計(jì)算流程中提供的10組質(zhì)量濃度和數(shù)量濃度，即10組不同的質(zhì)數(shù)比。為了便于觀察，把10組數(shù)據(jù)按照質(zhì)數(shù)比從小到大的順序排列。由于新方法需要幾何平均粒徑作為輸入?yún)?shù)，所以先使用原有計(jì)算方法進(jìn)行計(jì)算，得出平均粒徑后，再輸入到新計(jì)算工具中。若得出標(biāo)準(zhǔn)偏差和結(jié)果與原方法一致，即可驗(yàn)證該工具的適用性。

為了使對比結(jié)果更為直觀，采用相對偏差進(jìn)行觀察，相對偏差定義如下：

式中：POld——原計(jì)算方法所得結(jié)果；

PNew——改進(jìn)的新方法所得結(jié)果。

對比驗(yàn)證結(jié)果如表1所示。由表可得，除第1組數(shù)據(jù)之外，所得σ、質(zhì)量和數(shù)量校正因子的相對誤差均小于0.01%，可驗(yàn)證該工具有良好的準(zhǔn)確性?？ǚ秸`差所代表的是估算值與輸入真實(shí)測量值之間的相對誤差平方?？ǚ秸`差越小，所估算的值越接近真實(shí)測量值。在原計(jì)算方法中當(dāng)δ＜ε（ε=1×10-9）時(shí)，才能說明該計(jì)算結(jié)果有效。

表1 對比驗(yàn)證及誤差1）

因第1組數(shù)據(jù)在原計(jì)算工具中也出現(xiàn)問題，所以對該組數(shù)據(jù)存疑。觀察所得3個(gè)參數(shù)的結(jié)果可發(fā)現(xiàn)，在卡方誤差比原計(jì)算方法大的情況下，所得結(jié)果仍與原結(jié)果一致。推測該算法靈敏度更高，估算質(zhì)數(shù)比與測量輸入質(zhì)數(shù)比在誤差比原方法大的范圍內(nèi)，可得出相同結(jié)果。例如在第7組中新方法的卡方誤差大于1×10-9，但其結(jié)果的相對偏差在1×10-6左右，與原結(jié)果幾乎相同。由此，在不影響結(jié)果的情況下，降低判別新方法結(jié)果有效的條件，設(shè)定ε升高一個(gè)數(shù)量級，即ε=1×10-8。當(dāng) ＜ε時(shí)，新方法即可得出有效結(jié)果。由此可證明新計(jì)算工具通過驗(yàn)證。

2.2 結(jié)果對比

使用1.3節(jié)的12組CFM56-7B26/3型發(fā)動機(jī)數(shù)據(jù)，分別使用原計(jì)算方法和新計(jì)算工具進(jìn)行運(yùn)算，所得結(jié)果對比如下。

2.2.1 幾何平均粒徑

此處對平均粒徑的比較中，新計(jì)算工具使用的是真實(shí)粒徑（如表2所示），原計(jì)算方式為估算結(jié)果。如圖5所示，黑色實(shí)點(diǎn)為該型發(fā)動機(jī)實(shí)際平均粒徑。兩者趨勢相同只是原方法估算的平均粒徑偏小，最大絕對誤差達(dá)到15.3 nm，最小為5.3 nm。實(shí)際數(shù)據(jù)在推力為7%時(shí)為最小粒徑14.35 nm，且與3%處粒徑15.21 nm相差不大，所以在3、4組處有下凹趨勢，而原方法在第1組與第6組處出現(xiàn)最小值6 nm。說明原方法估算的平均粒徑與實(shí)際有較大差異。通過數(shù)據(jù)分析，原因?yàn)檫@兩組質(zhì)數(shù)比較小，原方法所體現(xiàn)的是平均粒徑與質(zhì)數(shù)比呈正相關(guān)關(guān)系，如圖6所示，但與實(shí)際情況不符。由此可猜想數(shù)據(jù)之間確實(shí)存在某種聯(lián)系，在下章將進(jìn)行詳細(xì)討論分析。

表2 CFM56-7B26/3型發(fā)動機(jī)不同推力下的幾何平均粒徑

圖5 幾何平均粒徑對比

圖6 平均粒徑與質(zhì)數(shù)比關(guān)系圖

2.2.2 nvPM數(shù)量濃度系統(tǒng)損失校正因子

數(shù)量濃度損失校正因子主要與平均粒徑相關(guān)，如圖7所示。因上述計(jì)算的平均粒徑兩者趨勢相同，所以兩者數(shù)量濃度損失校正因子總體趨勢也相同，但受原方法所估算平均粒徑偏小的影響，新估算的數(shù)量損失校正因子總體比原方法計(jì)算結(jié)果小，結(jié)果如圖8所示。在第1組和第6組處絕對誤差較大，最大為6.81，在第12組處最小為0.64，整體誤差趨勢是隨質(zhì)數(shù)比增大而減小的，可以說明工具在較大的質(zhì)數(shù)比處，結(jié)果更加穩(wěn)定。

圖7 數(shù)量校正因子受平均粒徑影響趨勢

圖8 nvPM數(shù)量濃度系統(tǒng)損失校正因子對比

2.2.3 nvPM質(zhì)量濃度系統(tǒng)損失校正因子

質(zhì)量濃度校正因子受平均粒徑和標(biāo)準(zhǔn)偏差影響不大，所以兩者的質(zhì)量濃度校正因子差異相對較小。在后幾組質(zhì)數(shù)比和發(fā)動機(jī)推力較大的情況下，部分校正因子重合，結(jié)果如圖9所示。在第6組處絕對誤差最大為0.14，在第12組處最小為0.007。說明該工具計(jì)算質(zhì)量校正因子在質(zhì)數(shù)比較大的情況下較為穩(wěn)定。

圖9 nvPM質(zhì)量濃度系統(tǒng)損失校正因子對比

3 再優(yōu)化與討論

3.1 數(shù)據(jù)分析與工具優(yōu)化

將掌握的12組CFM56-7B26/3型發(fā)動機(jī)輸入數(shù)據(jù)與更多文獻(xiàn)數(shù)據(jù)進(jìn)行對比分析，發(fā)現(xiàn)平均粒徑與推力呈現(xiàn)出正相關(guān)趨勢，如圖10所示。由此猜想，各數(shù)據(jù)與推力之間存在著某種特定關(guān)系。

圖10 掌握數(shù)據(jù)與Lobo[9]數(shù)據(jù)對比

經(jīng)數(shù)據(jù)分析后得到質(zhì)數(shù)比隨發(fā)動機(jī)推力增長曲線和平均粒徑與推力線性關(guān)系圖如圖11所示。

圖11 推力與質(zhì)數(shù)比及平均粒徑關(guān)系圖

運(yùn)用軟件對數(shù)據(jù)進(jìn)行分析擬合，可以得出兩個(gè)較好的擬合曲線，r2分別為0.958 42和0.986 91，推力與質(zhì)數(shù)比的殘差平方和較小，趨近于0，推力與平均粒徑的皮爾遜相關(guān)系數(shù)為0.99。

推力與質(zhì)數(shù)比關(guān)系式：

推力與平均粒徑關(guān)系式：

以推力為變量，代入式（2）、式（3）可得所對應(yīng)的質(zhì)數(shù)比和平均粒徑，將所得結(jié)果代入改進(jìn)的計(jì)算流程程序可得到一個(gè)只有推力作為單一變量的計(jì)算工具。由于只導(dǎo)入了CFM56-7B26/3型發(fā)動機(jī)的相關(guān)關(guān)系式，所以該計(jì)算工具只對應(yīng)于該型發(fā)動機(jī)的推力輸入和校正因子結(jié)果輸出。

3.2 結(jié)果分析

所得數(shù)量校正因子與推力的趨勢和原方法一致，且最大誤差為3.05，最小為0.58。與使用估算平均粒徑的原方法相比，新方法由于代入了真實(shí)測量的平均粒徑，極大地改變了在低推力和低質(zhì)數(shù)比下的數(shù)量校正因子。如圖12、圖13所示，隨著推力的增加，質(zhì)量和數(shù)量校正因子都呈下降趨勢，因質(zhì)數(shù)比和平均粒徑隨該發(fā)動機(jī)的推力增大而增大，導(dǎo)致較大推力下的尾氣中主要以粒徑較大的大顆粒為主，對某些損失機(jī)制來說，顆粒的粒徑較大意味著更高的穿透率，所以損失因子逐漸變小。

圖12 推力與數(shù)量校正因子關(guān)系圖

圖13 推力與質(zhì)量校正因子關(guān)系圖

數(shù)量濃度受粒徑變化導(dǎo)致的損失機(jī)制影響較大。因原方法使用的估算粒徑比實(shí)際小，所以新計(jì)算工具得到的數(shù)量損失因子整體比較小。在排放尾氣中，還可能存在多個(gè)小粒徑粒子粘黏聚合成較大粒子的情況，也會顯著降低數(shù)量濃度。

平均粒徑對質(zhì)量損失校正因子影響很小，且無論顆粒聚合與否，其總體質(zhì)量都不會改變，所以其變化量較小，最大絕對誤差為0.03，最小為0.006，與原方法計(jì)算結(jié)果趨勢相似。因此，此方法運(yùn)用更多的實(shí)際排放數(shù)據(jù)進(jìn)行運(yùn)算，能得到良好的效果，并且能與發(fā)動機(jī)推力相聯(lián)系，使其計(jì)算更加簡單便捷。

3.3 適用性討論

在表1中，第1組的數(shù)據(jù)經(jīng)兩種算法誤差都比較大。使用新計(jì)算方法反推發(fā)動機(jī)推力，結(jié)果發(fā)現(xiàn)在反推過程中由于該組質(zhì)數(shù)比過小，只有1.054 82×10-5，導(dǎo)致方程無解，即沒有過小的推力支持此項(xiàng)質(zhì)數(shù)比。由此可推測，兩種算法都不適合質(zhì)數(shù)比過小的情況，或者說該質(zhì)數(shù)比與CFM56-7B26/3型發(fā)動機(jī)不匹配，也可能該數(shù)據(jù)不真實(shí)存在。

4 結(jié)束語

1）掌握的多種發(fā)動機(jī)平均粒徑與標(biāo)準(zhǔn)偏差數(shù)據(jù)表明兩者呈正相關(guān)趨勢，由此原方法定義固定標(biāo)準(zhǔn)偏差是不合適的。而后對計(jì)算工具進(jìn)行了改進(jìn)。由于新計(jì)算工具使用了真實(shí)測量的平均粒徑數(shù)據(jù)，所以數(shù)量濃度校正因子受影響導(dǎo)致明顯差異。因真實(shí)平均粒徑總體都大于原方法所估算的平均粒徑，所以數(shù)量校正因子隨推力變化趨勢與原方法一致，但總體小于原方法的結(jié)果，最大絕對誤差為6.81，最小為0.64，且在推力為3%（慢車狀態(tài)）處校正因子最高，差異也較大。

2）新計(jì)算工具得出的質(zhì)量濃度校正因子受粒徑影響不大，并與原方法趨勢一致，與推力呈負(fù)相關(guān)，校正因子最高的地方同樣是在推力為3%處。與原方法結(jié)果對比，能夠得到良好的結(jié)果，絕對誤差最大為0.14，最小為0.007。通過結(jié)果的推力關(guān)系圖表明，在低推力下計(jì)算校正因子是測量校正工作的難點(diǎn)與重點(diǎn)。

3）測量是針對某一型發(fā)動機(jī)進(jìn)行的，所以針對CFM56-7B26/3型發(fā)動機(jī)數(shù)據(jù)進(jìn)一步分析，得出其推力與質(zhì)數(shù)比和平均粒徑的擬合函數(shù)曲線，兩者r2都大于0.95，由此得到該發(fā)動機(jī)下推力與校正因子的關(guān)系圖。為使計(jì)算更加簡便快捷，分發(fā)動機(jī)型建模計(jì)算與分析是合理且必要的。

4）CFM56-7B26/3型發(fā)動機(jī)是很典型的一種發(fā)動機(jī)，對其數(shù)據(jù)進(jìn)行分析研究與校正因子計(jì)算流程改進(jìn)，能夠?yàn)榘l(fā)動機(jī)排放測量損失校正團(tuán)隊(duì)提供一條與實(shí)際聯(lián)系更密切、計(jì)算更便捷的新思路。但針對其他型號發(fā)動機(jī)使用該方法計(jì)算時(shí)，應(yīng)先進(jìn)行一次全面的測試分析，得出準(zhǔn)確的推力與質(zhì)數(shù)比及平均粒徑之間的關(guān)系方程后，才能更便捷的使用該計(jì)算工具。