反應傳遞多尺度耦合的擬顆粒模擬

葛蔚，李成祥，陳飛國

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

引 言

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

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

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

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

1 擬顆粒模擬的基本思想與發展

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

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

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

圖3 反應-擴散-流動耦合模型示意圖Fig.3 Schematic diagram of reaction-diffusion-flow coupling model

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

2 多孔介質中吸附過程的擬顆粒模擬與優化

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

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

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

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

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

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

3 多級孔催化的擬顆粒模擬與優化

圖6 多級孔道結構模型[24]Fig.6 Hierarchical pore structure model[24]

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

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

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

4 復雜多相流動過程的擬顆粒模擬與反應過程模擬前景

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

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

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

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

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

5 總結與展望

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

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

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

符號說明

Ec——由反應/吸附活化能決定的臨界能量，J

Et——顆粒碰撞時相對速度對應的動能，J

e12——碰撞恢復系數

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

m1,m2——顆粒質量，kg

P1,P2——顆粒位置，m

t——時間，s

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

vc——與臨界能量Ec對應的臨界速度，m/s

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

ε——多孔介質孔隙率

下角標

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