冰 Ih 分子动力学对面心立方（FCC）金属的影响：界面熔化现象及抗结冰性能指数 亚历山大·布雷洛夫斯基（Alexander Brailovski）、 阿里·贝杜恩（Ali Beydoun）、 安德烈·格拉（André Guerra）、 亚历杭德罗·D·雷（Alejandro D. Rey） 以及 菲利普·塞尔维奥（Phillip Servio）

《Crystals》：Molecular Dynamics of Ice Ih Impacts on FCC Metals: Interfacial Melting and an Anti-Icing Index of Merit Alexandre Brailovski, Ali Beydoun, André Guerra, Alejandro D. Rey and Phillip Servio

　　摘要 冰在暴露结构上的附着仍然是一个主要的操作挑战，这促使人们寻找基于材料的被动防冰策略。分子动力学提供了一种受控的方法来研究冰与表面的相互作用，超越了实验装置的限制。在这项工作中，我们开发了一个模拟框架，用于模拟固态六角形冰滴对金属基材的影响。冰撞击的模拟涵盖了不同的速度（10-120米/秒）、温度（120-250开尔文）以及面心立方表面材料（金、铜、银、铝和镍）。使用LAMMPS、mW水力场、EAM/合金金属势以及Lennard-Jones水-表面相互作用，我们通过角序参数和准液层测量来量化相变过程，并结合OVITO中的CHILL+算法进行补充。通过隔离所有外部因素，我们发现熔化程度随着速度和温度的增加而增加，并且与基材属性相关：具有高热导率和低杨氏模量的金属倾向于减少碰撞后的冰熔化。前者与后者的比率，即一个衍生的性能指标Υ，与熔化百分比有显著相关性，并确定银是最有效的防冰材料。统计分析强烈表明，这些表面属性影响界面熔化，支持使用这种建模框架来筛选和设计防冰材料。 1. 引言冰附着在飞机机翼和机身表面是一个与航空安全和效率密切相关的工程问题[1,2]。这些结构上的冰积累不仅会造成空气动力学和飞行稳定性的损失，还可能导致重大事故[3,4]。这个问题不仅限于航空领域：其他外部结构也可能受到影响[5]。在冬季，由于冰的积累，太阳能板和风力涡轮机的生产效率会下降[6]，随着向替代能源的过渡，这个问题将变得更加频繁。目前的商业除冰策略主要集中在使用化学液体的表面处理上，其中丙二醇是最广泛应用的剂之一，特别是在航空领域[7]。丙二醇作为冰点降低剂，可以降低冰形成的温度，从而更容易地从表面去除冰和雪[8,9]。然而，尽管它被认为比乙二醇等替代品毒性更低，但其广泛使用仍对环境构成风险，因为它会流入水体并耗尽其氧气储备，降低其生物可用性并影响水生生物[9,10]。除了这些直接的生态影响外，大量除冰液体的运输和应用还会通过燃料消耗、排放和基础设施进一步增加环境成本[11]。鉴于这些挑战，人们已投入大量努力开发基于材料的疏冰表面，以防止冰积累或减少冰附着。近年来，关于防冰表面的研究有所增加[12,13]；然而，其中许多工作仍然停留在实验阶段，常常受到缺乏大气冰条件进行验证的限制。基于计算的研究通常依赖于静态模拟，测试水滴在表面冻结后的冰附着强度[14,15,16,17,18]。这些工作基本上是通过增加基材和水之间的接触角来测试超疏水材料促进冰去除的性能。由于接触角决定了表面的润湿性[19]，表面越疏水，与新形成的冰之间的附着强度就越小[14]。研究这种现象考虑的是液滴落在静止的冷表面上并冻结的情况，这种情况确实经常发生。然而，在飞机或风力涡轮机移动时以及液滴为固态的情况下（例如在冰雹风暴中），也会发生一些冰积累。由于气候变化的影响，近年来冰雹事件变得更加常见，尽管两者之间的联系尚未完全确立[20]。2023年，美国报告了近7000起冰雹事件[21]，据报道在2020年造成了140亿美元的损失[21]。因此，这一现象值得更多关注，我们认为建模是实现这一目标的正确工具。已经尝试使用大气型冰进行实验工作，例如使用超声波悬浮器[22]和放置在上面的冷冻枪[23]来冷冻水滴，然后吹冷 air [24]直到发生二次成核[24]并使液滴固化。然后关闭悬浮器，让液滴在重力的作用下加速直到撞击表面[25]。测试了几种表面，并用高速摄像机分析了撞击动态[25]。不幸的是，这些实验之间存在显著的变异性，因为很难同时控制表面和液滴的温度以及液滴下落和碰撞时的周围环境[26]。风 gusts、气流或悬浮器的不稳定性可能会引入旋转，这通常用高速摄像机难以检测到。温度或压力的任何变化也可能在碰撞发生之前引发提前融化[27]，严重影响结果的重复性。最后，由于高速摄像机的焦平面极其有限，再加上旋转的变化[28]，在合理的时间内获得可重复的实验结果非常具有挑战性。正如预期的那样，旋转会改变液滴的向下轨迹，使其经常脱离摄像机的焦平面。随着液滴高度的增加，这种效应会更加明显，这严重限制了实验工作中可以测试的速度，特别是因为在这种不耗费大量资源的设置中，液滴只能通过重力加速[25]。与之前的几项物理研究相比，模拟可以使用更高、更现实的速度[25,29]，并提供额外的环境控制，消除了许多实验限制。到目前为止，基于动态撞击的模拟工作要么严格使用过冷水[30]，要么测试冰滴是否会导致表面塑性变形[31]。虽然大部分冰积累来自过冷水冻结，但冰雹撞击也会在表面上积累，而这在建模中尚未得到探索。在我们的工作中，我们感兴趣的是开发一个系统，能够评估表面如何防止固态冰滴熔化。这是有用的，因为冰积累通常发生在冬季，表面通常是冷的，熔化的冰往往会重新冻结，从而增加冰的表面积（因此增加附着强度），从而使其更难以去除。一个能够模拟不同材料、速度和温度下的冰-表面相互作用的基于建模的框架可以为评估潜在表面成分的防冰性能提供一种成本效益高且高效的方法。先前的计算和实验研究已经确定了一个在固体冰撞击表面时需要考虑的关键参数：准液层（QLL）[32,33,34,35,36]。这一预熔化层形成在固态液滴的外部，是有序内部（六角晶体）和无序外部（液态或气态）之间的过渡状态[37]。这一层的厚度可以与液体滴的温度相关联，并可用于将实验数据与建模工作联系起来，进而作为模型验证的一部分。本研究的主要目标是（1）开发一个分子动力学框架，用于模拟固态六角形冰（ice Ih）液滴对金属基材的撞击；（2）确定冰的相变过程，并量化液滴在速度范围为10至120米/秒和系统温度范围为120至250开尔文的情况下与五种不同FCC金属表面的能量传递和温度变化；（3）统计评估材料属性对防冰性能和碰撞后冰熔化的影响；（4）通过实验验证模型。本文的组织结构如下。我们首先描述了用于设计液滴-基材系统的计算技术及其相关的力场。然后，我们详细讨论了温度、速度和金属类型对液滴在碰撞前、碰撞中和碰撞后的相变和能量平衡的影响。接着，我们通过对测试材料的机械和热属性进行统计分析，并将其与文献数据相关联，从而将结果与材料属性联系起来。最后，我们通过形态学和液滴结构比较以及准液层（QLL）的分析将模型与实验数据联系起来。 2. 材料和方法这些模拟使用了LAMMPS（大尺度原子/分子大规模并行模拟器，版本1.12.0）[38,39]和Moltemplate（版本2.22.5）[40]软件，可视化使用OVITO（版本3.15.3）[41]完成。x和y方向的边界条件采用收缩包裹，而z方向采用固定边界条件。选择收缩包裹是因为它允许横向盒子尺寸适应瞬时粒子位置，避免了半无限周期性单元的计算要求，同时允许在横向应力上进行微调[42]。这种动态调整有助于在这种固-固系统中稳定压力，因为系统不受严格固定边界的限制。文献也支持在这种系统中的使用[30]。模拟盒体的尺寸位于垂直于碰撞的平面上。金属基材的厚度为，包含80,200个原子，采用面心立方（FCC）晶格排列。FCC结构的总共有1%的空位，这对于这些金属来说是一个较高但合理的估计[43,44]。基材的晶格常数根据模拟的金属而变化，这会导致厚度的微小变化。液滴的中心距离基材顶部，由ice Ih组成。冰在x方向上的晶格常数为，在y方向上为，在z方向上为[45]，液滴的直径为。因此，这个液滴大约有3470个单元格，每个单元格包含8个分子，总计27,606个水分子。这提供了足够的样本量来得出有意义的结论[46]，而不需要巨大的计算资源。由于六角形冰的空位比例变化很大[47]，我们选择1%的空位作为一个合理的估计[48]。引入冰和表面缺陷是为了模拟不完美的真实晶结构，这些缺陷在分子动力学中经常被使用[49]。本研究中使用的模型系统及其组分的微观结构如图1所示。图1。（a）带有冰滴（蓝色）和金属表面（绿色）的系统模型，（b）ice Ih微观结构（蓝色）和H键（绿色），以及（c）金属表面微观结构（FCC，绿色）。 2.1. 水模型广泛使用的原子模型，如TIP4P/Ice [50]和TIP5P [51]在再现实验冰的许多属性（如熔点温度和密度）方面相对准确。然而，它们在计算上代价高昂，无法准确再现某些属性，如冰的热容量[52,53]或粘度[54]，需要校正。为了解决这些问题，可以使用粗粒度单原子水（mW）模型，因为它可以将计算时间减少多达99%[55]，同时与原子水模型一样准确甚至更准确[55,56]。虽然mW水在预测实验水的扩散系数方面不够准确（因为缺乏氢原子），但在其他参数（如熔点、密度等）方面表现更好[55]。然而，由于这项工作不涉及扩散问题，因此系统中不需要精确描绘这一参数。Stillinger–Weber（SW）势用于硅模型，它是mW模型的基础。像水一样，硅在常压下形成四面体晶体[57]。硅和水之间的相似性不仅限于此：它们的液态密度都比固态密度高[58]，并且它们的过冷液密度随着压力的增加而增加，这是一种异常行为[59]。原子水模型中存在的长程静电相互作用在mW水中不存在，但Molinero和Moore认为存在的短程相互作用足以准确表示实验数据中的异常行为[55]。SW势作为粒子对之间距离和粒子三元组之间角度的函数可以表示为[55,60] (1) (2) (3) 其中。变量（介电常数）、（粒子间距离）和（四面体性）被调整，以从基础硅SW模型产生单原子mW模型。为了获得mW水，调整了、和，这些参数已被证明能最佳拟合实验数据[55]。2.2 表面模型最常用来表示金属原子之间相互作用的模型是经验衍生的嵌入原子方法（EAM），它计算纯金属或合金中原子之间的相互作用。通过将这些金属模型与升华能、晶格常数、空位形成能和溶解热的实验数据拟合，我们可以得到一个能够准确再现所需金属或合金的各种体相和表面性质的模型[61]。EAM势能遵循以下表达式：（4）其中是原子类型α的嵌入能，是原子类型β的原子电子密度，是原子i和j之间的距离，是成对势能相互作用。成对势能定义为（5）其中和是原子类型i和j的有效原子电荷。由于EAM是一个多体势能，指定每种金属的势能需要一个包含相互作用数据的输入文件。在这个实验中，选择了五种具有不同机械和热性能的金属，以尽可能多地获得它们与冰的相互作用差异，从而测试该框架的适用性和实验相关性，即较软、导电性更好的金属（银[62]和金[63]）、中等强度和导电性的金属（铝[64]和铜[65]），以及较硬、导电性较差的金属（镍[66]），其参考资料附有模拟EAM数据的来源[67,68]。这些金属的机械和热性能分别列在表1和表2中。表1. 测试金属的机械性能[69,70,71,72]。表2. 测试金属的热性能[73,74,75]。这些金属虽然在其纯形式中不常用，但往往是结构合金[76,77]或防冰涂层的组成部分[78,79,80]。在这项工作中，金属的性能通过它们在表面上阻止冰滴融化的程度来衡量，这也可以定义为抑制冰的沉积。这是一个常用的量化防冰性的指标[81]，因为融化的冰在低温下会重新冻结，并且与表面的粘附力更强[82]。因此，更防冰的金属是能够最小化冰撞击后融化的金属，而不太防冰的金属则倾向于引发更多的融化。这些金属的共同点在于它们的晶格结构（面心立方FCC）、表面拓扑（平坦）以及存在的缺陷数量（随机分配的1%）。通过提供在不同条件下的可复现金属模型，我们分离出了基底的化学性质以及每种金属原子与水分子的相互作用作为金属类型之间变化的主要来源。2.3 表面-冰相互作用使用Lennard-Jones（LJ）势能与Lorentz-Berthelot（LB）混合规则来表征mW水分子与金属原子之间的相互作用。LB混合规则之前已在类似的冰-固态系统中使用过[83,84]。mW-金属相互作用在文献中可以近似为氧-金属相互作用[85]：（6）（7）（8），其中rc是截止半径（10 ?）[86]。介电常数和范德华半径根据金属类型的不同而变化，具体数值列在表3中。表3. 金属与水相互作用的Lennard-Jones参数。2.4 模拟程序实验的初始阶段是系统松弛过程，其中液滴和表面都受到1 ns内的正则（NVT）系综处理，时间步长为2 fs，每1000个时间步长进行采样。这种平衡过程允许压力稳定，并使温度在设定点附近振荡，振荡幅度可接受。表面的温度阻尼因子为100 fs，而液滴的阻尼因子为50 fs，这些参数被优化以缩短平衡时间。整个系统的温度设定范围在120至250 K之间。选择这个温度范围是为了得到从极小温度效应到自然温度振荡会导致液滴整体温度接近0 °C的各种冰行为。这种方法用于测试所测试框架的一致性，并了解在低温下QLL厚度变化是否显著影响观察到的熔化和能量传递。在此阶段，液滴和表面的质心被固定，以防止任何漂移。平衡阶段是标准LAMMPS工作流程的一部分，以确保系统松弛和收敛[92]。平衡阶段的持续时间取决于压力的稳定，由于NVT系综中没有压力计，因此压力不受控制，而温度由Nosé-Hoover恒温器控制[93]。一旦系统确定，压力会自然振荡，因为分子会移动并找到最低势能点，即稳态。由于计算时间非常宝贵，因此必须定义一个系统被认为已平衡的容忍度。这是通过运行足够长时间的模拟来确定的，以观察压力振荡是否达到一个恒定的幅度并且不偏离平均值（应该是大约0 bar，因为这是一个仅包含固体的系统）。对于温度而言，振荡在整个设定点范围内保持相同的幅度。因此，唯一需要考虑的参数是压力，它应该稳定在设定点附近。从图2可以看出，平衡时间的选择考虑了这些条件。图2. 在NVT（160 K）下系统的平衡过程，压力的变化（a）和系统温度的变化（b）。（a,b）共享同一个x轴（模拟运行时间）。这个阶段结束后，系统仍然保持NVT，但采用不同的时间积分方法。激活了微正则（NVE）系综，但温度仍然由两个恒温器控制：Langevin恒温器以100 fs的阻尼因子控制表面温度，而液滴的温度每200个时间步长重新缩放到目标值，同时减去质心速度。显式重新缩放允许更好地控制液滴的温度，而不考虑其净速度。在撞击前立即解除控制，因为我们正在跟踪能量传递。释放质心，对液滴施加z方向速度，碰撞时的速度范围为10到120 m/s。所选的速度范围对应于不同大小的冰雹的终端速度[94,95]。小冰雹（直径小于5 cm）的速度范围为10到40 m/s，而较大的冰雹的终端速度可能超过70 m/s[95]。太阳能板是静止的，因此可能会遇到较低的液滴速度[96]。然而，为了考虑飞机的移动或风力涡轮机的旋转[97]，范围已扩展到120 m/s，这是由于液滴与快速移动的外部表面碰撞时的撞击速度上限。为了评估碰撞前后结构的破坏（在本工作中简化为“融化”），我们在OVITO中使用CHILL+算法[98]，该算法可以识别参与冰Ih结构的粒子，并已在之前的mW水研究中成功应用[30]。2.5 验证过程为了验证模型对冰-表面相互作用的预测，第3.4节将碰撞的形态与实验数据进行了比较，而第3.5节通过角度序参数（AOP）量化了液滴的结构无序程度，AOP由Tung等人定义[99]，通过（9）其中N是粒子数量，是中心粒子与键j和k之间的角度。由于四面体性定义为所有角度都为109.5°，其余弦值为-1/3，因此AOP接近0的粒子群体表示体相冰，而AOP较高的群体表示更无序的排列。由于AOP量化了“偏离四面体性的程度”，它使我们能够数值上追踪QLL（准液层）的厚度[100]。由于这种现象已经被广泛研究，因此可以将结果与实验数据进行比较，从而验证模型的有效性。QLL的计算在3.5.1节中有详细说明。3. 结果与讨论3.1 温度对碰撞结果的影响3.1.1 温度对熔化的作用熔化与温度有根本联系，因为熔点是水分子的动能克服整体吸引力时的温度[101]。当接近熔点时，即使温度的微小增加也会大幅提高相变速率[102]，这对碰撞动力学有深远影响。由于熔化与局部界面温度有关，碰撞过程中与表面接触的液滴部分会经历显著的温度升高，如果整体温度足够高，可能会导致熔化。初始温度越低，预期的融化百分比应该越低。这些数据在图3中展示。通过使用CHILL+算法[98]比较碰撞前后液滴处于六角冰相的比例来量化熔化。图3. 在30 m/s碰撞前（蓝色）和碰撞后（橙色）液滴的冰Ih比例，在三个系统温度（160、190和250 K）下。前两个温度的初始体相冰含量大致相同（160 K时为76.7%，190 K时为76.4%）。然而，碰撞后观察到显著差异，因为在160 K时含量降低到72.2%（液滴融化了4.5%），而在190 K时降低到69.1%（液滴融化了7.4%）。因此，在190 K时由于平均粒子动能更高，导致的熔化增加了62%，而在250 K时这一差异更加明显。在该温度下，初始体相冰含量为71.1%，碰撞后降低到52.1%，即液滴融化了19%，相比190 K增加了158%。3.1.2 温度对材料依赖的熔化行为如预期，更高的温度使液滴更接近熔点，将更多的分子从固态中推出，因为碰撞向液滴注入了更多的动能。这种效应在不同表面之间的分布不均匀。更防冰的表面会使材料熔化得较少，从而使液滴更容易被剪切掉。相反，不太防冰的材料会熔化更多的冰，导致其重新固化并增加与表面的结合能[103,104]。测试的金属及其相关的熔化情况显示在图4中。图4. 不同温度下，每种金属类型碰撞前后液滴冰Ih百分比的减少情况，撞击速度为30 m/s。在较低温度下，表面的化学类型对碰撞后观察到的熔化没有重要影响。这是因为在这个温度下整体熔化量低，导致敏感性低，金属之间的差异不大。然而，随着温度的升高，观察到的熔化量增加，因为整体熔化量增加。由于每个表面的地形、晶体结构以及缺陷的数量和性质保持不变，因此可以确定金属的化学性质是这些观察结果的唯一原因。铝、金和银是在190 K以上温度下最防冰（导致较低熔化）的金属，而铜和镍在碰撞后熔化更多的冰。3.1.3 碰撞过程中的能量传递深入研究可以追踪每种金属类型在碰撞过程中冰的动能和分子势能的变化，如图5所示。图5. 在250 K和30 m/s的碰撞过程中，冰的分子势能（a）和动能及温度（b）的变化，碰撞时刻用黑色标出。由于能量是守恒的，碰撞过程中发生的能量转移可以表示为 (10)，其中 Ekin 是冰的总动能，GravitationalPotentialEnergy 是冰的重力势能，MolecularPotentialEnergy 是冰的分子势能（间接反映了冰晶体结构的变形），而 ConversionToSurfaceKineticEnergy 则是转化为表面动能的部分。当碰撞发生时，从重力势能到动能的转换几乎是瞬时的，因为冰的分子会向表面移动。随着越来越多的液滴分子被压缩到基底上，它们的动能增加，冰的温度也随之升高，温度呈指数级上升，直到达到一个平台值。一部分重力势能并没有转化为动能，因为它会在冰体内产生塑性变形 [105]，并导致融化，从而增加了液滴的分子势能 [106]，另一部分则转移到表面并转化为表面的动能（图6）。在计算冰的平均温度时，LAMMPS 中使用的公式是 (11)，其中 (12) 和 (13)，mi 是水分子的质量，vi 是它们的速度，kB 是玻尔兹曼常数，NDOF 是自由度数，它取决于空间维度的数量（3）、水分子的数量（27,606）以及限制粒子运动的固定点的数量（0）。知道了这些信息，就可以计算出冰体达到其熔点（272.15 K）以上的温度的限度：(14)。

