nvPM采樣輸運系統損失因子計算與改進研究

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

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

0 引 言

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

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

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

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

1 數據和方法

1.1 采樣標準方法

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

1.1.1 輸運損失

輸運損失是指尾氣中的顆粒在采樣管中發生的損失[6]，主要有擴散、彎曲、慣性、靜電、熱泳5種損失機制，每種機制都經過大量實驗和理論分析驗證，推導得出粒徑與穿透概率的關系式[7]。其中若顆粒100%穿過，則穿透概率為1。關系式主要變量為粒徑，也包含各種測量條件，如采樣管相關參數、載氣和環境條件參數。校正標準主要通過5種機制疊加相乘來估算粒子通過采樣管的穿透概率（如圖1中P Num和P Mass），也就是穿透率。總穿透率隨粒徑變化趨勢如圖1所示。

1.1.2 儀器損失

儀器損失主要是粒子穿過測量儀器時所發生的損失，標準中詳細描述了管路儀器布置，其中主要儀器有稀釋器（Diluter）、凝結粒子計數器（condensation particle counter，CPC）、旋風分離器（cyclone separator，Cyc）、揮發粒子去除器（volatile particle remover，VPR）。標準內詳細給出了儀器所對應的粒徑與穿透率關系式，如圖1所示。

圖1 粒徑與穿透率關系圖

1.1.3 收集損失

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

1.2 損失計算方法

SAE測量校正標準解釋了除收集部件處損失之外，估算nvPM質量和數量濃度損失校正因子的迭代計算方法。詳細迭代計算流程如圖2所示。

圖2 系統損失校正因子的迭代計算流程圖

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

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

1.3 航空發動機排放數據

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

圖3 多種發動機的標準偏差與平均粒徑關系圖

1.4 修改損失計算方法

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

圖4 改進的迭代運算流程圖

2 結果與分析

2.1 工具驗證

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

為了使對比結果更為直觀，采用相對偏差進行觀察，相對偏差定義如下：

式中：POld——原計算方法所得結果；

PNew——改進的新方法所得結果。

對比驗證結果如表1所示。由表可得，除第1組數據之外，所得σ、質量和數量校正因子的相對誤差均小于0.01%，可驗證該工具有良好的準確性。卡方誤差所代表的是估算值與輸入真實測量值之間的相對誤差平方。卡方誤差越小，所估算的值越接近真實測量值。在原計算方法中當δ＜ε（ε=1×10-9）時，才能說明該計算結果有效。

表1 對比驗證及誤差1）

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

2.2 結果對比

使用1.3節的12組CFM56-7B26/3型發動機數據，分別使用原計算方法和新計算工具進行運算，所得結果對比如下。

2.2.1 幾何平均粒徑

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

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

圖5 幾何平均粒徑對比

圖6 平均粒徑與質數比關系圖

2.2.2 nvPM數量濃度系統損失校正因子

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

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

圖8 nvPM數量濃度系統損失校正因子對比

2.2.3 nvPM質量濃度系統損失校正因子

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

圖9 nvPM質量濃度系統損失校正因子對比

3 再優化與討論

3.1 數據分析與工具優化

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

圖10 掌握數據與Lobo[9]數據對比

經數據分析后得到質數比隨發動機推力增長曲線和平均粒徑與推力線性關系圖如圖11所示。

圖11 推力與質數比及平均粒徑關系圖

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

推力與質數比關系式：

推力與平均粒徑關系式：

以推力為變量，代入式（2）、式（3）可得所對應的質數比和平均粒徑，將所得結果代入改進的計算流程程序可得到一個只有推力作為單一變量的計算工具。由于只導入了CFM56-7B26/3型發動機的相關關系式，所以該計算工具只對應于該型發動機的推力輸入和校正因子結果輸出。

3.2 結果分析

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

圖12 推力與數量校正因子關系圖

圖13 推力與質量校正因子關系圖

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

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

3.3 適用性討論

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

4 結束語

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

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

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

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