反應(yīng)傳遞多尺度耦合的擬顆粒模擬

葛蔚，李成祥，陳飛國

（1中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2中國科學(xué)院大學(xué)化學(xué)工程學(xué)院，北京 100049；3中國科學(xué)院綠色過程制造創(chuàng)新研究院，北京 100190；4中國科學(xué)院潔凈能源創(chuàng)新研究院，遼寧 大連 116023）

引 言

在化工過程中，化學(xué)反應(yīng)總是與反應(yīng)物、生成物或其他介質(zhì)的擴(kuò)散與流動(dòng)共存，而總體的反應(yīng)效果也受到這些過程的強(qiáng)烈影響，因此，反應(yīng)與傳遞的耦合是化學(xué)工程中反應(yīng)工程和傳遞過程共同關(guān)注的重要課題。但傳統(tǒng)理論對這種耦合的處理是松散而靜態(tài)的，即：微觀上根據(jù)活化能和反應(yīng)碰撞理論等可確定反應(yīng)的速率，而基于分子運(yùn)動(dòng)論和統(tǒng)計(jì)力學(xué)可確定擴(kuò)散和熱導(dǎo)率以及黏度等傳遞過程的物性參數(shù)，并且這些參數(shù)均為物質(zhì)宏觀狀態(tài)的函數(shù)，即本構(gòu)關(guān)系；宏觀上表達(dá)質(zhì)量、動(dòng)量和能量守恒的狀態(tài)演化方程組，如Navier-Stokes方程基于本構(gòu)關(guān)系可以封閉求解。這種宏微觀描述相對分離、通過物性相聯(lián)系的框架對一些簡單體系，如層流下反應(yīng)控制的均相反應(yīng)體系可能是良好的近似，但與工業(yè)上大多數(shù)多相反應(yīng)過程的真實(shí)圖景仍有很大的差別。

實(shí)際上，工業(yè)過程中反應(yīng)與傳遞的耦合大多緊密而復(fù)雜，難以截然區(qū)分為宏觀和微觀兩個(gè)尺度的描述[1]。比如，在廣泛應(yīng)用的多級(jí)孔催化劑中，微介孔中的擴(kuò)散由表面擴(kuò)散和努森擴(kuò)散主導(dǎo)，其速率強(qiáng)烈依賴于具體的孔徑和表面形貌等，而不是作為宏觀狀態(tài)函數(shù)的物性。同時(shí)實(shí)際的反應(yīng)網(wǎng)絡(luò)與速率也將依賴于這些結(jié)構(gòu)因素和受此影響的物質(zhì)分布。

為了應(yīng)對這種緊密耦合的復(fù)雜性，通常沿用傳統(tǒng)的描述框架，而以經(jīng)驗(yàn)的表觀動(dòng)力學(xué)和等效傳遞系數(shù)等獲得對特定系統(tǒng)的相對合理的描述。但這類描述因缺乏機(jī)理性，其通用性和預(yù)測性受到很大的限制。近年來，隨著對反應(yīng)過程精準(zhǔn)調(diào)控要求的不斷提高，以及過程強(qiáng)化[2-5]與微化工[6-7]等領(lǐng)域的迅速發(fā)展，對宏觀過程的微觀機(jī)理以及對微介觀過程本身（而不僅是它們對宏觀過程的影響）的認(rèn)識(shí)變得越來越重要。因此，前述傳統(tǒng)描述框架已不能滿足學(xué)科發(fā)展和工程應(yīng)用的需求。但另一方面，完全從微觀的原子、分子和基元反應(yīng)層面出發(fā)描述直至微秒和微米以上的時(shí)空多尺度行為也是不現(xiàn)實(shí)的。盡管嚴(yán)格的分子模擬與超級(jí)計(jì)算相結(jié)合已能復(fù)現(xiàn)億級(jí)原子體系的納秒級(jí)演化過程，但與上述時(shí)空尺度仍有約5個(gè)量級(jí)的差距[8]。而目前通過實(shí)驗(yàn)手段完整展示微介觀尺度的反應(yīng)傳遞演化過程面臨的挑戰(zhàn)更加巨大。因此亟需發(fā)展更加高效又具有足夠精度的微觀模擬和分析方法。

