基于LJ-GCEMC的炭化木纖維濾芯超微粒子吸附模擬

杜丹豐，馬巖，楊春梅，郭秀榮

（東北林業(yè)大學 交通學院，黑龍江哈爾濱150040）

近年來，國際上廣泛開展了柴油機超微顆粒排放方面的研究[1-2]。根據顆粒粒徑的不同，柴油機排放顆粒物（Particulate Material，PM）可分為核模態(tài)（粒徑小于50 nm）、積聚模態(tài)（粒徑小于100 nm）和粗糙模態(tài)（粒徑大于100 nm），其中核模態(tài)和積聚模態(tài)約占排氣總顆粒數的90%以上[3-5]，這些粒子（亦稱超微粒子）會對人體肺部和腦部造成重大傷害，還可能引發(fā)遺傳性基因突變[6]。近年的柴油發(fā)動機雖然在減少顆粒物質量濃度的設計方面取得了長足的進步，但是超微顆粒數濃度的上升卻成為困擾發(fā)動機繼續(xù)升級的障 礙[7]。微粒捕集器（Diesel Particulate Filter，DPF）被認為是解決柴油機微粒排放問題最有效的手段之一，已成為當前柴油車微粒排放控制技術的研究熱點[8]，但深受價格、再生、柴油品質以及過濾超微顆粒等諸多問題的困擾，在一些國家一直未得以良好的應用。

針對以上問題，利用特殊工藝對微米木纖維進行整體炭化和活化，制成綠色環(huán)保、成本低、凈化效率高、具有活性炭性質的炭化微米木纖維（Carbonized Micron Wood Fiber，CMWF）過濾體（濾芯）[9]。該過濾體的突出優(yōu)點是對超微顆粒吸附能力強。本文利用混合Lennard Jones 勢能——巨正則系綜蒙特卡羅（混合LJ-GCEMC）法模擬超微粒子在CMWF 孔隙中的吸附過程，以期對加工工藝起到一定的指導作用。

1 CMWF 濾芯制作工藝及特性

CMWF是順紋刨切木材制成的幾十微米厚、幾百微米寬、幾千微米長的高強度木纖維，經過膠合和熱壓重組成型，再經過特殊工藝炭化和活化而成的新型過濾材料[9]。文獻[10]以微觀力學和細胞學理論為基礎獲得了傳統(tǒng)木纖維端面細胞結構，文獻[11]指出，在木材微米切削時，細胞會被切開，可大幅度增加管胞容腔的接觸面積，提高PM的吸附速度及在細胞中的吸附量，圖1為傳統(tǒng)木纖維與微米木纖維端面細胞結構圖。

圖1 木纖維端面細胞結構圖

CMWF 過濾體制作過程如圖2所示。利用自行研制的數控大型微米刨片機加工微米級厚度刨花，并將刨花加工成微米木纖維，進而利用微纖絲模壓制品加工專用壓機制成微米木纖維濾芯，利用特殊炭化和活化工藝初步研制出CMWF 濾芯并組裝成可拆卸式DPF[9]。

圖2 CMWF 濾芯制作過程圖

2 柴油發(fā)動機排放顆粒物特性分析

文獻[12-13]指出，碳氫燃料的不完全燃燒所生成的碳煙是以碳原子作為主要成分并含有10%～30%氫原子的碳氫化合物所組成。在缺氧條件下，碳氫燃料中的重質烴在液態(tài)下直接脫氫碳化，成為焦碳狀的液相析出型態(tài)碳粒，粒度一般比較大；而輕質烴首先氣化，然后裂解成甲烷和乙烯之類的低分子碳氫化合物以及多環(huán)芳烴，它們不斷脫氫形成原子級的碳粒子并聚合成直徑2 nm 左右的碳煙核心，碳煙核心再碰撞凝聚并與其它有機物附著，碳核增大為直徑20～30 nm 左右的碳煙基元，繼而形成積聚態(tài)的超微粒子。

文獻[14]指出，柴油機排放顆粒物中約含60%～80%的元素碳和20%～30%的有機組分。文獻[15]指出，柴油機50%負荷下50 nm和76.7 nm 均相粒徑段粒子非揮發(fā)性組分質量分數大于90%，且排氣顆粒主要由10 nm 左右的初始粒子形成的聚團構成。文獻[16]發(fā)現稀釋參數（稀釋比例、稀釋空氣溫濕度以及稀釋混合氣在通道內的停留時間等）的變化顯著影響核模態(tài)粒子的數密度和粒徑分布，而對積聚模態(tài)粒子影響較小。文獻[17]通過試驗獲得柴油機在1400 r·min-1，50%負荷工況下超微粒子的顆粒數密度。文獻[18-19]指出，積聚態(tài)顆粒物的粒數濃度較穩(wěn)定，在實驗中表現出很好的可重復性。這些積聚態(tài)粒子絕大部分會被直接排入大氣，成為大氣中元素碳（Element Carbon，EC）的主要來源。

