颗粒撞击单颗粒覆层的数值计算

谢俊，李晨曦，朱正仁，马昊东，李润东

（沈阳航空航天大学能源与环境学院，辽宁 沈阳 110136）

颗粒碰撞过程广泛存在于生产生活中。如在喷涂领域，汽车喷漆时油漆颗粒与汽车表面或底层油漆颗粒的碰撞过程，类似的如利用陶瓷颗粒制备陶瓷薄膜等[1-2]。在石油化工中，滤床处理油田污水过程中也存在颗粒碰撞现象[3]。在工业生产中，火力发电厂中含尘气流流经换热面时，部分颗粒会停留在换热器的表面形成积灰，影响换热热阻，这会大大影响换热器的工作效率，严重时可能会导致发电厂停产[4]。因此，如何降低化工及工业生产中细颗粒物的危害是大家关注的热点问题。细颗粒物的沉积主要包括三个运动过程：颗粒向壁面运输、颗粒与壁面碰撞以及颗粒的黏附[5-6]。然而，在颗粒的运动过程中，并非所有的颗粒都会黏附在壁面上，在碰撞过后，颗粒也可能会反弹离开表面。颗粒与壁面不同的碰撞结果与颗粒的撞击速度、撞击角度、颗粒粒径以及颗粒的材料性质等因素有关[7]。若一旦生成积灰且不能及时清理就会发生颗粒与颗粒覆层的碰撞现象。

固体颗粒碰撞壁面后发生反弹行为的基本理论通常是基于Hertz[8]和Rayleigh[9]的开创性研究工作。针对不同的现象如塑性变形、黏附等，许多新模型对该基本理论进行了扩展。目前，国内外学者在颗粒碰撞方面已有许多研究，但对黏附和反弹现象的了解仍很有限[10]，特别是对于颗粒撞击覆层的基础研究有待深入。

现在有很多学者对颗粒撞击不同粗糙度的壁面过程进行了研究。王云萍[11]使用有限元软件Abaqus对二氧化硅颗粒撞击壁面的碰撞过程进行了模拟，详细研究了表面粗糙度、入射速度、颗粒粒径及颗粒的弹性模量对碰撞特性的影响。牛佳佳等[12]对粒径为600～800μm的砂粒撞击树脂涂层/铝/合金钢板的过程进行实验分析，研究了砂粒撞击不同壁面材料的反弹特性。基于此，有学者进一步研究了颗粒撞击颗粒覆层的研究。Van Beek 等[13]研究了50μm的玻璃颗粒在不同的碰撞角度、碰撞速度下对薄玻璃颗粒覆层的撞击情况，通过计算入射颗粒的法向恢复系数发现颗粒覆层存在很大的阻尼作用；并采用两体碰撞模型预测了颗粒的反弹角度。Konstandopoulos[14]研究发现，两体碰撞模型在碰撞颗粒具有较高动能时失效，导致预沉积颗粒重排或侵蚀。Abd-Elhady等[15]使用了与Van Beek等[13]相同的设备，用离散单元法代替两体碰撞模型，研究了覆层颗粒的相互作用和颗粒的反弹行为，并对临界捕集速度和反弹速度进行了预测。

颗粒撞击颗粒覆层与颗粒撞击光滑壁面有很大的差别。当撞击表面存在覆层时，覆层的变化情况也会影响颗粒的反弹特性。国内在微米颗粒碰撞过程，特别是撞击覆层的研究尚待深入，同时颗粒接触过程的动力学特性研究也需要进一步探索。因此，本文开展了颗粒撞击覆层的动力学特性研究，分析了颗粒撞击覆层的动力学过程。

1 JKR动力学方程

1.1 JKR理论

目前被广泛使用的各种黏附接触模型的基础为Hertz弹性碰撞模型。但Hertz模型只考虑了两接触颗粒在接触区域内的斥力，而忽略了接触球体间的引力[8]。Johnson 等[16]基于Hertz 模型，考虑接触区内部的接触压力和黏性的影响，如图1所示，利用储存的能量与表面能损失之间的平衡推导出了JKR弹性碰撞模型，得到了接触力FJKR和接触半径a之间的关系，如式(1)。