图6. 在250 K和30 m/s的碰撞过程中，表面动能的变化情况，碰撞时刻用黑色标记。因此，从图5中我们可以看出，在碰撞后的时间段内，五个表面中有两个表面的冰滴动能超过了这个熔化阈值，之后冰体应该开始融化。实际上，最不排斥冰的表面（Ag）和最不排斥冰的表面（Ni）之间的能量差为25 K。这意味着表面与液滴之间的热传递增加了，因为与铜和镍碰撞的液滴的平均温度高于表面温度（以及原始系统的温度250 K）。这种热传递并不像冰的能量转换那样瞬间完成，它还取决于材料的热属性以及冰与表面本体之间的温差。通过图6可以看出，驱动力和热导率决定了表面吸收液滴多余热量的速度：

如图6所示，最不排斥冰的金属铜和镍似乎比最排斥冰的金属吸收了更多的能量，这表明某些表面材料属性可能会导致较少的融化。温度在讨论热传递和融化时是一个重要因素，但研究速度可能会引发进一步讨论，以了解随着速度的增加，哪些材料属性会影响液滴的融化。

3.2. 速度对碰撞的影响
3.2.1. 作为速度函数的融化
速度测量时的温度为200 K，这提供了良好的灵敏度，同时将热效应降至最低，以便将速度作为参数单独考虑。实验中发现，较高的速度会导致冰体上出现更多的裂纹 [108]，这一观察结果需要通过建模数据来验证。然而，实际测量的是液滴中Ice Ih成分的变化，这间接反映了裂纹或融化的发生，如图7所示。

