从纳米到宏观：粗糙表面接触面积随正载荷的演化

摘要

粗糙表面真实接触面积随正载荷的演化具有重要基础和实际意义，涉及地震动力学到机器磨损。本文弥合了原子尺度到宏观尺度正常接触行为的空白。真实接触面积由大量离散接触簇形成，远小于表观表面积。离散接触簇的分布及其相互作用是揭示接触固体机制的关键。为此，采用格林函数分子动力学（GFMD）研究接触簇如何从原子尺度演化至宏观尺度及簇间相互作用。发现簇间相互作用对其形成有强烈影响。接触簇的形成和分布比凸峰模型预测复杂得多。忽略簇间相互作用会导致高估接触力。真实接触中，因凸峰间相互作用，接触簇更小、更离散。理解接触面积随正载荷的精确本质对后续摩擦研究至关重要。

背景

大多数宏观表面被认为粗糙且分形。粗糙表面间接触行为比完美光滑表面复杂得多。真实接触面积由大量离散接触区域（簇）形成，远小于表观表面积。正载荷与接触簇大小、形状、分布是揭示接触行为的关键，对后续摩擦研究必不可少。

自1960年代以来为获得接触面积-载荷关系提出了众多模型。其中凸峰模型最简单流行。Greenwood和 Williamson早期应用假设凸峰同半径不同高度描述接口粗糙度。此后凸峰模型盛行，大量文献涌现。Whitehouse和Archard引入凸峰尖端随机曲率半径发展G-W模型。Nayak将随机过程理论引入高斯粗糙度分析，后被Bush等人用于粗糙表面接触。

凸峰模型基本假设之一是忽略凸峰间相互作用，即潜在接触凸峰可预先由表面几何确定。但这可能导致接触力和面积估算不准。为获得接触簇演化及其相互作用，我们利用格林函数分子动力学（GFMD）研究分形粗糙表面。

本文旨在弥合原子到宏观尺度正常接触行为的空白。通过考虑凸峰相互作用的数值示例展示接触面积从原子到宏观尺度的演化。后续讨论先简述分形表面生成、GFMD模型、接触簇识别算法和数值实验设计，再聚焦接触簇形成发展及其对界面行为的影响。

方法

粗糙分形表面生成

研究粗糙表面接触行为需为数值模型生成表面。分形表面生成算法多种，本文采用傅里叶变换法生成分形粗糙表面（图1）。需四个参数确定分形粗糙表面几何：最大频率（w_H）、最小频率（w_L）、赫斯特指数（H）、振幅标准差（P）。表面基本统计参数如RMS高度（√M₀）、RMS斜率（√M₂）、RMS曲率（√M₄）是界面行为关键参数，M_i为表面谱第i径向矩。值得注意的是表面统计参数M_i与轮廓统计参数m_i关系：M₀=m₀, M₂=2m₂, M₄=4/3m₄。凸峰密度n（表面峰/谷）由公式（1）确定：n = 1/(6π√3) (m₄/m₂)。

表面总峰/谷数N由公式（2）表达：N = A₀×n = A₀/(6π√3) (m₄/m₂)，A₀为表观表面积。自相似分形表面统计参数与输入参数（w, H, P）通过公式（3）关联：m_i = ∫_{w_L}^{w_H} ω^i Φ_φ(ω) dω = ∫_{w_L}^{w_H} ω^i B ω^{-(1+2H)} dω，B为表面粗糙度常数，与P相关。公式（3）和（2）表明分形表面峰/谷数依赖波长和赫斯特指数。分形表面统计性质详细讨论见文献。

傅里叶变换算法中，典型设H=0.5，最大频率w_L=1/(24σ)，最小频率w_H=1/(256σ)，振幅标准差P=0.69，系统尺寸512×512原子（初始间距1.12σ）。这些参数生成表面统计参数：RMS斜率√M₂=0.077，RMS曲率√M₄=0.0077。公式（2）得表面总峰/谷数150，数值计数峰数158、谷数159，误差5%以内，系统尺寸统计上可接受。实际上增至2048×2048原子统计参数与小系统一致。

