柴油机工作产生的尾气中含有大量的碳烟颗粒群，这些颗粒群由基本碳粒子生长凝并而成，最终呈现出枝状、链状、团簇状和球状等结构. 为了降低柴油机中排放的碳烟颗粒，减轻对环境的污染，大量的研究者对碳烟的生成和演化过程进行实验研究和数值模拟.

在针对碳烟的实验测量方法中，利用粒径测量仪可以测得碳烟颗粒的质量、浓度分布和粒径分布等数据^[1-3]，利用颗粒采样技术加上高分辨率透射电子显微镜（HRTEM）能够获得碳烟颗粒的形貌信息^[4-6]. 实验方法只能测量有限的点和物理量，难以同时捕捉单个碳烟颗粒的团聚运动及碳烟颗粒群演化过程的问题. 为了解决该问题，研究者们开始运用分区法、蒙特卡罗法、朗之万动力学等数值方法，考虑初始颗粒的受力情况，建立碳烟颗粒成核、碰撞、表面生长、聚合、凝并、破碎^[7]等动力学仿真模型来进一步研究.

Moran等^[8]发现，利用蒙特卡罗方法模拟同时聚集和表面生长时，表面生长减小了初级颗粒和团聚体的几何标准差，但该方法的精度受到颗粒动力学简化和时间步长^[9]的限制. Ho等^[10]在欧拉拉格朗日方法的框架下提出颗粒的随机碰撞模型，考虑颗粒所受气流阻力、范德华力和重力，研究发现气流参数、颗粒粒径及数密度对颗粒碰撞有较大的影响.

当采用实验测量的方法研究碳烟颗粒的团聚运动时，不容易观测到一定微观区域内碳烟颗粒的生长变化过程. 在对碳烟的仿真建模研究中，难以获得准确的初始条件和参数，对碳烟团聚体的具体形貌和分形维数关注不够.

为了解决上述问题，本研究以一台直列式四缸四冲程柴油机为原型，应用三维CFD软件建立三维仿真模型. 考虑柴油机内碳烟生成过程质量和数密度的变化趋势，对碳烟颗粒进行分析，但是三维碳烟计算模型对碳烟颗粒的形状描述仅限于颗粒直径和总体积，不能体现碳烟颗粒及团聚体的具体形貌特征. 基于三维模型计算得到颗粒分布信息与气流环境参数，以单个碳烟颗粒为基础，运用颗粒动力学方法对碳烟颗粒的动力学演变过程进行模拟和分析，得到颗粒的形貌和分形维数，与实验进行对比.

三维仿真过程采用颗粒粒径模型（particulate size mimic, PSM）^[11-13]来构建详细碳烟模型，并耦合详细的化学反应机理^[14]，对碳烟的成核和生长过程进行模拟. PSM模型将颗粒尺寸分布离散为有限数量的分段，每个分段代表一个特定的碳烟颗粒体积区间，由相同体积的粒子组成^[15]，分段边界定义如下：

(1)$ \left.\begin{array}{l}{V}_{1,\mathrm{min}}={V}_{{\mathrm{MIN}}},\\ {V}_{1,\mathrm{max}}={V}_{{\mathrm{MAX}}},\\ {V}_{i,\mathrm{min}}={V}_{i-1,\mathrm{max}} \left(i > 1\right),\\ {V}_{i,\mathrm{max}}=({V}_{{\mathrm{MIN}}+}{V}_{{{\mathrm{c}}}{2}}){\left(\dfrac{{V}_{{\mathrm{MAX}}}}{{V}_{{\mathrm{MIN}}}+{V}_{{{\mathrm{c}}}{2}}}\right)}^{\frac{i-1}{{i}_{\mathrm{max}}-1}}.\end{array} \right\}$

式中：${i_{\max }}$为分段的数量；${V_{{\mathrm{MIN}}}}$为第1个成核粒子的体积，对应2个芘分子所占的体积； ${V_{{\mathrm{MAX}}}}$为最大的粒子体积；$ V_{{\mathrm{c}} 2}$为含有2个碳原子的碳氢化合物的体积；${V_{i,\min }}$和${V_{i,\max }}$为各分段左、右两端的颗粒体积. 不同分段的尺寸随着几何级的增加而增加.

碳烟颗粒的生长及氧化用以下方程表示：

(2)$ {Q}_{i}={Q}_{{\mathrm{inc}},i}+{Q}_{{\mathrm{coag}},i}+{Q}_{{\mathrm{sg}},i}+{Q}_{{\mathrm{ox}},i} .$

