葛蔚，李成祥，陳飛國

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

引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 孔體積占比不變，不同孔徑（D2）時的時空收率(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 時空收率(a)和轉(zhuǎn)化率(b)隨D3的變化Fig.8 Variation of the yield(a)and conversion(b)versus D3

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

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

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

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

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

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

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

5 總結(jié)與展望

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

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

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

符號說明

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

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

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

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

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

P1,P2——顆粒位置，m

t——時間，s

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

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

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

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

下角標

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