图1 JKR接触理论中颗粒间接触示意图

式中，接触半径a和颗粒重叠量δ的计算方法参考Chen等的研究[17]。

1.2 JKR动力学方程

根据上述JKR理论，分析两球体在法向碰撞过程中，量纲为1 法向接触力随接触位移的变化过程，具体见图2。两个球体相互接触时，接触位移为0，此时法向接触力由于范德华引力作用瞬间跳至A点。两球体不断靠近过程中，速度逐渐减小，形变随之增大，即接触位移增大，排斥力也同时增加。接触力在B点达到平衡，但由于速度还未达到0，球体继续相向运动直至速度为0。此时，总动能转化为球体的弹性能，接触位移在D点最大，斥力也达到最大。随后，压缩过程结束，卸载过程开始。在点A处，所有的动能都储存在这里，法向位移回到零。

图2 JKR模型的接触力随接触位移的变化曲线[18]

然而，两球体在A点上仍然处于黏附状态。随着球体继续分离，引力将在C点达到最大值。这个点上的力FC在JKR 模型中称为拉脱力（pull-off force）。当对于黏附球体施加的拉伸载荷大于临界载荷FC时，出现接触半径和接触位移突然消失，界面完全脱离，称为硬卸载[19]。当接触位移缓慢变化、无外加载荷的影响时，临界点C之后一段过程仍旧依据JKR 模型变化，直至S点，接触过程结束，称为软卸载。S点为斜率无穷大时的临界点。综上所述，两个球体碰撞后首先是初始点A点，然后到达平衡点B点，随后达到接触位移最大点D点，之后到达硬卸载点C点，最后是软卸载点S点。在本研究中，主要利用软卸载来进行颗粒撞击颗粒覆层的模拟计算。图2中阴影部分面积表示黏附力所做的功，叫作黏附剥离功Wst。Yin[20]和Thornton 等[21]对这部分能量损失进行了量化，如式(2)。

当仅考虑黏附剥离功时，根据能量守恒定律可以求得一个临界捕集速度，如式(3)。

黏附剥离功的存在降低了反弹动能，如果初速度较低或表面能较大，甚至可以完全阻止反弹动能的产生。柳冠青[19]指出一旦发生接触，Wst即已经产生，故也将其视为首次接触的能量损失。其次，黏弹性是固体材料的固有属性，由黏弹性所导致的阻尼作用将会使能量以弹性波的形式耗散。因此，将黏弹性效应与黏附效应相互耦合共同作用于颗粒与覆层的接触过程作为本研究的主要能量耗散途径，结合JKR静态接触模型，引入黏弹性引起的阻尼耗散，可得颗粒与壁面法向碰撞过程的动力学模型如式(4)。

接触力FJKR为JKR模型所描述的黏附和颗粒弹性变形的共同作用。接触力FJKR和接触位移δN的表达式为式(5)、式(6)。

式中，临界卸载力FC、硬卸载点C的接触位移δC以及平衡半径a0的计算方法参考Chen等[17]的研究。

假设固相耗散力与材料在形变中的变化率成正比，则有阻尼力Fnd，如式(7)。

其中，法向阻尼系数可由式(8)求得。

Tsuji的研究表明[22]恢复系数e仅与参数α有关，α可以用一个关于恢复系数e的函数表达式表达[23]。将式(5)和式(7)代入式(4)中得到动力学方程式表达式为式(9)。

式(9)中，K的计算如式(10)。

于是式(9)可以改写为式(11)。

1.3 模型求解方法及验证

常微分方程的求解方法有很多，对于实际问题中归结出的微分方程主要依靠数值解法求解。本文采用四阶龙格-库塔法，对参数δN进行数值求解。程序初始时刻（0 时刻）的接触位移为0，速度为初始入射速度，这部分数值可由实验数据计算获得。计算过程中所取的时间步长为0.1ns。若已知前一时刻i的接触位移和速度，即可通过龙格-库塔方法求得下一时刻i+ 1 的接触位移和速度，以此类推，确定整个碰撞过程的接触位移和速度的动态变化过程。JRK模型考虑了接触过程中的突跳现象，当接触位移小于零时，在黏附力作用下出现颈部而未完全脱离。因此，这里将软卸载点的位移值作为碰撞过程结束的判据。