GFMD模型

粒子间相互作用实验难捕捉。近年分子动力学模拟粒子间相互作用以探究接触/摩擦分子起源，但大规模计算开销大。引入GFMD模拟表面因其高效。GFMD用分子动力学模拟界面原子相互作用（此处两层），非界面层通常表现弹性行为，用格林函数模拟。从而将大原子系统简化为界面两层原子（图1），大幅降低计算开销。GFMD详细讨论见文献。GFMD模型中用Lennard-Jones（LJ）势模拟粒子间相互作用：u(r)=4ε[(σ/r)¹²-(σ/r)⁶]，ε为势阱深度，σ为粒子间势为零的有限距离，r为粒子间距。取ε、σ、ε/σ为能量、距离、力单位。按LJ势，r=2¹/⁶σ≈1.12σ时粒子间力为零；r>1.12σ吸引，r<1.12σ排斥。本文不考虑粘附，截断距离设1.12σ。原子层晶体结构为面心立方（FCC）。利用对称性仅取界面层形成表面几何（图1），下方弹性块用格林函数模拟。

接触簇识别方法

界面观察到三个尺度（图1）：(1)原子尺度，LJ势模拟；(2)凸峰尺度，接触原子群效应；(3)宏观尺度，接触簇群效应。接触簇大小、形状、位置、分布是连接分子行为与界面性质的关键桥梁。纳米尺度原子接触区难定义。我们定义接触原子为法向分力f_z>0。随后连通接触原子定义为接触簇。使用标记技术搜索接触簇，采用避免递归搜索的改进算法加速。算法流程图见图2，从GFMD模拟提取原子力数据。算法分八关键步骤：

步骤1. 开始行搜索获取新原子数据，逐行搜索原子。
步骤2. 判断原子是否接触。非接触回步骤1，接触进下一步。
步骤3. 比较当前原子与同行前一原子。若前一原子也接触，合并当前原子至前一原子所属簇，标记相同编号。若前一原子非接触，标记新编号（前一编号+1）。
步骤4. 判断是否最后原子；否回步骤1，是进列搜索。
步骤5. 开始列搜索获取新原子数据，逐列搜索。
步骤6. 判断原子是否接触。非接触回步骤5，接触进下一步。
步骤7. 比较当前原子与同列前一原子。若前一原子也接触且属不同簇，合并当前簇至前一原子所属簇，标记相同编号并存储。若前一原子非接触或同簇，进下一步。
步骤8. 判断当前原子是否最后；否回步骤5，搜索结束。

数值实验设计

双粗糙表面接触问题可通过引入等效弹性模量E*简化为一复合刚性粗糙表面与平坦弹性表面问题：1/E* = (1-ν₁²)/E₁ + (1-ν₂²)/E₂，E₁、E₂为上下表面弹性模量。为简化，考虑刚性粗糙表面接触弹性光滑表面，研究接触簇形成发展及其力-面积行为。以下讨论使用上述生成表面（上表面刚性粗糙E₁=∞，下表面光滑弹性E₂=3ε/σ³），ν₁=ν₂=0.5。系统尺寸512×512原子（初始间距1.12σ），x-y平面周期性边界条件。弹性块深度1024原子层（初始间距1.12σ）。常规分子动力学系统将含2.68亿原子；GFMD模型减至524,288原子（两层），见图1。逐步将粗糙表面（上）推入平坦弹性表面。刚性表面载荷由位移控制。每步位移载荷0.01σ，GFMD算法更新原子位置直到原子力满足收敛准则L¹-范数=0.01ε/σ。最大迭代数50,000避免无限循环。

结果与讨论

接触簇分布与发展

