反应传递多尺度耦合的拟颗粒模拟

葛蔚，李成祥，陈飞国

（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