在碰撞过程中，阻尼系数随着速度变化也时刻发生变化，参数α是关于法向恢复系数en的函数。对于确定的一个碰撞过程，参数α为定值，在初始时刻赋予一个值，通过给定的初始值（接触位移y0、入射速度z0等）计算获得碰撞结束时的反弹速度z，将计算所得的法向恢复系数与同入射条件下的法向恢复系数作对比，若相差过大则调整α的大小重新计算，直至相等或无限接近。根据已知的实验数据确定不同法向恢复系数下对应的α，拟合确定α与en的函数关系式。

根据实验数据确定参数α的表达式后，即可将法向恢复系数决定的α代入反向计算，不断调整入射速度以使计算获得的反弹速度与入射速度之比满足预设的法向恢复系数en，得到入射速度随法向恢复系数的变化关系。

Dahneke 的研究[24-25]详细报道了聚苯乙烯乳胶球以2.5～35m/s的速度冲击抛光石英、金/玻璃、聚苯乙烯表面的入射和反弹速度的测量结果。本文采用1.27μm 的聚苯乙烯颗粒撞击抛光石英表面的实验结果对模型进行验证。表1和表2显示了计算所用的重要参数，对应的JKR 模型关键点如表3所示。

表1 材料参数[26-27]

表2 黏附功与赫兹刚度

表3 JKR模型关键点

根据实验数据，采用上述离散方法，计算了不同法向恢复系数en所对应的参数α，如图3 所示，并对其关系作了最小二乘法拟合，得到拟合关系表达式如式(12)。

图3 参数α与恢复系数en的关系

临界捕集速度是指颗粒与平板表面发生碰撞时颗粒恰好发生反弹而不是黏附在表面上的速度，因此将法向恢复系数en等于0 代入式(12)，可求得颗粒的临界捕集速度。同样地，对于不同的法向恢复系数采用相同方法进行计算，得到恢复系数随入射速度关系变化曲线如图4 所示。从图中可以看出，模型计算结果与实验数据的吻合度较好，模型计算得到颗粒的临界捕集速度约为1.2m/s。因此，将该JKR模型的动力学方程应用于颗粒碰撞的数值计算是可行的。

图4 1.27μm的聚苯乙烯颗粒撞击抛光石英表面实验结果[23]和数值计算结果的比较

2 结果与讨论

2.1 单层颗粒法向碰撞

Abd-Elhady等[28]的研究表明对于矩形排列的颗粒覆层，接触时间与覆层数成正比。如图5 所示，当碰撞在目标颗粒的顶部时，力（波）在一个矩形颗粒排列中传播发生在单个碰撞颗粒链中。基于此，对颗粒撞击单层颗粒的模拟简化处理，分别研究颗粒-壁面和颗粒-颗粒的碰撞过程，将两者联系起来分析颗粒-颗粒(黏附)-壁面的碰撞过程，如图6所示。

图5 入射颗粒与二维颗粒覆层的碰撞[25]

图6 单层颗粒的法向碰撞

根据二氧化硅撞击不锈钢表面的已知物性参数[29]计算了JKR模型下颗粒与颗粒、颗粒与壁面碰撞过程中的关键点，并计算了黏附剥离功（即首次接触的能量损失），具体结果见表4。结果表明，对于颗粒与壁面的碰撞，由于二氧化硅颗粒的界面能与不锈钢相比要小得多，其黏附功会大于颗粒与颗粒之间的碰撞。颗粒与颗粒碰撞的临界捕集速度小于颗粒与壁面之间的碰撞。

表4 颗粒与壁面(或颗粒)间的物性参数[29]