式中： ${Q_i}$为i分段内的体积分数源项，下标inc、coag、sg和ox分别表示碳烟颗粒成核、凝结、表面生长和氧化过程，${Q_{{\mathrm{inc}},i}}$、${Q_{{\mathrm{coag}},i}}$、${Q_{{\mathrm{sg}},i}}$、${Q_{{\mathrm{ox}},i}}$分别表示i分段内的成核、凝结、表面生长和氧化源项. 碳烟成核由下式计算：

(3)$ {Q_{{\mathrm{inc}},j}} = 2{V_{{\mathrm{PAH}}}}{\beta _{{\mathrm{PAH}},{\mathrm{PAH}}}}N_{{\mathrm{PAH}}}^2 . $

式中：${V_{{\mathrm{PAH}}}}$为PAH组分的体积，$ {\beta _{{\mathrm{PAH}},{\mathrm{PAH}}}} $为2个PAH组分相互碰撞的碰撞系数，${N_{{\mathrm{PAH}}}}$为PAH的数量密度.

碳烟表面反应过程中的表面生长体积分数和氧化体积分数通过Marchal方程^[16]计算：

(4)$ {Q_{{\mathrm{sg}},i}} = \alpha {V_{{{\mathrm{c}}2}}}{K_{{\mathrm{sg}}}}{\int_{{V_{i,\min }}}^{{V_{i,\max }}} {\left( {\frac{V}{{{V_{{{\mathrm{c}}2}}}}}} \right)} ^{{\theta }/{3}}}{n_i}(V){\mathrm{d}}V , $

(5)$ {Q_{{\mathrm{ox}},i}} = \alpha {V_{{{\mathrm{c}}2}}}({K_{{{\mathrm{O}}_2}}}+{K_{{\mathrm{OH}}}}){\int_{{V_{i,\min }}}^{{V_{i,\max }}} {\left( {\frac{V}{{{V_{{{\mathrm{c}}2}}}}}} \right)} ^{{\theta }/{3}}}{n_i}(V){\mathrm{d}}V . $

式中：$ \alpha $为反应的表面位置分数，K_sg 为碳烟颗粒表面生长过程的反应动力学速率，n_i(V)为第i段碳烟颗粒粒子的体积分数，$K_{{\mathrm{O}}_2} $和K_OH分别为氧化过程中O₂和OH⁻的反应动力学速率，$\theta $为表示碳烟颗粒不规则性质的表面指数因子.

利用详细的碳烟模型计算得到的单个碳烟颗粒数量太多，且伴随气流湍流运动、温度的变化，计算每一个颗粒的完整生长过程需要巨大的计算资源和时间. 为了保证对颗粒运动模拟的可行性，仅对成核碳烟颗粒的团聚过程逐个进行模拟计算，作出如下假设.

1）计算域中的气流场是均质各向同性的湍流场.

2）颗粒仅受到重力、范德华力及气流对其的阻力.

3）模拟中的颗粒为刚性小球.

本研究的计算流程图如图1所示. 将三维仿真模型中得到的气流参数和碳烟颗粒参数输入颗粒动力学模型，按计算流程将粒子编号并逐个进行计算，当每个粒子（或凝并后的粒子）都计算完时，一个计算循环完成. 一个循环后，凝并后的新粒子和未凝并的粒子都将进入下一次循环，直到整个计算过程结束.

在颗粒动力学的计算过程中考虑了碰撞效率，采用欧拉-拉格朗日方法建立动力学模型，将颗粒视为离散相，气体视为连续相，追踪每一个颗粒的位置、速度、湍动能、湍流耗散率等参数. 其中，颗粒运动轨迹、流体流动的具体计算方法参见文献[10].

实际运动中大颗粒产生的气流会使小颗粒偏离原始的运动轨迹，因此实际碰撞的颗粒数量将小于理论碰撞的颗粒数量，实际碰撞次数与理论碰撞次数的比值称为碰撞效率^[17]. 判断2个颗粒是否发生碰撞，在一个计算域内，一个计算时间步长内颗粒的碰撞效率公式如下：

(6)$ {P_{{\mathrm{coll}}}} = (\pi /4){({d_{{\mathrm{p}}1}}+{d_{{\mathrm{p}}2}})^2}|{{\boldsymbol{u}}_{{\mathrm{p}}1}} - {{\boldsymbol{u}}_{{\mathrm{p}}2}}|{N_{\mathrm{p}}}\Delta t . $

式中：d_p1和d_p2为碰撞颗粒的粒径，${\text{|}}{{\boldsymbol{u}}_{{\mathrm{p}}1}} - {{\boldsymbol{u}}_{{\mathrm{p}}2}}|$为2个颗粒相对速度的绝对值，N_p为计算域内颗粒的总数量，Δt 为时间步长.