图7. 在200 K的系统温度下，液滴在三种冲击速度（30、60和120 m/s）碰撞前（蓝色）和碰撞后（橙色）的Ice Ih成分。不同速度下的初始冰体含量大致相同，均为液滴的74%，因为系统温度保持恒定在200 K。然而，我们注意到碰撞后有显著差异：在30到60 m/s之间，液滴的融化比例从5.6%增加到6.7%，增加了19.5%；在60到120 m/s之间，这一比例从6.7%增加到9.7%，增加了45.3%。

3.2.2. 速度对材料依赖的融化行为的影响
正如预期的那样，速度的增加会导致更多的冰融化，因为速度中的额外能量部分转化为动能，提高了液滴的温度并引发了融化。虽然不同金属在相同速度下的表现存在一些差异，但总体变化不大，如图8所示。

图8. 对于每种金属类型，在200 K的系统温度下，碰撞前后液滴Ice Ih百分比的减少情况与冲击速度的关系。在低速度下，最不排斥冰的金属和最不排斥冰的金属之间的差异较小，因为镍和铜在10 m/s和30 m/s时对冰滴结构的融化影响更大。然而，随着速度的增加，所有金属的融化现象都加剧了，尤其是那些在较低速度下融化较少的金属，因为冰结构在较高碰撞能量下吸收冲击的能力较弱。在120 m/s时，不同金属之间没有显著差异，尽管总体趋势是在最高速度下有更多的结构破坏。