为清楚地了解黏附剥离功在碰撞过程中的影响，避免黏弹性效应的耦合影响，考虑了无阻尼耗散的动力学模型，即式(13)。

数值计算结果如图7 和图8 所示。图7 显示的是颗粒撞击壁面（或颗粒）的接触位移随时间的动态变化曲线。颗粒发生碰撞后，在入射阶段接触位移逐渐增大，直至最大接触位移后开始减小，即进入反弹阶段，最大接触位移随着入射速度的增加而增加。在0.031m/s 和0.014m/s 时分别在颗粒与壁面、颗粒与颗粒间碰撞时发生黏附。当入射速度稍大于黏附速度时，接触位移停止随接触时间的变化，即发生反弹。对比黏附剥离功的理论计算结果，可以发现数值计算的临界捕集速度与理论值基本吻合。随着入射速度的增加，碰撞过程的颗粒接触时间不断减小。因此，较高的入射速度能缩短接触时间。当撞击速度低于临界捕集速度时，颗粒由于反弹时的能量不足以支持其完全脱离壁面，在反弹速度达到0以后又开始向壁面（或颗粒）运动。

图7 接触位移随接触时间的变化关系(无阻尼)

图8 碰撞速度随接触位移的变化关系（无阻尼）

图8显示了颗粒的碰撞速度随接触位移的变化过程。图中黑实线表示在临界捕集速度以下时的碰撞速度变化情况，不同于图7，此次对较长接触时间内的变化过程作了绘制。可以发现，低于临界捕集速度入射的颗粒会以一个循环的形式呈现碰撞速度随接触位移的变化过程。值得注意的是，该循环圈的面积会有缩小，这主要是由于数值计算过程中的截断误差所导致。当截断误差不存在时，临界捕集速度以下的颗粒的黏附剥离功一直在系统中循环，由于没有其他能量损耗，将会在接触面上一直往复运动。这种情况在实际过程中是不存在的，因此，添加阻尼耗散使动力学方程更加符合实际。

图9 参数α与恢复系数en的拟合关系

如图10 和图11 所示，计算过程中考虑了黏弹性阻尼的影响，研究了颗粒与壁面（或颗粒）的接触过程中的动态变化，与不考虑阻尼耗散的情况相比呈现出不同的变化情况。在考虑阻尼耗散时，阻尼力与绝对速度成正比，随着接触位移的变化不断损失整体的能量。在第一次反弹阶段速度变为0，接触位移未超过临界位移（软卸载点S）时，颗粒将会由于接触力表现为引力作用又开始向接触面运动，此后由于阻尼耗散的影响，碰撞速度会在更早的位置变为0，颗粒完全脱离壁面的情况将不可能出现。在入射速度大于临界捕集速度时，碰撞速度随接触位移变化，在接触位移达到临界值时，速度仍不为0，颗粒脱离壁面，曲线变化过程也在这一刻停止。当小于临界捕集速度时，碰撞速度随接触位移呈现一个环形变化过程，由于能量耗散，颗粒的最大动能不断减小，环形半径逐渐减小。该过程类似于阻尼振动，在第一次反弹结束未能脱离壁面时，颗粒与接触面直接进入一个靠近-反弹的循环之中，往复运动中动能耗散，最终在平衡位置B点停止运动。颗粒与壁面(或颗粒)的碰撞都表现出类似的特征，但由于接触对象的不同将导致黏附功、有效杨氏模量、有效质量等的不同，影响临界捕集速度。对比发现，颗粒与颗粒碰撞的临界捕集速度为0.394m/s，低于颗粒与壁面的碰撞，这表明颗粒更容易在壁面上发生黏附。相对于仅考虑黏附剥离功的情况，阻尼耗散的存在也使得临界捕集速度增加，黏附更容易发生。对于颗粒与壁面的碰撞情形，在考虑阻尼耗散情况下，临界捕集速度在0.86m/s附近，远高于理论计算值0.032m/s。

图10 接触位移随接触时间的变化关系

图11 碰撞速度随接触位移的变化关系