擬顆粒模擬（pseudo-particle modeling,PPM）[9-10]正是為此目的而提出的。它結(jié)合了分子模擬中簡化的硬球模型算法上的高效與精確性和軟球模型物理上的靈活性與算法上的并行性，利用當(dāng)前的超級(jí)計(jì)算系統(tǒng)即可實(shí)現(xiàn)微米-毫米與微秒-毫秒量級(jí)的反應(yīng)-傳遞過程直接微觀模擬，為描述連接宏微觀的多孔介質(zhì)和表-界面介尺度上復(fù)雜的反應(yīng)傳遞緊密耦合過程提供了一種有效手段。

1 擬顆粒模擬的基本思想與發(fā)展

分子動(dòng)力學(xué)模擬是在經(jīng)典物理框架下描述原子和分子微觀運(yùn)動(dòng)的最基本方法，其中大多數(shù)采用連續(xù)光滑的勢函數(shù)和時(shí)間驅(qū)動(dòng)的同步更新算法，或簡稱軟球（soft sphere,SS）方法。其優(yōu)勢是能靈活接受多種粒子間作用，適用面廣，易于并行編程，可擴(kuò)展性好，但也存在顯著的積分?jǐn)?shù)值誤差（截?cái)嗾`差）。而采用階躍勢函數(shù)和事件驅(qū)動(dòng)的異步更新算法的硬球（hard sphere,HS）方法正好相反，其粒子間碰撞可解析求解，沒有截?cái)嗾`差，并能準(zhǔn)確描述氣體運(yùn)動(dòng)，且計(jì)算效率非常高，但并行編程困難，可擴(kuò)展性差。如圖1所示，PPM則采用階躍勢函數(shù)，但允許粒子間微量的重疊，并采用時(shí)間驅(qū)動(dòng)的同步更新算法，從而結(jié)合了軟球和硬球方法各自的優(yōu)勢又避免了它們的部分缺點(diǎn)(圖1中a、b代表t+1時(shí)刻分兩步：a、b兩個(gè)事件)[10]。之后在建立擬顆粒與硬球物性間映射關(guān)系的基礎(chǔ)上，通過硬球與擬顆粒方法的空間分區(qū)耦合，即HS-PPM[11-12]，又進(jìn)一步提高了模擬的速度和可擴(kuò)展性（圖2[11，13]）。目前對大多數(shù)狀態(tài)下的氣體擴(kuò)散與流動(dòng)模擬，PPM可以比傳統(tǒng)分子模擬快3～4個(gè)量級(jí)，而模擬精度基本相當(dāng)。另外根據(jù)本征反應(yīng)動(dòng)力學(xué)確定粒子生成、刪除與重組的規(guī)則后，該方法還可以在微觀上自然耦合基元反應(yīng)和擴(kuò)散與微流動(dòng)過程，從機(jī)理上復(fù)現(xiàn)這種耦合帶來的復(fù)雜反應(yīng)網(wǎng)絡(luò)和表觀反應(yīng)動(dòng)力學(xué)（圖3）。

圖1 擬顆粒碰撞過程示意圖[11]Fig.1 Schematic diagram of the PPM collision process[11]

圖2 HS-PPM空間分區(qū)耦合[11,13]Fig.2 Coupling of spatial partitions of HS-PPM[11,13]

圖3 反應(yīng)-擴(kuò)散-流動(dòng)耦合模型示意圖Fig.3 Schematic diagram of reaction-diffusion-flow coupling model