3.2.3. 碰撞过程中的能量传递
为了比较不同金属在较高速度下的表现，观察60 m/s时的能量传递情况是有用的，在这个速度范围内金属之间的差异最为明显。
总体而言，由于系统温度较低，图9中的(a)部分显示冰的分子势能增加较少（后者最大增加200 eV，而前者为100 eV），因此融化也较少。在(b)部分，与前一节类似，与镍和铜的碰撞比与铝、银和金的碰撞更能提高冰滴的动能。这也意味着温度的升高幅度更大。实际上，与镍的碰撞使液滴温度升高了50 K，而与银的碰撞使冰的温度升高了20 K，两者相差30 K。这种温度的升高增加了液滴与表面之间的温差；在与不太排斥冰的金属碰撞时，热传递的驱动力更大，从而在足够的时间后使表面本体的动能更高，如图10所示。

图10. 在200 K和60 m/s的碰撞过程中，表面动能的变化情况，碰撞时刻用黑色标记。

3.3. 将材料属性与抗融化性能联系起来
在3.1节和3.2节中，我们讨论了温度和速度如何影响融化过程以及能量如何在五个测试表面（银、金、铝、镍和铜）之间的传递。因此，研究这些金属的机械和热属性对于理解实验结果至关重要，以便确定建模可以用来预测抗冰性能的具体属性。