在颗粒与壁面（或颗粒）研究的基础上，本文进一步讨论了颗粒撞击黏附在壁面上的颗粒的情况，如图6(c)所示。这里作了一定简化处理，对于入射颗粒P1在碰撞过程中仅考虑与黏附颗粒P2的接触作用，忽略壁面影响。对于颗粒P2，考虑了颗粒P1 和壁面的共同作用，并假设两者引起的形变不相互影响。采用前述动力学模型及求解方法求解，结果如图12和图13所示。

图12 接触位移随接触时间的变化关系

图13 碰撞速度随接触位移的变化关系

图12 显示了单颗粒撞击覆层时入射颗粒P1 的接触位移随接触时间的变化情况。在0.394m/s的入射速度以下，颗粒P1的接触位移不会低于-0.386nm（临界软卸载点S），由于阻尼耗散的影响，颗粒往复运动直至在接触面上黏附。而大于临界捕集速度时，入射颗粒P1 在第一次反弹阶段即脱离壁面，且随着入射速度的增大，最大接触位移峰值增加，并向右移，接触时间也随之减少。相应地，在入射颗粒P1 的接触位移变化的同时，黏附颗粒也在发生变化。在0.69m/s以下，颗粒P2在接触过程中位移随接触时间变化，始终不低于临界值-0.73nm，并最终仍旧黏附于壁面上。由于一开始设定颗粒P2 以平衡状态黏附于壁面，假定在初始时刻位移即为平衡位移值0.884nm，当颗粒运动结束时，该值不变。在图13中可以观察到颗粒P2由于阻尼耗散的影响，碰撞速度在平衡位移点从0开始增加并在同一点结束动态变化。当入射速度大于0.7m/s时，颗粒P2 在第一次远离壁面运动时的接触位移超过临界值发生反弹，该黏附颗粒将从壁面脱离。

3 结论

本文基于JKR 静态接触理论，结合阻尼耗散理论，对单颗粒撞击覆层的碰撞过程进行了数值计算。首先对黏附效应和黏弹性效应进行了分析，建立了基于JKR 理论的动力学方程。其次，采取四阶龙格－库塔法对方程离散求解，用Dahneke的实验结果对模型进行验证。最后，简化单颗粒撞击覆层的碰撞过程，分别讨论有无阻尼耗散时的颗粒-壁面和颗粒-颗粒的碰撞过程，综合考虑并分析了颗粒-颗粒(黏附)-壁面的碰撞过程。主要结论如下。

（1）获得了参数α和法向恢复系数en的拟合关系式，将其代入动力学模型，结合Dahneke的实验数据验证了拓展的JKR模型应用于颗粒碰撞过程的有效性。

（2）分别对有无阻尼耗散时的颗粒撞击颗粒和颗粒撞击壁面的情况进行了模拟计算，结果表明存在阻尼耗散时颗粒所需的临界捕集速度更大，且更加符合实际情况。

（3）单颗粒撞击覆层表面的过程简化为颗粒-颗粒(黏附)-壁面的碰撞过程，主要研究了法向碰撞过程的动态变化。在0.394m/s 的入射速度以下，颗粒P1 会直接在接触面上黏附。在0.69m/s 以下，颗粒P2 会仍旧黏附于壁面上。但当入射速度大于0.7m/s时，黏附颗粒P2将从壁面脱离。

符号说明

a—— 接触半径，m

a0—— 平衡半径，m

E*—— 有效弹性模量，Pa

en—— 法向恢复系数

FC—— 临界卸载力，N

FJKR—— 接触力（JKR模型），N

Fnd—— 膜阻尼力，N

m*—— 有效质量，kg

R*—— 有效半径，m

vst—— 临界捕集速度，m/s

Wst—— 黏附剥离功，J

w—— 黏附功，J/m2

y0—— 接触位移，m

z—— 反弹速度，m/s

z0—— 入射速度，m/s

δ—— 接触位移，m

δc—— 硬卸载点C的接触位移，m

δN—— 法向重叠量，m

ηN—— 法向阻尼系数，kg/s

ρp—— 颗粒密度，kg/m3