计算开始时颗粒随机生成在计算域空间内，三维仿真模型模拟结束获得参数后，根据图1所示的流程计算所有的粒子顺序，粒子是否凝并由临界速度v_cr判断^[18]，公式如下：

(7)$ {v_{{\mathrm{cr}}}} = \frac{1}{{{d_{\mathrm{p}}}}}\frac{{(1 - {e^2})}}{{{e^2}}}\frac{A}{{\pi z_0^2\sqrt {6{p_{{\mathrm{pl}}}}{\rho _{\mathrm{p}}}} }}. $

式中： e为能量恢复系数，A为Hamaker常数，ρ_p为颗粒密度，${z_0}$为颗粒的接触距离， p_pl为材料的极限接触压力. 当2个颗粒碰撞的相对速度大于临界速度时，两颗粒仅发生碰撞；反之，则凝并成新的颗粒.

以YC4FA 116-40型水冷增压柴油机为基础建立数值模型，主要的技术参数见表1. 柴油机采用7孔喷油器，使用1/7燃烧室作为模型的几何结构，提高计算效率. 对温度梯度和速度梯度进行网格加密，对壁面和喷雾区域进行局部网格加密，提高计算精度，总网格数为152 454，活塞处于上止点时的燃烧室网格划分示意图如图2所示. 仿真选用的主要计算模型如表2所示.

模型验证采用与实验^[19]相同的工况（50%负荷）. 图3给出数值模拟与实验的缸压及放热率曲线. 图中，p为缸内压力，H_r为瞬时放热率，$\theta_{{\mathrm{ca}}} $为曲轴转角. 缸压的最大误差为3%，模拟和实验中所用的喷油规律有所差异，因此放热率误差略大于缸压误差，但误差在可接受范围内，说明构建的数值模型具有较好的准确性.

选取与实验^[19]相同的初始条件，从进气门关闭时刻（−130°）到排气门开启时刻（126 °），对10%、30%、70%、90%负荷下柴油机缸内的碳烟生成过程进行模拟. 图4给出排气门开启时不同负荷下缸内碳烟质量的分布情况. 图中，w_B为碳烟质量分数. 从图4可以看出，随着负荷的增加，喷油量上升导致碳烟的生成质量逐渐升高，但每种工况下碳烟颗粒都集中在气缸中部，气缸底部和气缸壁附近的区域碳烟生成量较少. 在中高负荷下，靠近气缸顶部的区域碳烟质量分数较高.

图5给出排气门开启时不同负荷下气缸顶部的碳烟颗粒数密度N分布图. 可以看出，低负荷下气缸顶部的碳烟颗粒平均数密度较小，碳烟颗粒密度随着负荷的增加而增加. 当负荷超过70%时，碳烟颗粒数密度快速上升且中心区域数密度最高，这是因为喷油量增大导致油气混合不均，离喷嘴越近，局部缺氧情况越严重. 在碳烟颗粒形成的过程中，碳烟成核主要集中在燃烧刚结束时，此时多环芳香烃生成的数量较多，90%负荷下的碳烟数密度峰值约为7.23×10²⁰个/m³，其他工况下的局部碳烟数密度峰值均超过3.5×10²⁰个/m³.

图6给出50%负荷下缸内碳烟颗粒数密度和总质量m的变化曲线. 可以看出，后燃期碳烟数密度迅速增大，当数密度达到最大时，碳烟的总质量快速增大. 图7给出50%负荷下碳烟颗粒成核、表面生长、凝结和氧化对质量的影响. 可以看出，燃烧刚结束时碳烟颗粒的生成速率远大于氧化速率，但后续的氧化速率逐步增大，远高于表面生长和凝结速率，成核过程基本结束后碳烟颗粒的总质量不断减小.

图8给出50%负荷下不同时刻缸内碳烟颗粒的总质量分数. 可以看出，受气流的影响，碳烟主要在气缸中部区域生成，气缸顶部、气缸壁及活塞处的碳烟质量分数均较小.

通过详细的碳烟模型，可以得到碳烟形成中各过程质量及数密度的变化情况，但对于碳烟颗粒的具体生成过程不够详细，在该计算结果上对单个碳烟的凝并过程进行模拟. 考虑到碳烟颗粒的采样实验往往在排气管中进行，为了尽量接近排气管中气流的环境参数及团聚体的状态，选取排气门开启时刻气缸顶部的气流参数作为颗粒动力学模拟的环境参数，此时排气门已打开，气缸上平面与排气歧管相连通，气缸上平面的气流环境最接近实验工况^[19].