3.3.1. 机械属性
在比较金属与固态冰滴的相互作用时，首先要验证它们的机械属性是否存在显著差异。由于碰撞发生在非极端速度下，因此预计金属不会出现塑性变形，正如McElligott等人的研究 [25] 中反复冲击后未观察到的现象。因此，比较与弹性变形和能量传递相关的机械属性是有意义的，例如杨氏模量和硬度 [109]，如第2.2节中的表格1所示。
从表格1可以看出，最不排斥冰的金属（金、银和铝）的共同特点是它们的杨氏模量和维氏硬度都低于其他两种金属。为了进行更显著的比较，应该对表格1中显示的平均融化百分比与每种机械属性进行线性回归。为了将机械属性与热属性区分开来，融化百分比是在相同温度（200 K）下，不同速度下金属之间的比较所得到的。

从表格4可以看出，平均融化百分比与机械属性之间的相关性相当高，但由于三种测试金属在融化百分比方面存在相似性（Al、Ag和Au），因此没有达到统计显著性所需的阈值（）。尽管如此，这些结果仍然令人鼓舞，因为文献支持以下观点：硬度较低和杨氏模量较低的金属（1）通过弹性变形吸收更多的冲击能量，从而减少了冰的势能转化为动能的比例，降低了观察到的融化量 [110,111]；（2）产生的接触压力较小，减少了摩擦加热和界面温度 [112,113]；（3）减少了界面处的应力集中，减少了液滴的变形 [114,115]。

