由水合离子产生的水合斥力与双电层力、van der Waals力一起被统称为表面力，可以用Derjaguin-Landau-Verwey-Overbeek(DLVO)理论描述[12]。表面力的明显作用范围在纳米量级，是一种短程力[13-14]，可以提供高压来支撑外部负载[15]。虽然水溶液中表面力对改性聚合物材料润滑性能的明显改善已经在多项实验中得到验证[16-17]，但在摩擦过程中两表面间的表面力分布情况和膜厚情况还没有得到直接观察确认，只能通过数值计算结果来尝试阐明。一些针对高硬度材料的数值计算研究表明，水合力可以显著增大摩擦副间的膜厚，减小膜内的速度梯度，降低直接接触比例，从而降低摩擦[18]。但在低弹性模量条件下，考虑表面力的接触情况还处于未知状态。因此，有必要建立考虑表面力的低弹性模量聚合物材料接触计算模型。
1 考虑表面力的接触模型
1.1 接触模型中使用的公式和方程
本文的算法模型针对的是球-盘模型中聚合物球与蓝宝石盘在盐溶液中的接触过程，如图 1所示。
该模型中适用的表面间距离表达式如下：
其中：hr表示刚体位移; $R-\sqrt{R^2-x^2-y^2}$ 表示理想球体表面与平面的间隙，R为上试样球的半径；v表示表面上分布的压力造成的弹性形变；δ表示表面的粗糙度；(x，y)是在图 1所示直角坐标系下，接触区中某压力作用点的位置，原点o位于接触平面上和上试样球心的正下方。
弹性形变使用Boussinesq积分来进行计算：
其中：Ω是整个接触算法的计算域；ξ和η也是接触区中另一个压力作用点的位置坐标；p是计算域中的压力分布，E′是摩擦副材料的等效弹性模量。
其中：E1、E2分别为上、下试样的弹性模量，ν1、ν2分别为上、下试样的Poisson比。由于式(2)只能应用于小变形假设成立的条件下，因此本文讨论的聚合物材料弹性模量不低于0.2 GPa，以防止变形过大。
1.2 表面力的描述方程
根据DLVO理论，表面力为
其中：pvdW是van der Waals力，pdl是双电层力，phf是水合斥力。
其中：AH为摩擦副的Hamaker常数，是两摩擦副材料的Hamaker常数的几何平均值；h是节点上两表面间的距离。
其中：ε和ε0分别是真空中的介电常数和溶液的相对介电常数，本文选取的溶液是质量分数为0.9%的NaCl水溶液；κ是Debye长度的倒数；ψ1和ψ2分别是摩擦副上、下表面的表面电势。
其中：ΔGh是阳离子的水合能，对于本文中的Na+离子取410 kJ/mol；q是元电荷；σ1和σ2是上、下表面的表面电荷密度。σ1和σ2可以通过Graham方程由表面电势计算得出：
其中：kB是Boltzmann常量；T是热力学温度；c是溶液中远离表面处的钠离子物质的量浓度。分别代入ψ1和ψ2即可得到σ1和σ2。表面电势可以由实验测定。
由式(4)—(9)可以计算出上、下表面分别为UHMWPE、蓝宝石表面时，在质量分数0.9%的NaCl水溶液中的表面力与表面间距的关系。代入的物理量数值为：上、下表面的Hamaker常数分别为AH1=0.6×10－20 J、AH2=5.3×10－20 J，两者的几何平均为AH=4.1×10－20J；ε=78.4；T= 298 K；Debye长度κ －1=0.77 nm；ψ1=20 mV，ψ2=30 mV。经式(8)—(9)计算，p0=799 MPa。计算结果如图 2所示。
计算时，对λ0<0.2 nm的部分，表面力值取为常数psf, max。这样取值可以防止在计算中出现某个节点上的间距收敛到h为[0, 0.2) nm的情况，即可以在计算结果中规避不足一个水分子大小的两表面间距，使计算结果符合物理规律限制。计算中，模型取式(4)的计算值作为该表面节点与对偶面节点之间的所有表面作用的大小的精确值。
1.3 接触计算方法
将计算域对式(2)进行离散和无量纲化可以得到
其中：Vi, j是节点(i, j)上的无量纲化形变，i和j分别是节点沿x和y方向的编号，节点在接触区内沿x和y方向等距分布，节点编号沿轴的正向从1递增；P k, l是节点(k, l)上的无量纲化压力；C i－k, j－l是离散后的节点(k, l)上的无量纲化压力对节点(i, j)上的无量纲化形变的影响系数。张量 C的具体计算方法可见文[24]；变形矩阵 v和压力矩阵 p的无量纲化方式为：
其中：H为无量纲化标准长度，ph是球形接触的Hertz接触压力。
其中：b是Hertz接触半径，w是外加接触载荷。
计算P与 C的卷积时，使用了快速Fourier变换(FFT)加速计算[25]。
在计算中假设下表面是一个不动的平面，接触过程是球形的上表面向下表面接近，直至两表面间作用使得接触计算区域内各点的压力和形变均满足式(2)的压力-形变关系的过程。其中，如果下降过程使得2个表面直接接触，则采用文[21]的处理方式，迭代的关系为
其中：Pi, jold和 P i, jnew分别是迭代中节点(i, j)上本次和上次的压力；H i, jold是迭代中节点(i, j)上上次的膜厚；C 0, 0表示离散后的节点上的无量纲化压力对自身无量纲变形的影响系数。
此迭代过程的物理含义是：两表面间的某节点上如果存在互相穿透，则认为此节点处于接触压力的作用范围，迭代将令该点的压力增加，直至两表面间的穿透程度足够小或无穿透。此外，若两表面不直接接触，而是进入了表面力的作用范围，迭代则通过表面力和表面间距的关系式(4)和(7)进行控制。结合之前的研究结果，表面力的作用范围认定为2 nm以内。这样的简化忽略了表面力中双电层力在超出2 nm范围时的长程作用，这是因为双电层力的长程作用比2 nm以内的van der Waals力和水合斥力弱许多，在接触中不能明显体现出来。因此某节点上两表面间的作用力情况如图 3所示。
由此，本文建立了图 4的数值计算方法。整个数值计算过程分为2个循环：循环1是压力迭代的循环，作用是在固定的刚体位移下，使得压力与形变的关系既满足表面力公式，又满足两表面不得互相穿透的边界约束；循环2是载荷平衡的循环，作用是使得表面力和边界反力的总值与载荷值相平衡。输入数据后，首先计算C中的各个元素并保存，初始化表面形貌并将压力矩阵置为0；然后将两表面之间的最小间距调整为无作用力的最小值2 nm；最后再逐渐缩小两表面的间隔距离。两表面间隔距离的缩小方法如下。
定义载荷误差为
其中px, y是本次迭代获得的各节点(x, y)上的压力的带量纲数值。这里的节点坐标是节点对应的实际位置坐标值。
定义压力迭代变化量为
其中px, yold是上次迭代获得的各节点(x, y)上的压力的带量纲数值。
每当压力迭代变化量满足εp<ε1时，判断压力迭代收敛。ε1是一个设定的小量，用以规定最终结果中压力迭代收敛的程度。压力迭代收敛后，将两表面之间的距离缩小量规定为
其中α为调整系数，用于建立间距调整值与载荷误差之间的关系，可以通过它来调节收敛的快慢和稳定性。α越大，程序收敛越快，但越可能进入不稳定状态；α越小，程序收敛越慢，但状态越平稳。
两表面逐渐接近到一定程度时，如果有εw<ε2，则认为载荷平衡也已经收敛，程序结束。ε2也是一个设定的小量，用以规定最终结果中载荷平衡迭代收敛的程度。
ε1、ε2和α均为控制程序迭代速度与判断收敛的无量纲小量，是调节程序迭代的重要参数。
2 计算条件和结果
由于低弹性模量材料表面形貌在接触区尺度上不会有非常密集的峰谷结构，计算中使用了64×64个节点的稀疏网格以提高计算效率。以接触区正中心为原点建立图 1的平面直角坐标系o-xy，来标记计算区域中各节点的位置。沿着x和y轴的正方向，节点序号从1到64递增。后文用到的y=0线上的数据结果，实际上就是y方向上节点序号为32的所有点的数据集合。迭代参数分别设置为ε1=10－5和ε2=10－4，α=10－6。
2.1 表面的构建
2.2 表面力对接触区域状态的影响
选用文[9]使用的UHMWPE的数据代入本算法进行计算，即：E1=0.7 GPa，ν1=0.46，R=6.35 mm；E2=410 GPa，ν2=0.3；w=3 N。此条件下得到：b=200.53 μm，pH=35.622 MPa。
为了探究两表面的完全分离与表面力的关系，在计算中设置表面力分别为原表面力的0、0.01、0.02、0.05、0.1、0.2、0.5、1倍，如图 7所示。
计算的结果如图 8所示。从图 8a可以看出，在相同载荷下，随着表面力幅度的增加，接触两表面的距离逐渐增加，这表明水合层的法向承载能力不断增强。其中，0.05倍原表面力幅值线出现了明显的一部分被水合层承载、一部分不能被水合层承载的情况。图 8b是表面力取0.02、0.05和0.1倍原表面力幅值时，接触区域y=0线上的压力分布与对应的表面力最大值的对比。从图 8b可以看出，在最大值线以上的部分，就是两表面直接接触的部分。图 8c和8d分别是0.02倍和0.05倍原表面力幅值下直接接触的部分形状，0.1倍原表面力数值条件下没有直接接触区域，故没有画出。由图 8c和8d可以看出，随着表面力幅值的增加，表面之间接触的部分越来越少。从图 8e和8f可以看出，随着表面力幅值的增大，直接接触承载比例逐渐下降，表面力承载相应上升，最小表面间距在变为正值后，随表面力幅值的上升逐渐上升。由此可以确定，摩擦副两表面的完全分离是由于表面力的幅值足够大造成的。表面力最大值小于实际接触压力最大值和两表面完全分离之间是等价的关系。
2.3 弹性模量对接触区域状态的影响
2.4 弹性模量和表面形貌对表面力接触的联合影响
本节定性地讨论弹性模量和表面形貌对接触的联合影响，不考虑接触区的尺度，以一个统一表面方程描述上试样表面。即讨论接触区内有相同的粗糙峰数量和分布时，表面力对不同弹性模量上试样接触的影响情况。
使用图 5描述的表面进行接触计算，接触区域y=0线上的接触情况如图 10所示。图 10a中的纵坐标范围为－5~700 nm，而图 10b则为图 10a中纵坐标0~3 nm部分的放大。可以看到，在同样的表面力幅值取原值0.5倍的条件下，低弹性模量摩擦副的表面被明显分离，而高弹性模量材料则没有被明显分离。图 10a显示，在弹性模量较高时，表面的粗糙结构没有被压力消除，导致接触区内真实的接触区域呈现点状分散。图 10b则表明，在0.5倍原表面力幅值的表面力作用下，接触区域y=0线上，高弹性模量材料在直接接触区域直接穿透了水合层，而低弹性模量材料则完全被水合层分离。直接接触承载比例的变化也说明了这一结论，在计算结果中，上试样弹性模量为0.7 GPa时，增添0.5倍原表面力幅值的表面力，即可使直接接触承载比例从100%下降到0%，表面完全分离；而上试样弹性模量为70.0 GPa时，增添0.5倍原幅值的表面力，仅能使直接接触承载比例从100%下降到83.93%。可以看出，在表面有相似的粗糙度时，水合力作用于高弹性模量材料的效果明显差于它作用于低弹性模量材料时的。因此，在设计水合效应减摩的改性聚合物材料时，在满足承载要求的前提下，建议选择弹性模量较低的材料。
3 结果讨论
3.1 聚合物中的塑性形变分析
对于表 1给出的数值计算条件，计算结果显示，沿压力方向的最大变形量达到12.94 μm。由于下试样弹性模量较大，因此变形主要发生在上试样。比较大的形变蕴含着塑性形变的可能，故需要探讨上试样体内的塑性形变情况。在上试样内部建立一个在x、y、z方向上64×64×32个节点的立方网格划分，使用文[25]中的由表面压力分布计算应力的公式：
其中：σα, β, γij表示上试样内部节点(α, β, γ)处的应力，而i, j=x, y, z，两者遍历所有x、y、z可能的组合，表示应力的不同分量，如σα, β, γxx表示该点沿x方向的正应力，σα, β, γxy表示该点x方向上沿y方向的剪应力，等等；pζ, ξ表示接触计算域内平面节点(ζ, ξ)上的压力值；Dα－ζ, β－ξ, γij是pζ, ξ对σα, β, γij的影响系数。使用von Mises屈服准则即式(21)是否成立判断上试样内部是否发生了塑性形变。
其中：σi和τij为各个方向的应力分量；σs是材料的屈服极限，本文中对UHMWPE取23 MPa。经过计算，在2.2节的计算条件下，上试样内部的von Mises应力最大值为σmax=19.6 MPa，满足σmax<σs。
因此可以确定，在整个试样内部，均没有发生塑性形变。如此大的形变中还没有发生塑性形变，主要是因为低弹性模量导致整个接触区域中的接触压力均较小。另一方面，从应变的角度来讲，在屈服强度固定的前提下，低弹性模量导致了较大的屈服伸长率。
3.2 节点间细节对结果的影响
这两个表面在表面力的作用下在y=0线上的、在表面力作用下的接触结果展示如图 12所示。图 12a和12c中，“省略表面”指的是未叠加正弦波动的表面的计算结果，“叠加表面”则是叠加后的结果；图 12b是用叠加正弦波动后的表面间距分布减去未叠加正弦波动的表面间距分布得到的；图 12d同样是用叠加正弦波动后的压力分布减去未叠加正弦波动的压力分布得到的。可以看出，无论是图 12a、12b的膜厚分布，还是图 12c、12d的压力分布，2个表面的计算结果都非常接近。值得注意的是，叠加了波动的表面在接触区域内的间隙明显比未叠加表面的大，这说明低弹性模量表面形貌的次级结构越粗糙，表面力的作用反而更显著，水合层的承载作用反而能得到加强。
将按256×256个节点划分的网格做2次抽除奇数点的操作，生成按64×64个节点划分的网格。对这个64×64网格进行表面力作用下的接触计算，与256×256网格的计算结果对比如图 13所示。比较发现，这两者的计算结果是比较接近的，只是稀疏的网格可能略微高估了膜厚值。有趣的是，稀疏网格造成的膜厚增加反而正好弥补了次级表面形貌造成的膜厚上升。
由此可以断定，节点间信息的省略对结果的影响很小，64×64网格下得出的结论是可信的。
3.3 弹性模量对计算过程的影响
4 结论
本文建立了考虑表面力的低弹性模量聚合物材料的接触模型数值计算方法，在聚合物弹性模量为0.2~20.0 GPa的条件下给出了一些算例，探究了表面力大小和弹性模量大小对接触区域的状态分布的影响情况，为今后的聚合物材料在表面力作用下的混合润滑计算奠定了基础。主要结论如下：
1) 静态接触条件下，低弹性模量材料表面被表面力完全分离。计算结果表明，低弹性模量材料接触压力较小，可以完全低于表面间作用力最大值，造成两表面间完全分离，表面间载荷全部由水合层承担，在表面间形成一层贯通的纳米级间隙；此外，由于接触区域大、弹性模量低，聚合物材料更不易发生塑性形变。
2) 低弹性模量材料的表面形貌对表面力的作用效果影响不大。计算结果指出，对于弹性模量低的情形，聚合物表面形貌对接触压力分布和表面力的作用效果影响很小，高弹性模量的下情形则相反。因此，希望应用水合效应的硬质材料摩擦副必须注意表面的加工，软质材料则不必要。对于低弹性模量材料，更粗糙的次级表面形貌会增加水合层的承载能力。
3) 低弹性模量显著增加了迭代算法收敛的难度。研究中发现，低弹性模量聚合物材料表面的次级细节不会显著影响接触计算的结果，因此可以采用较低的网格密度计算，以提高计算速度；材料弹性模量显著影响表面力迭代算法的收敛能力，弹性模量越低，算法就越难收敛，需要更严格的参数约束来辅助收敛。
下一步将研究算法收敛性、表面力接触与Reynolds方程的结合以及表面形貌分形。