采用PPM已成功復(fù)現(xiàn)了多個(gè)經(jīng)典的低速[10]和高速[14]流動(dòng)中、以及單相[15]與多相[16]體系的傳遞現(xiàn)象，并在沒有引入經(jīng)驗(yàn)參數(shù)的條件下與實(shí)驗(yàn)或理論結(jié)果吻合良好，證明了該方法的有效性。同時(shí)還發(fā)現(xiàn)了反應(yīng)與擴(kuò)散耦合可引起反應(yīng)位點(diǎn)附近非隨機(jī)的濃度漲落[17]，初步展示了該方法揭示反應(yīng)與傳遞耦合機(jī)理的能力。此外，采用HS-PPM還研究了催化劑孔道中積碳對氣體擴(kuò)散的影響，解釋了實(shí)驗(yàn)中觀察到的不同積碳模式導(dǎo)致的催化劑快速和緩慢失活的不同表現(xiàn)[13]。下面的三個(gè)實(shí)例將進(jìn)一步說明該方法的應(yīng)用前景與發(fā)展?jié)摿Α?/p>

2 多孔介質(zhì)中吸附過程的擬顆粒模擬與優(yōu)化

多孔介質(zhì)吸附氣體廣泛應(yīng)用于多種化工過程，如多相催化[18]、氣體凈化和分離[19]、環(huán)境保護(hù)[20]、CO2捕集[21]與新型儲(chǔ)氫[22]等。由于不同尺度孔道內(nèi)傳輸機(jī)制的差異，其吸附過程的隨機(jī)性和不均勻性顯著，尤其是孔道內(nèi)的擴(kuò)散阻力導(dǎo)致氣體向顆粒內(nèi)部擴(kuò)散的深度較小，制約了吸附效率的提高，也對其優(yōu)化設(shè)計(jì)提出了挑戰(zhàn)。本節(jié)為采用HS-PPM研究低壓下氮?dú)庠趶?fù)雜孔道內(nèi)的擴(kuò)散-吸附耦合過程，并提出對孔道結(jié)構(gòu)優(yōu)化的建議。

模擬基于氣體碰撞吸附理論，假設(shè)氣體碰到吸附位點(diǎn)后，如其動(dòng)能超過吸附熱，則自發(fā)物理吸附。而已經(jīng)吸附的氣體分子如果與自由運(yùn)動(dòng)的氣體分子碰撞且獲得的能量超過吸附熱，則發(fā)生脫附。

如圖4所示，所考慮的多孔吸附劑以活性炭為模型，主要由石墨微晶結(jié)構(gòu)組成，氣體沿X方向擴(kuò)散進(jìn)入吸附劑。模擬區(qū)域沿X方向分為氣體分子數(shù)密度（0.1 /nm3）以及溫度（77K）恒定的控制區(qū)和靠近顆粒表面的復(fù)雜孔道吸附區(qū)。模擬中在Y和Z方向施加周期性邊界條件，而氣體在X方向離開孔道結(jié)構(gòu)即被刪除。構(gòu)建活性炭孔道結(jié)構(gòu)時(shí)，首先生成石墨的亂層結(jié)構(gòu)，然后根據(jù)孔隙率ε以及孔徑分布等參數(shù)生成隨機(jī)孔道結(jié)構(gòu)[11]。前期模擬結(jié)果表明，當(dāng)氣體從X方向擴(kuò)散進(jìn)入孔道時(shí)，100nm深度處的吸附速率已經(jīng)遠(yuǎn)低于外表面的初始吸附速率。因此，就微觀模擬目前能有效考慮的快速吸附過程而言，X、Y和Z方向的模擬區(qū)域邊長大約為100、20和20nm基本可以滿足模擬精度要求。同時(shí)為了在Y和Z方向施加周期性邊界條件，邊長分別圓整為石墨相應(yīng)方向晶格長度的400、100和60倍，即96.72 、20.94 和20.1 nm。此外X方向還設(shè)厚度為20nm的控制區(qū)以有效抑制進(jìn)入孔道的氣體溫度和數(shù)密度等的漲落。吸附劑平均ε=0.396 ，孔徑在0.4 ～2.0 nm之間隨機(jī)分布，平均為1.0 nm。吸附位點(diǎn)隨機(jī)分布于孔道內(nèi)表面，數(shù)密度為3個(gè)/nm3，氮?dú)夥肿拥挠睬虍?dāng)量直徑為0.3777 nm，其在石墨上的吸附熱為20.71 kJ/mol。擬顆粒模擬的時(shí)間步長為72fs，模擬總時(shí)間為215.14 ns。