3.3.2. 热属性
在比较金属与其与冰的相互作用时，第二个步骤是研究它们如何有效地在其内部传递和分布热量。为了测试热属性对冰滴融化的影响，必须在不同的温度下进行比较，并考虑驱动力的变化。这可以隔离出感兴趣的变量，同时控制所有其他条件。比较的热属性包括导热性、比热和热扩散率。导热性和热扩散率决定了热量从撞击区域重新分配到整个金属本体的速度，控制液滴经历的温度升高 [116,117]。比热决定了金属在撞击过程中的温度升高程度，从而影响热传递的驱动力，并决定了液滴会经历的热量 [118]。第2.2节中的表格2比较了测试金属的热属性。

与表格4不同，表格2并没有明确表明热属性是否显著影响冰的融化。因此，需要进行统计分析（表格5）：

表格5. 熔化百分比与热属性的相关性。比热基本上没有显示出对观察到的冰融化的预测能力。这可能是由于交换的能量较少，所以金属的温度升高不足以影响结果。热导率与融化百分比之间的相关性更为显著，但计算出的p值（0.323）仍远高于传统的0.05阈值，这种差异可能无法仅通过样本大小来解决。最后，热扩散率的R2值（0.690）高于导热性和比热，相关的p值更接近0.05阈值。这表明融化与热扩散率之间存在弱相关性，但这种相关性不能轻易地与其他变量（如基底的机械属性）区分开。理解冰滴与金属基底之间的相互作用需要一个综合考虑这些相互作用的框架，而不仅仅依赖于单一属性来预测融化。

