由于放热量大，因此铝(Al)粉常作为助燃剂添加于火炸药中，而单质Al在燃烧过程中会发生团聚现象，进而降低铝的燃烧效率[1]。Al粉和聚四氟乙烯(PTFE)粉末通过高能球磨工艺制备的机械活化含能材料可以解决团聚问题，该球磨工艺大多为湿法，对于Al-PTFE复合粉体在机械活化过程中的变形机制的研究，有助于改良球磨工艺以及球磨机结构[2-4]，因此对球磨过程中粉体变形机制的研究具有深刻意义。陶俊等[5]通过自制高能球磨机制备了机械活化含能材料，并对其微观性能表征，解释了微米级的Al粉和PTFE粉末表面结合的过程，且高能球磨带来的仅仅是Al-PTEE复合材料的微观物理结构变化，但该方法无法追踪整个粉体的变形过程，且成本较高。相比于实验观察，仿真分析在记录各时间点的磨球及粉体状态上有明显优势，另有学者通过离散元仿真方法来模拟磨球运动规律[6]，分析磨球碰撞过程中影响球磨效率的主要因素，对于湿法球磨，可通过流体力学-离散元(computational fluid dynamics-discrete element method，CFD-DEM)耦合方法来模拟液体环境[7]。
在球磨细微粉体仿真建模过程中，由于粉体与磨球、球磨机尺度跨度过大会导致无法计算，因此目前尚未开发出成熟的模型将微米级粉体与磨球、球磨机置于同尺度下进行仿真。李旭等[8]在行星式球磨机微米级粉体的建模过程中将粉体模型进行放大处理，对于磨球冲击粉体时的相对速度分析了不同速度区间内的碰撞次数分布频率，并用威布尔分布表征分布频率，而实际球磨过程中相对速度的法向与切向分量对Al-PTEE复合粉体的作用效果不同，且并未进行逐一分析。此外，通过磨球对粉体的冲击能量分析了粉体的破碎率，但该方法并不适用于湿法球磨，粉体的体积放大会影响流体对粉体的作用力，进而改变原尺寸粉体的运动规律。
对于湿法球磨，Beinert等[9-10]提出了一种多尺度建模的思想，利用CFD-DEM耦合方法将磨球、粉体分别建立不同尺度下的模型，以解决尺度跨度过大的问题，在介观仿真中建立磨球与粉体模型，分析了不同磨球碰撞相对速度下的粉体位置变化，进而计算出粉体的应力概率。该研究基于磨球与粉体尺度差异较小的前提条件下进行建模，对于不同相对速度建立了8种分析模型，其中包含2种不同法向相对速度以及4种不同切向相对速度的组合，分析得到相对速度切向分量是影响粉体应力概率的主要因素之一[9]。该建模方法虽不适用于磨球与粉体体积差异过大的工况，但多尺度的建模思想在构建应力概率分析模型中具有良好的指导意义。在球磨仿真分析过程中，粉体的应力概率是粉体破碎率计算或后续建立微观粉体变形与破碎分析模型的数值基础[11]，因此对于粉体应力概率分析的需求十分迫切。本文中借鉴了多尺度建模的思想，通过流体体积法(volume of fluid)-离散元法(VOF-DEM)对磨球碰撞时的相对速度进行统计与分析，建立稠密颗粒两相流分析模型确定粉体在球磨过程中的分布区域，最终建立粉体应力概率分析模型，分析粉体应力概率与磨球碰撞相对速度的切向、法向分量的关系。
1 基于CFD-DEM耦合的高能搅拌球磨机宏观建模与分析
1.1 多尺度建模方法
高能球磨多尺度建模是由宏观到微观的过程。宏观仿真建立磨球与球磨机模型，通过CFD-DEM耦合模拟得到磨球碰撞前的相对速度分布，相对速度分布用于计算应力能量分布，基于接触开始时的相对速度应力能量近似为动能[9,11]。
在建立介观模型前，需要确定粉末在流体中的分布区域。介观仿真基于稠密离散相模型(dense discrete phase model，DDPM)建立液体环境下磨球撞击粉体模型，利用不同相对速度下的模拟，得到粉体应力概率与磨球相对速度之间的联系，通过碰撞应力能量分布与粉体应力概率计算出作用于单个粉体上的应力能量分布。
微观仿真通过离散元方法建立粉体的聚集体结构，通过大量小颗粒堆积形成粉体模型，颗粒与颗粒之间通过粘结(bond)键连接。在粘结颗粒模型两端施加固定约束迫使部分粘结键断裂，进而模拟粉体变形与破碎过程，粘结键断裂吸收能量与前文得到粉体应力能量分布对比，最终实现粉体微观形貌分布的预测。
1.2 高能搅拌球磨宏观仿真模型建立
1.2.1 立式搅拌球磨机模型构建与材料参数
通过EDEM-FLUENT耦合对高能搅拌球磨机工况进行数值模拟，仿真模型基于学者陶俊自制立式高能搅拌球磨机，搅拌器为由上到下5层棒式结构，其中内衬与搅拌棒材质均为氧化锆，内部高度为191 mm，底面半径为101.5 mm。磨球选用直径为3 mm的氧化锆(ZrO2)球，ZrO2参数如表1所示。实验中添加磨球的总质量为3 kg、Al粉与PTFE粉末总质量为0.6 kg，其中Al粉与PTFE质量比为6∶4，粒径为5 μm，研磨助液为1.2 L的无水乙醇。
表1 ZrO2球参数
Tab.1 Parameters of ZrO2
密度ρ/(kg·m-3)泊松比μ杨氏模量E/Pa碰撞恢复系数e静摩擦系数μs滚动摩擦系数μk6 0500.32.1×10110.50.50.01
球磨机的工作容积为3 L，在工作时部分搅拌棒与液体无接触，模型涉及气液固三相流耦合导致计算量较大，为减少计算域，加快计算速度，将球磨机上部部分空间以及3层不与液体直接接触的搅拌棒去除，流体计算区域上边界置于球磨机工作状态时的自由液面以上，简化后模型能够满足创造自由液面环境，可确保仿真结果与实际情况一致，简化后导入EDEM中，模型如图1所示。
图1 立式搅拌球磨机模型
Fig.1 Model of vertical stirring ball mill
1.2.2 CFD-DEM耦合计算设置
本研究中采用Hertz-Mindlin接触模型，在该模型中法向力分量是基于Hertzian的接触理论，切向力模型基于Mindlin的研究提出，该模型能够准确且高效的计算出颗粒之间的接触力[12]。在设置颗粒工厂后，开始生成直径为3 mm的球形颗粒，待颗粒静置后设置初始时间为0，搅拌器转速设置为300 r/min，仿真时间步设置为固定步长，数值为瑞利时间步长的35.9%，设置颗粒碰撞追踪后等待耦合接口接入。
CFD-DEM耦合接口采用学者吴亮[13]开发的VOF-DEM模型接口，即含自由液面的气体液体固体三相流耦合，并通过颗粒沉降实验验证了VOF-DEM模型能够正确计算颗粒-流体相互作用力。该接口将VOF模型扩展到用Reynolds湍流模型来模拟湍流流动，并提出了一种新的虚拟双网格空隙度模型，用于描述流体与颗粒间的相间体积置换。空隙率计算公式为
εf,i=1-∑iφj,i，
(1)
式中，φj,i为颗粒j的体积对实际计算网格i的贡献值。
曳力计算公式为
(2)
式中：εf为流体局部空隙率；u为流体速度，m/s；v为颗粒速度，m/s；Vp为颗粒体积，m3。
系数β表达式为
(3)
式中：ρf为流体密度，kg/m3；dp为颗粒直径，m。
曳力系数CD表达式为
(4)
式中，Rep为颗粒的雷诺数[13]，表达式为
(5)
宏观仿真建立球磨机、磨球以及流体仿真模型，不添加粉体，而在实际球磨过程中微米级粉体会改变流体力学性能和动力学性能，进而影响磨球的运动规律[14-15]。为了避免球磨机空转情况，减小误差，通过实验测试结果对流体参数进行修正，混有复合粉体的乙醇悬浮液密度和黏度测试结果如表2所示，该参数是影响流体力学性能和动力学性能的最主要因素，进而影响流体对颗粒的曳力、升力等作用力。
表2 复合粉体悬浮液的参数
Tab.2 Parameters of composite material suspension
材料密度ρ/(kg·m-3)黏度υ/(MPa·s)无水乙醇 7901 200复合粉体-乙醇悬浮液9693 580
曳力模型选择Ergun&Wen-Yu模型，多相流模型选择VOF模型隐式算法，主相对应乙醇，第2相对应空气，湍流模型选择线性湍流模型。在流体域设置中选择动坐标系方法来计算网格运动，同时加载UDF源项与多孔介质域模型，通过patch功能对计算域内的乙醇和空气体积分数进行设置。时间步长设置为5×10-6 s，总计算时间设置为0.7 s，开始耦合计算。Kim等[17]在立式搅拌球磨实验中设置搅拌轴转速为300 r/min，运行至时间为0.16～0.2 s时磨球高度几乎不变，说明立式搅拌球磨机可短时间内到达稳定或准稳定状态，此外，李旭等[8]在行星式球磨机仿真中观察到磨球在0.7 s达到稳定状态，同样对上述观点进行验证。
1.3 宏观仿真磨球与流体运动规律分析
图2所示为流体和球磨机内壁对单个磨球沿z轴方向的作用力的平均值随时间变化曲线，其中磨球与内壁接触产生的作用力包括底面对磨球沿z轴的法向力以及竖直壁面对磨球沿z轴的切向力，图3所示为0.05 s时流体对磨球作用力矢量分布图。由图2、3可知，在0～0.17 s流体对磨球作用力法向沿z轴向下，由于初始时刻流体产生的离心力使磨球向球磨机壁面处运动，同时磨球不断向上堆积，上部磨球沿z轴瞬时速度超过了流体上升速度，受到流体的作用力整体向下，而下方磨球沿z轴瞬时速度小于流体上升速度，受到流体的作用力整体向上，但整体磨球受到流体向下作用力大于向上作用力，因此该时间段内流体对磨球沿z轴作用力向下。将所有磨球视作整体的磨球运动系统，该系统在z轴方向受重力、流体作用力及球磨机内壁提供作用力，每个磨球的重力为8.55×10-4 N，在0～0.2 s时间段内，流体与球磨机内壁对磨球作用力矢量和大于磨球自身重力，磨球开始向上运动。在0.2～0.7 s时间段内，流体作用力趋于稳定，而磨球与球磨机内壁接触力波动较大，取球磨机内壁与磨球接触力平均值为4.85×10-4 N，流体作用力波动区间为(3.40～3.80)×10-4 N，二者矢量和分布区间约为(8.25～8.65)×10-4 N，与磨球重力数值接近，可保证磨球系统在该工况下竖直方向受力平衡，磨球高度保持稳定。
图2 磨球所受外力沿z轴方向分量随时间变化
Fig.2 Component of external force on ball along z axis varies with time
图3 0.05 s时流体对磨球作用力矢量分布
Fig.3 Vector distribution of force exerted by fluid on grinding ball at 0.05 s
图4所示为磨球碰撞相对速度法向分量与切向分量平均值随时间变化曲线，图5所示为磨球碰撞力法向分量与切向分量平均值随时间变化曲线，图中曲线变化趋势一致，在初始阶段升高后趋于稳定。磨球整体运动规律是反映球磨是否达到稳定状态的重要标志，在0.7 s时球磨达到稳定状态，满足后续分析磨球在碰撞过程中的相对速度分布规律的条件。
图4 磨球相对速度随时间变化Fig.4 Relative speed of grinding ball varies with time图5 磨球碰撞力随时间变化Fig.5 Impact force of grinding ball varies with time
在构建粉体应力概率分析模型中，在构建粉体应力概率分析模型中，磨球的相对速度法向分量与切向分量会影响二者接近过程中磨球间的流场运动，进而直接影响粉体的应力概率分布，因此对磨球撞击过程中的相对速度分析是必要的。对球磨机运行达到稳定状态后的数据进行分析。图6所示为0.68～0.7 s内磨球碰撞的相对速度法向分量vn频率分布图，横坐标表示相对速度法向分量，将0～1 m/s分为50个速度区间，纵坐标n表示对应碰撞次数、ψ表示碰撞次数累计占比，柱状图呈现的规律为幂律分布，整个速度区间分布广泛，碰撞次数主要集中在相对较低的速度区间，随着相对速度法向分量的增加，碰撞次数迅速减少。其中相对速度法向分量平均值为0.14 m/s，小于平均值的碰撞次数占总体的67%，而相对速度法向分量大于0.32 m/s的碰撞次数占比仅为10%，在概率统计分析中，针对各区间的占比对速度较大区域进行精细划分，对速度分布占比较小区域进行粗略划分。
图6 磨球相对速度法向分量分布规律
Fig.6 Distribution of normal component of grinding ball’s relative velocity
图7所示为相对速度法向分量低于0.14 m/s时各区间内对应的相对速度切向分量vt频率分布图，纵坐标n表示对应碰撞数。由图可知，随着相对速度法向分量的增加，相对速度切向分量分布频率曲线逐渐趋于平缓，且峰值向右移动。说明相对速度法向分量对切向分量的分布存在间接影响，整体大相对速度法向分量区间内对应的相对速度切向分量更大，此外，大相对速度法向分量区间内的切向分量数值分布较为平均。
图7 磨球相对速度切向分量分布规律
Fig.7 Distribution of tangential component of grinding ball’s relative velocity
由于相对速度分布区间较为广泛，且针对各种可能性进行仿真计算量过大，因此将各区间内不同相对速度出现概率进行统计，结果如表3所示。
表3 相对速度频率分布统计
Tab.3 Statistics of relative velocity frequency distribution
切向分量分布区间/(m·s-1)法向分量区间/(m·s-1)[0, 0.04)/(m·s-1)[0.04, 0.08)/(m·s-1)[0.08, 0.14)/(m·s-1)[0.14, 0.22)/(m·s-1)[0.22, 0.32)/(m·s-1)[0.32, 1.00]/(m·s-1)累计[0, 0.20)19.79 12.06 7.533.561.480.9645.37%[0.20, 0.40)4.275.325.904.332.832.0424.70%[0.40, 1.00)2.082.833.883.993.294.8320.89%[1.00, 2.00]1.281.151.371.210.763.269.03%累计27.42%21.36%18.67%13.09%8.36%11.09%100.00%
1.4 粉体在球磨机内的运动与位置分布分析
对于粉体的应力概率分析，需要对磨球碰撞时处于冲击区域内的粉体数量进行计算，由于Al粉和PTFE粉末的密度大于乙醇的，且其在液体中的悬浮特性，在搅拌过程中粉体在流体中的并非均匀分布[16]，因此在进行粉体应力概率分析前，需要对粉体在流体中的分布区域进行分析。
取搅拌棒所处的液体层建立流体区域2D模型，进行稠密颗粒流体两相流分析模型，多相流模块选择欧拉(Eulerian)多相流，并勾选DDPM，选择隐式算法，湍流模型选择标准的k-ε模型。在DDPM模型中，通过作用在颗粒的各种力平衡(在拉格朗日参考系中的)来预测计算离散相的运动轨迹。在牛顿第二定律下控制粒子运动的常微分方程为式(6)：
(6)
(7)
式中：m为粒子的质量，kg；x为粒子的位移，m；t为粒子运动的时间，s；v表示颗粒的速度，m/s；Fd、Fv、Fg与Fo分别表示颗粒所受的曳力、虚拟质量力、重力与其他作用力，N。
为了还原搅拌球磨前粉体在流体中的随机分布，在沿底面半径方向设置离散相初始生成位置为均匀分布，同时设置离散相喷射速度为1×10-6 m/s，使其初始速度接近0，离散相粒径设置为5×10-3 mm，密度设置为2.7×103 kg/m3。搅拌棒的转动通过动网格实现，网格划分为非结构化网格，动网格计算方法使用弹簧光顺法(spring-base smoothing)以及局部网格重划法(local remeshing)，通过编译并导入profile文件来描述搅拌棒的转动。
设置时间步长为1×10-5 s、总计算时间为1 s，得到粉体速度矢量及分布如图8所示。由图可知，靠近转轴中心的粉体逐渐汇入外部粉体的运动路径中，由于微米级粉体密度大于乙醇液体，致使其在离心力的作用下向罐壁周围靠拢，在稳定时粉体的分布区域与磨球近似。在3D模型中计算该区域内流体的体积为6.87×10-4 m3，实验中实际添加粉体颗粒数为3.70×1012，在此基础上计算得到粉体质心之间的距离为5.7×10-3 mm，该参数用于后续建立粉体应力概率分析模型。
图8 粉体运动分布规律
Fig.8 Distribution law of powder motion
2 粉体应力概率分析
2.1 基于DDPM模型建立粉体应力概率分析模型
在建立粉体的应力概率分析模型时，由于磨球与粉体尺度差异大，将磨球模型建立为运动的几何体，粉体为离散相，磨球撞击区域以及周围流体为计算域。建立磨球撞击2D模型，由宏观仿真得到碰撞前磨球质心之间的平均距离为3.355 mm，作为磨球运动初始距离，由于磨球与附近的粉体和流体运动速度基本一致，在建立稠密离散相模型时将流体初始速度设置为0，离散相喷射速度设置为1×10-6 m/s，离散相生成位置间隔设置为5.7×10-3 mm，离散相粒径设置为5×10-3 mm，密度设置为2.7×103 kg/m3，将磨球冲击路径区域及其周围区域设置为离散相生成区域，确保粉体分布区域能够覆盖流场运动变化区域。多相流模型选择为Eulerian多相流模型下的DDPM，磨球的运动通过动网格实现，可以准确地反映磨球在接近过程中的流场变化以及对粉体运动的影响，网格为非结构化网格划分，动网格计算方法使用弹簧光顺法以及局部网格重划法。通过编译profile文件对磨球运动进行描述，设置时间步长为1×10-5 s，在磨球运动至碰撞点之间的间隔为5×10-3 mm时停止，对冲击区域内的受撞击的粉体数量进行统计。
2.2 流场运动与粉体应力概率分析
磨球的相对速度法向分量分别选取0.02、0.04、0.06、0.08、0.10、0.12 m/s，相对速度切向分量选取0.10、0.20、0.30 m/s进行分析，以上共计18种组合，取自磨球碰撞相对速度频率分布内概率集中区域。图9所示为不同相对速度法向分量对流场运动以及粉体应力概率的影响分析，云图绘制取自0.10 m/s切向相对速度下法向相对速度为0.02、0.04、0.06 m/s的仿真结果，流场速度云图如图9(a)、(c)、(e)所示，粉体运动速度分布如图9(b)、(d)、(f)所示。对比流场速度云图可知，相对较高的法向相对速度会导致更加明显的流场速度梯度变化，同时会弱化切向相对速度对流场以及粉体运动规律的影响，进而导致在相对速度切向分量保持不变的情况下，较大相对速度法向分量仿真中受应力作用的粉体数量相对更多。在磨球冲击区域内，以撞击时的接触点为圆心，以8.60×10-2 mm为半径划分出的区域为粉体受应力变形区域，由于磨球接触时该区域高度小于粉体的粒径，导致区域内所有粉体均会受到应力作用，对应力区域内的粉体数量进行统计得到单次碰撞受应力作用的粉体数量。
(a)vn=0.02 m/s、 vt=0.10 m/s时流场速度云图(b)vn=0.02 m/s、 vt=0.10 m/s时颗粒位置分布(c)vn=0.04 m/s、 vt=0.10 m/s时流场速度云图(d)vn=0.04 m/s、 vt=0.10 m/s时颗粒位置分布(e)vn=0.06 m/s、 vt=0.10 m/s时流场速度云图(f)vn=0.06 m/s、 vt=0.10 m/s时颗粒位置分布图9 不同相对速度法向分量下的磨球间流场速度云图与粉体位置分布Fig.9 Velocity cloud diagram of flow field between grinding balls and powder location distribution under different normal component of relative velocity
图10所示为不同相对速度法向分量对流场运动以及粉体应力概率的影响分析，云图绘制取自0.06 m/s法向相对速度下相对速度切向分量为0.10、0.20、0.30 m/s的仿真结果，流场速度云图如图10(a)、(c)、(e)所示，粉体运动速度分布图如图10(b)、(d)、(f)所示。对比流场速度云图，随着相对速度切向分量与法向分量的比值增大，流场的速度分布梯度变化逐渐减弱，说明在相对速度法向分量保持不变的情况下，相对速度切向分量增加导致碰撞接触点两侧的流场速度水平方向分量增加，进而导致受撞击的粉体数量减少。其中图10(e)中的流场速度云图与图9(a)中的规律一致，且二者受撞击的粉体数量相等，由于二者是基于相同相对速度切向与法向分量比值(比值为5)下的仿真结果，因此粉体的应力概率与相对速度切向与法向分量比值存在联系。
(a)vn=0.06 m/s、 vt=0.10 m/s时流场速度云图(b)vn=0.06 m/s、 vt=0.10 m/s时颗粒位置分布(c)vn=0.06 m/s、 vt=0.20 m/s时流场速度云图(d)vn=0.06 m/s、 vt=0.20 m/s时颗粒位置分布(e)vn=0.06 m/s、 vt=0.30 m/s时流场速度云图(f)vn=0.06 m/s、 vt=0.30 m/s时颗粒位置分布图10 不同切向相对速度下的磨球间流场速度云图与粉体位置分布图Fig.10 Velocity cloud diagram of flow field between grinding balls and powder location distribution under different tangential component of relative velocity
对上述18种组合的仿真结果进行分析，其中相对速度切向与法向分量比值分布区间由0.833至15，对于其中数值跨度较大的区域进行补充仿真分析。图11所示为不同相对速度切向与法向分量比值下的受撞击粉体数量分布以及拟合曲线。曲线类型为指数型函数曲线其中判定系数R2为0.93，数值接近于1，表明该拟合曲线的预测效果较为理想。
图11 不同相对速度切向与法向分量比值下受撞击粉体数量分布的拟合曲线
Fig.11 Fitting curve of quantity distribution of impacted powder under different relative velocity tangential and normal component ratios
2.3 相对速度切向与法向分量比值-受撞击粉体数量拟合曲线普遍适用性验证
在前文中进行的粉体应力概率分析中，磨球相对速度法向分量取值区间为0～0.14 m/s，进而得到不同相对速度切向与法向分量比值下受撞击粉体数量分布的拟合曲线。为了验证该拟合曲线普遍适用于各速度区间，在0.14～1.00 m/s速度区间随机选取相对速度法向分量来进行验证。在相对速度选取为vn=0.20 m/s、vt=0.20 m/s时得到的粉体速度分布图像如图12所示。由图观察得到受撞击的粉体数为76。
图12 粉体速度分布
Fig.12 Distribution of powder velocity
分别选取切向与法向相对速度比值为1、2.5、10进行验证，代表低、中、高水平的数值，统计不同相对速度下的受撞击粉体数量，如表4所示。通过拟合曲线对以下3种情况中受撞击粉体数量的预测结果分别为77、38、2，实际仿真测试结果与拟合曲线推测结果数值接近，验证了拟合曲线的普遍适用性。
表4 受撞击粉体数量统计
Tab.4 Statistics of impact powder quantity
相对速度法向分量vn/(m·s-1)相对速度切向分量vt/(m·s-1)受撞击粉体数n0.200.20760.280.70400.404.004
3 结论
1)高能搅拌球磨微米级粉体过程中，通过CFD-DEM耦合建立宏观仿真模型对磨球碰撞中的相对速度法向和切向分量进行统计和分析，发现不同相对速度法向分量区间内的碰撞次数分布呈现幂律分布规律，随相对速度法向分量增大，对应的切向分量整体呈现增大趋势。
2)在稠密颗粒两相流分析中确定粉体在球磨过程中的分布区域，通过单位体积内的粉体数量，计算得到相邻粉体质心之间的距离为5.7×10-3 mm，用于建立粉体应力概率分析模型。
3)基于以上仿真结果，建立粉体应力概率分析模型，发现粉体应力概率与磨球碰撞相对速度切向、法向分量比值存在直接联系，拟合出二者关系曲线为指数型函数曲线，其中判定系数R2为0.93，数值接近于1，说明拟合曲线的预测效果较为理想。拟合曲线结合磨球碰撞中的相对速度分布确定出有效碰撞频率分布，为后续建立粉体塑性变形模型或计算粉体破碎率提供数值基础。
[1]SIPPEL T R, SON S F, GROVEN L J.Altering reactivity of aluminum with selective inclusion of polytetrafluoroethylene through mechanical activation[J].Propellants, Explosives, Pyrotechnics, 2013, 38(2): 286-295.
[2]ADAMENKO N A, KAZUROV A V, SAVIN D V, et al.Study of structure formation in aluminum filled polytetrafluoroethylene after explosive compression[J].Inorganic Materials: Applied Research, 2022, 13(1): 68-74.
[3]OKHLOPKOVA A A, SLEPTSOVA S A, PARNIKOVA A G, et al.Triboengineering and physicomechanical properties of nanocomposites based on PTFE and aluminum oxide[J].Journal of Friction and Wear, 2008, 29(6): 466-469.
[4]DING T, GUO W C, CAO W, et al.Experimental study of reaction properties of aluminum/polytetrafluoroethylene powder under laser ablation[J].AIP Advances, 2021, 11(8): 085010.
[5]陶俊, 王晓峰, 韩仲熙, 等.铝粉/聚四氟乙烯机械活化含能材料的制备及其微观性能研究[J].材料导报, 2018, 32(6): 894-898.
[6]UHLMANN E, EULITZ A, DETHLEFS A.Discrete element modelling of drag finishing[J].Procedia CIRP, 2015, 31: 369-374.
[7]CHANG Q, DI S B, XU J, et al.Direct numerical simulation of turbulent liquid-solid flow in a small-scale stirred tank[J].Chemical Engineering Journal, 2021, 420: 127562.
[8]李旭, 刘彦, 安丰江, 等.行星式球磨颗粒流场分布与形貌变化规律[J].兵工学报, 2022, 43(4): 876-891.
[9]BEINERT S, FRAGNIERE G, SCHILDE C, et al.Multiscale simulation of fine grinding and dispersing processes: stressing probability, stressing energy and resultant breakage rate[J].Advanced Powder Technology, 2018, 29(3): 573-583.
[10]BEINERT S, KWADE A, SCHILDE C.Strategies for multi-scale simulation of fine grinding and dispersing processes: drag coefficient and fracture of fractal aggregates[J].Advanced Powder Technology, 2018, 29(3): 707-718.
[11]BEINERT S, FRAGNIERE G, SCHILDE C, et al.Analysis and modelling of bead contacts in wet-operating stirred media and planetary ball mills with CFD-DEM simulations[J].Chemical Engineering Science, 2015, 134: 648-662.
[12]TSUJI Y, TANAKA T, ISHIDA T.Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[J].Powder Technology, 1992, 71(3): 239-250.
[13]吴亮.含固体颗粒的两相流界面变化的数值研究[D].天津: 天津大学, 2018.
[14]WANG Z K, TENG Y J, LIU M B.A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and hilbert curve based searching strategy[J].Journal of Computational Physics, 2019, 384: 151-169.
[15]BHARADWAJ R, KETTERHAGEN W R, HANCOCK B C.Discrete element simulation study of a freeman powder rheometer[J].Chemical Engineering Science, 2010, 65(21): 5747-5756.
[16]SUN Z N, ZHU J.A four-quadrant flow regime map for two-phase liquid-solids and gas-solids fluidization systems[J].Powder Technology, 2021, 394: 424-438.
[17]KIM S, WOO S C.Analysis of ball movement for research of grinding mechanism of a stirred ball mill with 3D discrete element method[J].Korean Journal of Chemical Engineering, 2008, 25(3): 585-592.
[18]WANG Z, TENG Y, LIU M.A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and hilbert curve based searching strategy[J].Journal of Computational Physics, 2019, 384: 151-169.