通過以上分析可知，在柴油車排放物中，從數量上看，大部分為由碳煙基元凝并而成的超微粒子，所以研究CMWF 濾芯對超微粒子的吸附容量意義重大。本文主要利用混合LJ-GCEMC 法對CMWF 濾芯吸附超微粒子的顆粒數量加以統(tǒng)計，在分析過程中，將含超微粒子的尾氣視為流體進行研究。根據上述文獻，本文假設積聚態(tài)粒子完全由元素碳組成，其平均直徑為76.7 nm，到達CMWF 濾芯處的粒子數密度為9×106個/cm3[20]，其中元素碳部分主要是由原子級的碳粒子碰撞凝聚而成，可將其看做是碳原子簇。

3 CMWF 濾芯吸附過程模擬

吸附實驗一般只能對吸附現象進行宏觀上的測量和分析，無法對吸附的微觀過程進行描述，而利用計算機模擬具有活性炭性質的CMWF 濾芯的吸附，可為了解吸附的微觀過程和活性炭吸附機理提供信息。混合LJ-GCEMC 法是模擬吸附過程的一種有效方法。目前，國內外研究者已利用GCEMC方法對碳孔墻吸附氣體分子及金屬納米粒子特性進行了研究[21-23]，本文在此基礎上對CMWF 濾芯吸附超微碳煙顆粒的相關特性進行初步探討。

3.1 勢能模型

在分子動力學模擬中，粒子之間相互作用力計算普遍采用的勢能模型是Lennard-Jones（LJ）勢能模型[22]。ψLJ為Lennard-Jones（LJ）勢能，可表示為

式中，r為粒子間的距離；σ為尺寸參數；ε為能量參數。為了利用分子動力學模擬粒子之間的相互作用，本文采用在巨正則系綜中廣泛應用的截斷漂移LJ勢能（Cut -Shifted Lennord - Jones（CSLJ）potential）來表征流體超微粒子之間的相互作用。

式中，ψff為流體粒子之間的勢能，J；rc為截斷半徑，rc＝3.8σff，σff為超微粒子尺寸參數。

流體相與單個狹縫炭孔墻之間的勢能采用平均場理論中的10-4-3勢能模型[24]。圖3為流體在CMWF濾芯狹縫孔中的示意圖，圖中狹縫孔孔徑H是炭孔墻的中心原子到正對面炭孔墻中心原子的距離。

圖3 狹縫孔中超微粒子示意圖

定義與狹縫炭孔墻垂直的方向為z軸，則勢能可表述為

式中，ρw為炭孔墻的數密度，ρw＝114 nm-3；Δ表示石墨的晶面間距，Δ＝0.335 nm；z表示流體粒子和炭孔墻之間的距離，nm；εfw和σfw分別為超微粒子和炭孔墻的交互作用參數，它們可以從Lorentzberthelot 混合規(guī)則獲得[25]：

式中：εff為超微粒子的能量參數；σww為炭孔墻的尺寸參數；εww為炭孔墻的能量參數，σww＝0.34 nm，εww/kB＝28 K[26]。

對固定的孔徑H，流體相在活性炭孔中的總勢能ψT可以表示為兩種相互作用勢能之和，即流體相粒子之間的勢能ψff和流體相粒子與炭孔墻之間的勢能（ψfw（z）+ψfw（H-z））之和：

3.2 超微粒子勢能參數的確定

為了便于計算，只考慮超微粒子中最近碳原子之間的相互作用[27]。假設碳原子簇由大量配位數為6的簡立方結構組成，即每個中心碳原子周圍有6個最近的碳原子，2個碳原子可形成一個鍵，這樣，每個碳原子會與最近的碳原子形成6個鍵；然而，由于每個鍵屬于2個碳原子，所以，每個碳原子的有效鍵數應為總鍵數的一半，對于簡立方結構來說，一個碳原子的有效鍵數是3，一個內部碳原子的結合能等于這3個鍵的總能量，這樣，超微粒子的結合能即為粒子中所有鍵能之和。