3.3.3. 热属性和机械属性的结合
从3.3.1节和3.3.2节中，我们确定了两种似乎与碰撞后表面冰融化量有关的属性是杨氏模量和热扩散率。由于这两者结合对本研究极为重要，因此创建一个将两者联系起来的Ashby图是有用的，如图11所示。

图11. 本研究中研究的金属（蓝色；原始数据来自其他来源 [69,70]）以及尚未研究的额外金属（橙色；原始数据来自其他来源 [119,120,121,122,123,124]）的热扩散率与杨氏模量的关系；虚线对应于优值设计指南的指标。正如我们所知，这项研究并未证明这些属性与冰融化减少之间存在明确的相关性。然而，结合文献数据，可以推断出这些性质确实会影响材料的抗冰性能。杨氏模量似乎与抗冰性能成反比，而热扩散率则成正比。因此，我们创建了一个评估指标，用一个变量来排名材料的理论抗冰能力：(15) 其中 是抗冰性能指数， 是热扩散率， 是杨氏模量。图11中不同材料的评估指数比较结果可以在图12中看到。图12显示了图11中所有金属的评估指数（蓝色表示已研究，橙色表示未研究）。根据图12，银在防止融化方面表现最佳，其评估指数为2.30 mm2/(s·GPa)。相比之下，310不锈钢的评估指数仅为0.02 mm2/(s·GPa)，比银低115倍以上。在测试的材料中，最不具备抗冰性能的是镍，其评估指数为0.12 mm2/(s·GPa)，比银低19倍。由于该指数中的每个参数都与性能有密切关联，因此统计分析该指数与融化百分比之间的关系是值得的。这一分析在表6中进行了说明。表6显示了研究材料的评估指数与融化百分比之间的相关性。尽管单独看热扩散率和杨氏模量都与碰撞后的融化没有显著相关性，但两者的比率显示出很强的相关性，其R2值为0.905，p值低于0.05的阈值（0.013）。这表明融化响应是由这些性质的比率控制的，而不是由任何一个性质单独决定的，但需要通过更广泛的材料集合进一步验证。尽管如此，这个指数提供了一个用于根据抗冰能力对金属进行排名的变量，但它只是决定选择哪种涂层材料的众多因素之一。例如，在航空航天应用中，通常还会考虑密度、强度、价格和耐久性等因素[125]。本文提供了另一个变量，用于在比较具有相似机械性能的金属时进行验证。虽然这个系统使用的模型纯粹是理论上的，但它可以通过类似工作的实验数据进行验证，这增加了其可信度。

3.4 实验验证
文献中存在研究固体冰与表面碰撞的实验工作[25,126]；这些实验通常涉及微米级别的液滴，碰撞时间可以在毫秒级别测量。由于我们在尺度上受到限制，因此使用分子动力学很难实现这一点。因此，将纳米尺度下的系统和纳秒级别的滴落时间与宏观尺度的碰撞联系起来是一个挑战。然而，Liu等人（2022年）[127]和Wang等人（2025年）[30]的研究表明，在实验工作和分子动力学模拟中，液滴碰撞在形态上有相似之处。同样，从关于固体碰撞的研究中，我们也能发现McElligott等人[25]的研究与我们的分子动力学模拟之间存在形态上的相似性，如图13所示。这些相似性不仅体现在液滴碰撞过程中的行为（没有反弹），也体现在碰撞后的液滴形状上。

3.5 准液层（QLL）的结构和作用
3.5.1 准液层的起源和性质
在建模结果中，我们观察到液滴外部有一层预熔化的层（图14），这可以被识别为引言中提到的准液层（QLL）。QLL在温度低于材料熔点时形成，因为形成该层的分子与材料内部的六方结构之间的键合较少[128]。QLL成为了材料内部有序状态与外部无序状态之间的过渡状态，在这种状态下，氢键的缺失使水分子进入能量较低的状态[37]。这给粒子提供了更多的自由度，从而减少了对热效应的抵抗[32]。这可能起到了缓冲作用，使能量在液滴内部更加均匀地传递。这种QLL的存在可以解释为什么本研究测试的温度范围内，半径为60 ?的纳米液滴在碰撞后永远不会反弹。