圖4 氣固吸附過程模型示意圖Fig.4 Schematic diagram of the gas-solid adsorption model

模擬得到的孔道內(nèi)吸附氣體密度（Γ）沿進(jìn)入孔道深度（δ）的分布如圖5(a)所示，主要吸附區(qū)靠近外表面，并在深度約27.1 nm處吸附量趨近于0。究其原因，理論上單分子層吸附后的孔道內(nèi)能夠繼續(xù)擴(kuò)散的最小孔道直徑為3個(gè)氣體分子直徑，約1.1 nm，而本例中平均孔徑只有1.0 nm，導(dǎo)致孔道堵塞、吸附劑利用率低。該現(xiàn)象在快速吸附中更加明顯。相應(yīng)地，如圖5(b)所示，當(dāng)生成的平均孔徑為2.0 nm時(shí)，氣體幾乎可以吸附在整個(gè)孔道區(qū)域?？椎榔骄睆捷^大時(shí)，孔徑分布也較寬，而受孔徑與內(nèi)表面積復(fù)雜影響的吸附量波動(dòng)也相應(yīng)地較大，因此，圖5(a)、(b)僅作定性對比，如需獲得吸附氣體密度隨深度更光滑的分布，需統(tǒng)計(jì)大量相同條件下模擬樣本的平均值。最近相應(yīng)的實(shí)驗(yàn)研究也證實(shí)了模擬的預(yù)測并已應(yīng)用于一種快速吸附裝置的改進(jìn)[23]。

圖5 孔道內(nèi)吸附氣體密度沿X方向的分布Fig.5 Distribution of adsorbed gas density in the pores along the X direction

3 多級(jí)孔催化的擬顆粒模擬與優(yōu)化

圖6 多級(jí)孔道結(jié)構(gòu)模型[24]Fig.6 Hierarchical pore structure model[24]

為優(yōu)化孔道直徑及其體積占比，進(jìn)行了一系列模擬對比。首先，保持孔徑D1=1.23 和D3=24.63 以及孔體積占比RD1=29%、RD2=50%、RD3=21%不變，調(diào)控D2在4～15之間變化。如圖7模擬結(jié)果顯示，C2和C3的臨界值的時(shí)空收率（yield,y）和反應(yīng)物轉(zhuǎn)化率（conversion,CVR）隨著ε顯著增加，然而對于同一ε，y和CVR幾乎不受D2影響，說明只改變中等大小的孔徑對催化劑整體性能的影響不大。為了進(jìn)一步獲得不同孔徑對催化性能的影響，在保持ε=0.36 不變且D1=1.23 ，D2=12.32 時(shí)調(diào)控最大孔徑D3（受到計(jì)算規(guī)模的影響，孔體積占比也相應(yīng)改變）。結(jié)果如圖8所示，當(dāng)D3大于12.32 時(shí)，Y和CVR隨D3的減小而增大，表明孔隙率不變時(shí)，適當(dāng)減小最大孔徑及其體積占比有利于提高催化性能。當(dāng)進(jìn)一步減小D3時(shí)，將出現(xiàn)D3＜D2，為了保持其他量基本不變而實(shí)現(xiàn)對最大孔徑的獨(dú)立調(diào)控，保持ε=0.30 ，D1=1.23 不變，而調(diào)控D2=D3同時(shí)變化。由上面結(jié)果可知D2的影響可以忽略。結(jié)果如圖8所示，y和CVR在D3約為7.39 時(shí)獲得最大值，并在D3＜7.39 后明顯降低。