根據上面的討論可獲得超微粒子的結合能。為了簡化問題，假設超微粒子是一個理想的立方形，其內部在三維空間中具有相同數量的炭原子單元。如果這個立方形超微粒子的邊界尺寸是D，那么，納米粒子的體積和表面積分別是D3和6D2。由于簡立方、體心立方和面心立方的體積變化相差很小，所以本文忽略不同晶格結構之間的差異，假設所有的晶格都是簡立方結構，那么，立方超微粒子的總碳原子數n為D3/d3，表面碳核數N為6D2/d2，這里d表示碳原子的尺寸。這樣，便有N為6n2/3，內部炭原子數為（n-6n2/3）。

令α為一個內部碳原子的鍵數，考慮到表面弛豫，在晶格中每個表面炭原子的鍵數超過一半是懸掛鍵，本文取表面炭原子鍵數為內部炭原子鍵數的1/4，則超微粒子的結合能可表示為

式中，Ebond是鍵能。將式（6）進行簡化，得：

式中，E0為超微粒子的結合能

根據E0及文獻[28]可以算出超微粒子的能量參數 εff/kB＝1620K，其中kB為Boltzmann 常數，J·K-1。另外，超微粒子的尺寸參數可近似為其邊界尺寸，本文取σff為76.7 nm。

3.3 巨正則系綜模擬

3.3.1 巨正則系綜模擬理論

巨正則系綜的配分函數[29]可表示為

式中，N是系綜體系的分子數；ri是第i個分子的向量坐標；λ為de Broglie 波長是化學勢，μ＝kBT lnλρ；h為Planck 常數。因此，熱力學量X的系綜平均可表示為

式中＜＞表示系綜平均。方程（8）～（10）是進行模擬的基本方程。為了方便和簡化計算，文中模擬所涉及量都以流體相粒子的尺寸參數和能量參數作為對比量進行無量綱化，具體形式如下：

式中，上標有星號的量均為無量綱量；μ為化學勢，J；V為體積，nm3；T為溫度，K；p為壓強，Pa；ρ是數密度；H是狹縫孔徑。

在μVT系綜熱力學平衡時，孔中流體相與主體流體相μ 相同。由熱力學基本關系式可知，μ與主體流體的壓力p 直接相關。由于p是操作變量，比較直觀，因此可采用文獻[30]提出的MBWR 狀態(tài)方程來表述主體流體的壓力。

3.3.2 巨正則系綜模擬過程

模擬初始，給定系綜的壓強p和溫度T，取狹縫孔中l(wèi)×l×H的模擬盒子，對其吸附流體相的過程進行GCEMC 模擬，l為模擬盒子的長如圖3所示，與炭孔墻垂直的方向為z軸，當狹縫孔的孔徑H 一定時，在x方向和y方向采用周期性邊界條件。

在每一步GCEMC 模擬中，模擬盒子里的流體超微粒子存在3種等幾率的行為[27,30]，分別是移動、產生和刪除。其接受概率：

1）移動概率ζm，模擬盒子里的超微粒子從一個位置移動到另一個位置：

2）產生概率ζc，在模擬盒子里產生一個新粒子：

3）刪除概率ζd，從模擬盒子里刪除一個粒子：

式中，N為流體粒子數，△ψ為經過移動、產生和刪除某一種行為后系統(tǒng)總勢能的變化量。

模擬時，化學勢、體積和溫度固定，以相同的概率進行插入、刪除和移動3種操作。在初始構型中，假設在模擬盒子內有60個超微粒子。每次模擬都產生2×106個構象，最初的1×106個構象用來使體系達到平衡而被丟棄，只對后面的1×106個構象進行統(tǒng)計平均，計算所需的物理量。本文計算采用Accelrys公司Cerius2分子模擬軟件包，在Origin300服務器上進行。

4 結果與討論

在加工制造CMWF 濾芯的過程中，可以通過控制微米木纖維的填充率及炭化活化主要影響參數確定其比表面積及孔容孔徑分布。在CMWF 濾芯中，有大量的微孔、中孔和大孔存在，這些孔可有效過濾和吸附不同粒徑的PM。本研究的主要目的是通過混合LJ-GCEMC 法模擬獲得CMWF 對積聚態(tài)粒子的吸附特性，從而對確定CMWF 濾芯的工藝參數起到一定的指導作用。

模擬中，孔中粒子平均數密度ρ*T可表示為[22]

式中，ρ*（z*）為局域密度，可通過式（10）獲得。

式中，S*為狹縫墻面積；〈N（z*）〉為孔中微元胞腔（S*·Vz*）內粒子數的系綜平均。

根據柴油車尾氣流經DPF時的實際壓力與溫度，確定GCEMC 模擬中壓力和溫度變化范圍分別為102～120 kPa和300～500 K。

4.1 CMWF 濾芯吸附數密度隨孔徑的變化