按照上述方法，对10%、30%、70%、90%负荷下（对应的平均有效压力为0.18、0.38、0.88、1.13 MPa，平均有效压力用p表示）的碳烟质量和密度进行仿真和分析，获得不同工况下的颗粒参数和气流参数. 5种负荷下碳烟颗粒的平均直径分别为20.01、24.27、23.87、31.49、49.26 nm，得到的部分气流参数(气流速度v_g、湍动能k、湍流耗散率ɛ)如表3所示.

在不同负荷下，取排气门开启时气缸顶部计算得到的气流参数作为计算域内碳烟颗粒的环境参数，对碳烟颗粒的凝并过程进行计算. 在气缸顶部选取边长为5 μm的正六面体作为计算域. 由于颗粒数密度越大，碰撞效率越高^[20]，结合计算所得的碳烟颗粒数密度，计算域内数密度取4×10²⁰个/m³，共5×10⁴个颗粒. 计算时长为2.5×10⁻⁴ s，时间步长Δt = 5×10⁻⁶ s，其余的计算参数如表4所示. 表中，${u_{\text{f}}} $为摩擦系数.

计算域中颗粒的初始位置由随机数产生，仿真得到的平均碳粒子直径较小. 结合实验^[21]所得的粒径分布，采用粒子直径d_p = 25 nm的颗粒作为模拟对象，计算域内颗粒的初始空间分布如图9所示. 可以看出，颗粒位置基本符合均匀分布.

由于计算域内颗粒较多，单个颗粒相对于计算域太小，不便于观察，在计算域内取距离底面1 μm、厚0.5 μm的薄层截面域进行观察. 如图10、11所示为在50%负荷下计算开始和结束时计算域的截面图. 从图10可以看出，计算开始时基本粒子均匀分布在计算域中. 如图11所示为计算结束时计算域截面内的碳烟颗粒分布. 对比图10、11可知，大部分碳烟颗粒经过动力学过程后凝并(见图11中方框内的颗粒)，其中完全未发生凝并(见图11中圆圈内的颗粒)的单个颗粒数量较少，由几十个基本粒子凝并成大积聚体的数量不多，由几个或者十几个基本粒子组成的团聚体占总团聚体的比例大.

为了进一步观察碳烟颗粒由基本粒子凝并形成积聚态的具体形貌，在计算结束时随机选取几个包含基本碳粒子数较多的积聚体，将其放大如图12所示. 将图12与实验工况下在排气管内所采样的碳烟形貌(见图13)进行比较^[21]. 可以看出，模拟计算所得的颗粒呈枝状与簇状，与实际采样的形貌相似.

实验测量中的碳粒子数较多，基本碳粒子直径较小，模拟选用的粒径相对较大，且是积聚态，因此具有分型特征. 文献[22]中提供了一种可以在二维层面分析分形维数D_f的方法，该方法引入重叠系数来减少因采样角度不同导致空间中不同团聚体的重叠误差，而动力学模拟计算的结果可以在三维空间内单独显示每一个团聚体，以此避免不同团聚体的重叠误差，因此在利用文献[22]中分形维数D_f的分析方法时得到的结果误差更小. 在每种工况下的计算结果中随机选取30个团聚体计算分形维数后取均值，并与实验结果对比，如图14所示. 随着负荷（等效为平均有效压力p_e）的增加，喷油量和进气速度随之增加，产生的基本碳粒子更多，碰撞频率更高，形成更大的团聚体，且碳烟颗粒的分型维数增大，分形维数为1.68~1.89，与文献[21]中1 800 r/min时不同负荷下在排气管内（采样点1）所测的分型维数趋势相同且数值相近，误差约为1.5%. 考虑到气缸顶部和排气管内颗粒尺寸和气流环境有差异，此误差在可接受范围内. 通过模拟所得的颗粒形貌、分形维数的结果和实验结果具有较高的一致性，验证了颗粒动力学模型的准确性.

（1）燃烧结束时碳烟颗粒的生成速率远大于氧化速率，但后续的氧化速率逐步增大. 气缸顶部、气缸壁及活塞处的碳烟质量分数均较小，碳烟颗粒的氧化作用越强，氧化时间越长，则碳烟颗粒的总质量越低.

（2）在不同负荷下，利用构建的颗粒动力学模型，能够计算得到碳烟颗粒的凝并过程. 所得的碳烟颗粒积聚体形貌均呈枝状和簇状，与实验结果相似.

（3）与实验相比，随着负荷的增加，计算得到的分型维数在低负荷时略大，在高负荷时略小，整体上呈上涨趋势，为1.68~1.89，最大误差约为1.5%.