圖7 孔體積占比不變，不同孔徑（D2）時(shí)的時(shí)空收率(a)和轉(zhuǎn)化率(b)隨孔隙率的變化Fig.7 Variation of the yield(a)and conversion(b)versus the porosity with different pore diameter(D2)under constant proportion of pore volume

圖8 時(shí)空收率(a)和轉(zhuǎn)化率(b)隨D3的變化Fig.8 Variation of the yield(a)and conversion(b)versus D3

以上的討論表明，在孔隙率和孔徑分布給定的多級(jí)孔催化劑中，中孔徑對催化劑整體性能的影響不大，而大孔的孔徑對總時(shí)空收率存在最優(yōu)值。目前正與相關(guān)企業(yè)合作改進(jìn)相應(yīng)催化劑制備工藝，以期獲得實(shí)際的催化效果提升。

4 復(fù)雜多相流動(dòng)過程的擬顆粒模擬與反應(yīng)過程模擬前景

作為簡化的分子動(dòng)力學(xué)方法，PPM可自然處理多相體系相界面處的微觀作用。Chen等[27]耦合PPM和軟球分子動(dòng)力學(xué)，分別表達(dá)氣體分子大部分時(shí)間自由飛行而偶有瞬時(shí)碰撞和液體分子間存在持續(xù)的相互作用等不同特征的行為，建立了氣液兩相流動(dòng)的微觀模擬方法。對納微尺度下兩平板間氣液運(yùn)動(dòng)的模擬表明，在不同體積力作用下可依次出現(xiàn)泡狀、彈狀、環(huán)狀和霧（散）狀流型（圖9，圖中以氬原子直徑σ=3.4 ×10-10m為長度單位），該趨勢與經(jīng)典的納微流動(dòng)實(shí)驗(yàn)[28]和分子模擬[29]結(jié)果符合良好，而在模擬速度與規(guī)模上優(yōu)勢明顯，展示了PPM在復(fù)雜多相流動(dòng)研究中的巨大潛力。

圖9 不同體積力下平板間氣液兩相流型[27]Fig.9 Steady state flow patterns under different intensities of gravity between parallel plates[27]

作為氣體運(yùn)動(dòng)的離散模擬方法，PPM還可與光滑粒子流體動(dòng)力學(xué)方法（smoothed particle hydrodynamics，SPH）等其他流體模擬的無網(wǎng)格方法耦合，形成多相流動(dòng)的高精度多尺度離散模擬方法。與基于網(wǎng)格的常規(guī)計(jì)算流體動(dòng)力學(xué)相比，純粹拉格朗日方法的多相流動(dòng)全離散模擬具有以下主要優(yōu)點(diǎn)：①易處理大變形及快速移動(dòng)的相界面；②可自發(fā)捕捉和追蹤間斷面；③易耦合多相傳質(zhì)傳熱和反應(yīng)等過程；④易實(shí)現(xiàn)可擴(kuò)展性強(qiáng)的并行計(jì)算。在PPM-SPH耦合模擬中，擬顆粒與SPH粒子碰撞時(shí)即時(shí)改變速度，而碰撞過程產(chǎn)生的沖量轉(zhuǎn)化為碰撞時(shí)間內(nèi)的持續(xù)作用力施加于SPH粒子。PPM-SPH模擬的氣液兩相霧化過程如圖10所示。在8mm×2mm×2mm的模擬區(qū)域中，液體從頂部以一定速度注入，受橫向氣流作用液柱偏離入射方向，并發(fā)生周期性剝落破碎成特定尺寸分布的液滴群。橫流霧化模式?jīng)Q定于液氣動(dòng)量比q和氣體韋伯?dāng)?shù)We[30]，上述模擬條件為q=19.5 和We=30，對應(yīng)于袋狀破碎霧化，形態(tài)合理。與采用流體體積法（volume of fluid，VOF）或水平集法（level set，LS）等界面捕捉技術(shù)的連續(xù)方法模擬相比，PPM-SPH模擬的霧化結(jié)構(gòu)更為精細(xì)（可達(dá)微米級(jí)），計(jì)算過程也更為高效。