為了研究孔徑對吸附數密度的影響，本文對溫度為350 K時不同壓力條件下（105 kPa、110 kPa、115kPa）CMWF 濾芯的吸附數密度進行了模擬。從圖4可以看出，在孔徑為380 nm和壓力為115 kPa時吸附數密度最大，可達到3.7×106個/cm3。同時可看到當孔徑小于200 nm時，吸附數密度隨著孔徑增大而增加；當孔徑從200 nm向320 nm 過渡時，吸附數密度急劇增加；當孔徑從320 nm向420 nm過渡時，吸附數密度相對較大；當孔徑從420 nm向650 nm 過渡時，吸附數密度隨孔徑增大急劇減少；當孔徑大于650 nm時，吸附數密度隨孔徑的增大變化不明顯。以上分析說明孔徑過小或過大都會阻礙吸附數密度的增加，在本研究設定的條件下，孔徑范圍應為320～420 nm。另外，從圖4還可以看出，隨著壓力的升高吸附數密度會有所增加。

圖4 不同壓力條件下粒子吸附數密度隨孔徑的變化曲線

4.2 CMWF 濾芯吸附數密度隨壓力的變化

本文計算了孔徑為400 nm時不同溫度條件下（350 K、400 K和450 K）CMWF 濾芯吸附數密度隨壓力的變化關系（圖5）。在相同溫度下，吸附數密度隨壓力升高而增加。壓力在110 kPa 前，吸附數密度隨著壓力的升高而迅速增加；之后隨壓力增加，吸附數密度的上升趨勢有所減緩，當大于115kPa時，吸附數密度基本保持不變。另外，從圖5還可看出，隨溫度升高吸附數密度會有所減少。

圖5 不同溫度條件下粒子吸附數密度隨壓力變化曲線

5 結論

利用綠色環(huán)保、成本低、凈化效率高、具有活性炭性質的CMWF 濾芯可有效凈化柴油車尾氣中的超微碳煙顆粒。超微粒子組成復雜，目前利用GCEMC 進行吸附過程模擬的相關研究較少，本文只是一個探索性的理論研究，具體結論如下：

1）柴油機排放顆粒物中，超微粒子占排氣總顆粒數量的90%以上，且絕大多數由碳煙基元凝聚而成，可利用混合LJ-GCEMC 法模擬粒子在CMWF 濾芯中的吸附過程。

2）當溫度、壓力和孔徑分別為350 K、115 kPa和380 nm時，CMWF 濾芯的吸附數密度最大，可達到3.7×106個/cm3。隨著孔徑的變化，吸附數密度呈現先增大，然后減小，最后趨于平緩的趨勢。模擬結果表明，CMWF 濾芯孔徑過小或過大都會阻礙超微粒子吸附數密度的增加，在本研究設定的條件下，孔徑范圍應為320～420 nm。

3）當溫度保持不變且孔徑為400 nm時，吸附數密度隨壓力升高而增加。其中，壓力小于110kPa時，吸附數密度升高迅速；當壓力在110kPa和115kPa之間時，吸附數密度上升趨勢有所減緩，當壓力大于115 kPa時，吸附數密度基本保持不變。另外，模擬結果還表明，當孔徑不變時，隨溫度升高吸附數密度會略有減少。

[1]Shi J P,Mark D,Harrison R M.Characterization of particles from a current technology heavy-duty diesel engine[J].Environ Sci Tech,2000,34（5）∶748-755.

[2]Mathis U，Ristim ki J，Mohr M，et al.Sampling conditions for the measurement of nucleation mode particles[]，2004，38（12）：1149-1160.

[3]M.Matti Maricq.Chemical characterization of particulate emissions from diesel engines：A review[J].Journal of Aerosol Science 2007，38（11）：1079-1118.

[4]李新令，黃震，王嘉松，等.二甲醚發(fā)動機超細顆粒排放屬性實驗研究[J].科學通報，2007（7）：1706-1713.

[5]Brown J E，Clayton M J，Harris D B，et al.Comparison of the particle size distribution of heavy-duty diesel exhaust using a dilution tailpipe sampler and an in-plume sampler during on-road operation[J].J ournal of the Air & Waste Man-age Asso ciation，2000，50（8）：1407-1416.

[6]SARNAT J A，MARMUR A，KLEIN M，et al.Fine particle sources and cardiorespiratory morbidity：an application of Chemical Mass Balance and factor analytical sourceapportionment methods[J].Environmental Health Perspectives，2008，116（4）：459-466.

[7]陳慶平，余明洲.柴油發(fā)動機尾氣顆粒物研究進展[J].浙江水利水電?？茖W校學報，2010（3）：41-43.