凸峰模型认为凸峰为球形或椭球形且忽略接触凸峰间相互作用。本文凸峰模型凸峰从上述生成表面提取。凸峰模型中潜在接触凸峰可预先由表面几何根据高度确定；即表面峰/谷按高度形成接触簇。然而现实中凸峰形状不规则，通常多个相邻凸峰随接触力增加合并为大簇（图3）。初始观察到六个独立凸峰，随接触力增加最终合并为大接触簇。这表明假设凸峰间距足够大互不影响可能导致不准确结果。

图4显示随接触面积增加，簇数先增后减，而表面凸峰数随接触面积单调增加。这是图3解释的合并效应所致。

簇合并效应在凸峰模型和GFMD模型均观察到。但在相同接触面积下，GFMD模型簇数远多于凸峰模型（图5）。GFMD模型簇数几乎为凸峰模型两倍。主因是凸峰模型忽略凸峰间相互作用。GFMD模型中接触簇相互影响，簇产生的位移场在整个表面连续，大量簇集合位移导致弹性表面新几何，影响新簇形成。因此簇形成不仅基于刚性粗糙表面高度，也受光滑弹性表面变形影响。图6展示两模型不同接触面积下簇分布：5%接触面积时，凸峰模型簇数17，GFMD模型34；10%时分别为24和52。表明GFMD模型簇更离散。GFMD模型平均簇尺寸更小，但大多簇与峰/谷重合（图7）。凸峰模型视峰或谷为潜在凸峰（取决哪面接触）。但图8显示随接触面积增加，峰和谷均可形成接触簇。小接触面积时多数接触簇为表面谷；大于10%表面积时越来越多峰也形成簇。

接触面积-载荷关系

正载荷下力-面积关系是接触行为核心。早期模型常假设凸峰为球形椭球形。但真实接触簇复杂得多。本节比较三模型力-面积关系：(1) GFMD模型；(2) 凸峰模型（AM），凸峰直接从表面提取再用GFMD压入平坦弹性表面（确保接触无意外簇形成）；(3) Greenwood-Williamson模型（G-W），凸峰转等效球。球半径由公式（6）得：1/R = 8/3 (m₄/π)¹/²。GFMD模型和提取凸峰的凸峰模型，界面总力由GFMD提取各簇力求和。G-W模型用Hertz理论计算各凸峰力（同GFMD材料属性），总力F = Σ(4/3 E* R¹/² (d-z_i)³/²)，z_i为凸峰高度，d为刚性表面施加位移，f_i为Hertz接触理论凸峰接触力。

图9比较三模型力-面积关系，均呈线性。观察到GFMD总力远小于凸峰模型和G-W模型。凸峰模型F为GFMD预测1.80倍，G-W模型为1.54倍。这可由接触簇RMS斜率解释。正载荷与RMS斜率成正比，即L∝√M₂。GFMD模型接触面积由更多簇组成，穿透深度比凸峰模型浅。因凸峰尖端斜率较小，GFMD模型接触簇RMS斜率也较小。图10显示三模型接触簇RMS斜率随接触面积变化。可见GFMD模型接触簇RMS斜率小于表面RMS斜率0.077，而其他两模型簇RMS斜率大于表面RMS斜率。

结论

为探明接触面积从原子到宏观尺度演化，利用GFMD模型研究分形粗糙表面接触问题。以法向力>0定义原子接触，在同一系统研究三个长度尺度。发现凸峰间相互作用对接触簇形成至关重要。部分簇足够接近可合并为大簇。真实接触区因弹性光滑表面变形比表面几何预测复杂得多。大多数接触簇位置与表面峰/谷重合，但簇尺寸更小，形成不由表面凸峰高度决定。随接触面积增加，峰和谷均可形成簇。GFMD模型力远小于凸峰模型，但簇数远多。GFMD模型接触簇RMS斜率小于凸峰模型，解释为何凸峰模型导致更高压力。研究表明真实接触面积不能简单由表面几何预测。真实接触面积随正载荷的演化对后续摩擦研究至关重要。