圖10 橫流霧化過程模擬Fig.10 Atomization of a liquid jet in a crossflow

對上述不同尺度多相流動(dòng)的全離散模型還可方便地引入傳熱傳質(zhì)和反應(yīng)模型。擬顆粒的均方位移可自然表達(dá)多相體系中不同組分的傳質(zhì)過程；而表示分子熱運(yùn)動(dòng)的擬顆粒脈動(dòng)速度反映了體系的溫度，并通過擬顆粒的脈動(dòng)、相互作用和整體運(yùn)動(dòng)分別表達(dá)熱傳導(dǎo)與對流傳熱。另外也可以通過求解能量方程獲得SPH粒子溫度分布情況，從而形成多相復(fù)雜體系的傳遞和反應(yīng)模型。

5 總結(jié)與展望

擬顆粒模擬提供了反應(yīng)傳遞耦合過程微觀模擬的一種高效方法，并具有較高的準(zhǔn)確性和通用性。采用該方法能系統(tǒng)獲得各種復(fù)雜條件下的反應(yīng)與傳遞的速率關(guān)系，而不是反過來通過經(jīng)驗(yàn)地建立這些關(guān)系描述耦合的效果，因此它是一種機(jī)理性和預(yù)測性的模擬手段。通過超級(jí)計(jì)算的應(yīng)用，能有效覆蓋從納米到毫米和從皮秒到毫秒的反應(yīng)傳遞過程的介尺度區(qū)間，也有利于探索其中的介科學(xué)原理。

已有工作和本文的實(shí)例還表明，通過PPM可以深入全面地優(yōu)化反應(yīng)傳遞耦合過程，包括反應(yīng)條件、物性和反應(yīng)器與介質(zhì)結(jié)構(gòu)等，當(dāng)然由于本征反應(yīng)動(dòng)力學(xué)和原子、分子模型的準(zhǔn)確性與非確定的問題，其優(yōu)化結(jié)論往往只有定性意義，但配合相應(yīng)的實(shí)驗(yàn)研究可以更快地獲得具有工程意義的結(jié)論。同時(shí)隨著對反應(yīng)和微觀分子間作用的描述水平的不斷提高，特別是人工智能與量子化學(xué)計(jì)算方法的改進(jìn)與耦合應(yīng)用，PPM的準(zhǔn)確性也會(huì)相應(yīng)提高，并實(shí)現(xiàn)工程問題的定量設(shè)計(jì)與優(yōu)化。

目前由于硬球方法的限制，PPM對非常稠密的氣體、液體或流體相中有復(fù)雜反應(yīng)的體系還難以準(zhǔn)確描述。今后將發(fā)展多體和可變徑硬球的擬顆粒模型，并結(jié)合更靈活多樣的反應(yīng)規(guī)則，實(shí)現(xiàn)更通用準(zhǔn)確的PPM。

符號(hào)說明

Ec——由反應(yīng)/吸附活化能決定的臨界能量，J

Et——顆粒碰撞時(shí)相對速度對應(yīng)的動(dòng)能，J

e12——碰撞恢復(fù)系數(shù)

g——重力加速度，m/s2

m1,m2——顆粒質(zhì)量，kg

P1,P2——顆粒位置，m

t——時(shí)間，s

V1,V2——顆粒碰撞前的速度，m/s

vc——與臨界能量Ec對應(yīng)的臨界速度，m/s

vt——顆粒碰撞時(shí)的相對速度，m/s

ε——多孔介質(zhì)孔隙率

下角標(biāo)

1,2——分別表示顆粒1和顆粒2