[8]龔金科，伏軍，龍罡，等.汽車尾氣微粒過濾器噴油助燃再生控制策略研究[J].環(huán)境工程學報，2011（7）：1603-1608.

[9]Guo Xiurong，Du Danfeng，Wang Fenghu，Ma Yan.Study on test instrument and filtration theory of the carbonized micron wood fiber DPF[J].Microporous and Mesoporous Materials，2011，142（2-3）：655-660.

[10]馬巖.木材橫斷面六棱規(guī)則細胞數學描述理論研究[J].生物數學學報，2002（1）：64-68.

[11]馬巖.利用微米木纖維定向重組技術形成超高強度纖維板的細胞裂解理論研究[J].林業(yè)科學，2003（5）：111-115.

[12]L.Morawska，Z.Ristovski，E.R.Jayaratne，D.U.Keogh，X.Ling.Ambient nano and ultrafine particles from motor vehicle emissions：Characteristics，ambient processing and implications on human exposure[J].Atmospheric Environment，2008，42（35）：8113-8138.

[13]張洪濤，莫建文，邱卓丹.柴油機碳煙的生成及控制[J].廣西工學院學報，2000（12）：73-76.

[14]Kleeman M J，Schauer J，Cass G R.Size and composition distribution of fine particulate matter emitted from motor vehicles[J].Environmental science and technology，2000，34（7）：1132-1142.

[15]H.Burtscher.Physical characterization of particulate emissions from diesel engines：a review[J].Aerosol Science，2005，36（7）：896-932.

[16]付娟，寧智，蘇麗萍，等.汽車排氣尾流中微粒分布及其變化特性的研究[J].環(huán)境科學，2007（5）：952-957.

[17]李新令，黃震.柴油機排氣尾流中核模態(tài)顆粒數濃度和粒徑分布變化特性[]中國科技論文在線，（8）：613-616.

[18]姚志良，劉榮.柴油車排放顆粒物研究進展[J].科技導報，2011，29（8）：75-79.

[19]Harley R A，Marr L C，Lehner J K，et al.Changes in motor vehicle emissions on diurnal to decadal time scales and effects on atmospheric composition[J].Environmental Science and Technology，2005，39（14）：5356-5362.

[20]劉志華，葛蘊珊，丁焰，等.柴油機和LNG 發(fā)動機排放顆粒物粒徑分布特性研究[J].內燃機學報，2009（6）：518-522.

[21]S.V.Nedea，A.J.H.Frijns，A.A.van Steenhoven，A.P.J.Jansen，A.J.Markvoort，P.A.J.Hilbers.Density distribution for a dense hard-sphere gas in micro/nano-channels：Analytical and simulation results[J].Journal of Computational Physics，2006，219（2）：532-552.

[22]Xiang Zhang，Xiaohong Shao，Wenchuan Wang，et al.Molecular modeling of selectivity of single -walled carbon nanotube and MCM -41 for separation of methane and carbon dioxide[J].Separation and Purification Technology，2010，74（3）：280-287.

[23]W.H.Qia，M.P.Wanga，W.Y.Hu.Calculation of the cohesive energy of metallic nanoparticles by the Lennard-Jones potential[J].Materials Letters，2004，58（11）：1745- 1749.

[24]Piotr Kow alczyk，Hideki Tanaka，Robert Holyst，et al.Storage of Hydrogen at 303K in Graphite Slitlike Pores From Grand Canonical Monte Carlo Simulation[J].J Phys.Chem.B，2005，109（36）：17174-17183.

[25]張翔，馬聰聰，張學軍，等.空氣中微量1，1，1，2-四氟乙烷吸附分離的分子模擬[J].物理化學學報，2010，26（6）：1637-1642.

[26]Xuan Peng，Jinsong Zhao，Dapeng Cao.Adsorption of carbon dioxide of 1-site and 3-site models in pillared clays：A Gibbs ensemble Monte Carlo simulation[J].Journal of Colloid and Interface Science，2007，310（2）：391-401.

[27]Kittel，C.Introduction to Solid State Physics[K].USA：John Wiley & Sons，Inc.，1976：173.

[28]王風云，檀革江，陳民生.勢能參數的確定與估算[J].華東工學院學報，1992（3）：35-42.

[29]Daan Frenkel.Understanding Molecular Simulation -From Algorithms to Applications[M].Amsterdam：Academic press，2002：127-135.

[30]J.K.Johnson，J.A.Zollweg，K.E.Gubbins.The Lennard -Jones Equation of State Revisited[J].Molecular Physics，1993，78（3）：591-618.