3.5.2 QLL厚度：分子动力学与实验的对比
通过TIP4P/Ice建模，Conde等人[36]发现，在-30°C到-10°C的温度范围内，QLL的厚度在0.3到0.8纳米之间。Michaelides等人[32]也证实了这些值。由于QLL占液滴半径的10–15%，因此它足以显著减弱碰撞冲击。相比之下，McElligott等人[25]实验中的液滴半径为0.7毫米，这意味着计算出的QLL仅占其半径的0.0002%，极大地影响了其缓冲能力。在分子动力学模拟中增加液滴直径会降低QLL相对于半径的比例，从而减少碰撞的缓冲效果，导致整体结构变形加剧。然而，由于液滴体积随其半径的三次方增加，使用当前的计算工具来降低QLL与半径的比率会带来过高的计算成本。不过，由于实验的控制系统较为简单，液滴的实际温度可能更接近其熔点，从而导致QLL的值较大（高达100纳米[32]，即其半径的0.02%）。不同实验之间QLL厚度的变化可以解释为什么有时液滴会反弹，而有时会粘附在表面上。在当前模型中，如果温度保持不变，所有模拟运行中的QLL厚度都是相同的。这种可重复性是建模的一个主要优势，但它无法自然再现实验中存在的变异性。因此，模型中固体液滴的碰撞机制只能再现液滴在碰撞后粘附在表面上的情况，而对于其他情况，则假设液滴具有更厚的QLL（温度高于-5°C）。

使用分子动力学模拟，研究了固体液滴撞击金属表面的动力学。为了使液滴达到热力学最优状态，需要在恒定温度和体积下进行平衡，让QLL有足够的时间在厚度和形态上稳定下来。回想一下，计算QLL的方法在第2.5节中已经有介绍，即使用AOP来量化液滴中心不同距离处的冰相。当我们从液滴中心向外移动时，会发现AOP突然升高的点，这标志着QLL的开始（因为界面冰的AOP值高于内部冰）。

图15显示了平衡完成后，QLL厚度随材料温度变化的演变情况：(a) 不同温度下AOP随距离液滴中心的距离的变化。(b) 不同温度下球体内部的平均QLL厚度。如图15a所示，温度的升高导致AOP提前增加，表明QLL厚度增加，从而增强了碰撞过程中的缓冲效果。此外，图15b中的数值与Michaelides等人[32]之前的实验和建模结果一致。这是该分子动力学系统的一个固有特性，只有通过更大的系统才能克服这一限制，但这需要更多的计算资源。不过，这种框架在比较不同材料时仍然有用，因为分子动力学中与纳米液滴相关的限制在各种基底之间是相同的。因此，可以使用任何一种材料的实验结果与建模结果结合起来，预测未来的实验数据，这是一个有用的功能。

4. 结论
本研究提出了一个可复制的分子动力学框架，用于评估不同材料表面的动态抗冰能力，通过分析机械和热性能如何影响碰撞后的冰融化。我们使用LAMMPS软件和OVITO可视化工具，在120到250 K的温度范围内，10到120 m/s的撞击速度范围内，以及五种不同的金属表面（金、银、铜、铝和镍）下，结合mW水模型、EAM/Alloy模型和LJ势能，研究了固-固-冰-表面碰撞。在NVT平衡后，跟踪了冰滴的势能和动能、温度，以及表面的动能和温度在碰撞前、碰撞中和碰撞后的变化，以确定能量传递和融化百分比。然后利用这些数据进行了统计分析，以选择影响融化程度的表面材料属性。最后，通过比较实验和模拟中的撞击形态和QLL形成情况，使用AOP随距离液滴中心的径向距离的变化作为QLL的主要指标，进行了实验验证。我们的主要发现是，具有较高热扩散率和较低杨氏模量的金属倾向于减少碰撞后的冰融化，这一点通过引入一个评估指数得到了支持，该指数将这两个参数联系起来。我们发现这个指数与融化比例之间存在统计上的显著关系，其中银是研究中最有效的材料。这项研究证实了先前的实验研究，这些研究将碰撞能量传递与表面的热和机械性能联系起来，但更重要的是，它结合并定义了一个新的参数——一个特定于应用的指数，用于测量表面的抗冰性能。这项研究填补了在将表面的整体性质与其抑制多个温度和速度下融化能力联系起来方面的知识空白，并限制了可能影响结论的附加因素。

本研究中的模拟使用了分子动力学（MD），但它存在一些实现上的局限性：基于力场的假设（例如使用mW粗粒化的水模型，忽略了显式的氢键作用，以及水分子与表面金属之间的Lennard-Jones相互作用不完美），模拟时间可能不足以使系统完全平衡，以及非平衡MD模拟中常见的采样误差。通过仔细的模拟设计，试图减少这些局限性。此外，该系统忽略了表面粗糙度、氧化层和系统中的杂质，这限制了与真实世界表面的直接比较。最后，由于该指数基于组成相似（纯金属）和性质相似的金属，可能存在过度拟合的问题。未来的工作应该通过引入不同的表面、更真实的表面组成和地形，并对本文开发的评估指数进行实验验证，来进一步验证该系统。