基于分子动力学的密封副润湿特性对泄漏量影响研究

杨丽红 程 强

(上海理工大学机械工程学院 上海 200093)

基于应对能源危机以及环境保护等问题，在液压传动领域中，以海水或者淡水来替代传统的油液介质，大力发展无污染、绿色、清洁的水液压技术成为该行业领域研究的前沿方向。但由于在实际应用中水介质与油液介质的理化性能存在较大差异，无法直接采用油压设备的设计理论，极大制约了其发展，因此迫切地需要对水压相关理论研究进行完善。在液压传动中，采用间隙密封的传统油压设备依靠油液介质的黏度及表面张力等性质能够实现一定程度的自密封[1-2]。但当流体介质换成水后，由于黏度和表面张力等理化性质的改变，将对设备内泄漏造成重要影响。流体介质的内泄漏将会导致设备工作精度以及传递效率降低，最终影响设备工作的可靠性与安全性。

随着新型材料研究与制备的进展，研究人员发现不同润湿特性的壁面可以改变固-液界面上所产生的边界滑移长度，进而影响流体介质的流动阻力及流量大小[3-5]。CELATA等[6]通过在玻璃管道的内壁镀膜研究了疏水性通道对流体流动行为的影响，指出当壁面不可润湿(疏水)时，将会产生滑移条件。许少锋等[7]利用边界模型模拟研究了平板间的Couette流动，研究结果表明，当壁面与流体之间的斥力越大，即疏水性越强，所观察到的边界滑移现象越明显。狄勤丰等[8]利用纳米SiO2涂层制造出强疏水性壁面，证明岩石微孔道中产生了流体滑移，能大幅减小水流的流动阻力，增加水的流速和流量。虽然壁面润湿性对边界滑移有着重要影响在学术界已达成共识，但其作用机制和影响规律尚不明确，一些学者认为边界滑移只能在疏水性的壁面产生，而另一些学者却发现在完全润湿的壁面上也会产生细微的边界滑移现象[9-10]。

在微纳尺度下的间隙密封中，流体是否产生滑移对泄漏量有着重要影响。考虑密封副通常为金属材料，亲水性较好，本文作者从分子动力学的角度出发，通过改变力场中分子间范德华力的固液相互作用力参数，以不同润湿性的壁面进行剪切模拟，研究了在润湿状态下剪切作用对流体的边界滑移以及泄漏量的影响，旨在探索不同润湿性壁面对于密封防泄漏效果的影响，这对丰富水压领域以及间隙密封的理论研究与提高产品的可靠性、促进安全生产均具有重要意义。

1 模拟相关理论

1.1 分子动力学基本原理

分子动力学的研究方法最早可追溯至20世纪50年代。1957年ALDER和WAINWRIGHT采用硬球模型进行了最早的分子动力学模拟，这是一种基于牛顿运动定律来进行分子体系运动模拟的研究方法，通过在原子或分子层面上的模拟，从微观结构出发对实验现象进行推测和解释[11]。近几十年来计算机运算能力的不断提升以及势函数的日益完善，分子动力学研究方法被越来越广泛地应用于各种相关领域[12-15]。

牛顿运动方程是进行分子动力学模拟的核心，模拟的基本流程为：根据能量最低原理，首先给定模拟体系的初始构型，确定原子的初始位置；然后进行动力学弛豫，确保系统在开始模拟之前达到平衡状态；最后赋予原子速度和原子间的相互作用形式，即势函数，通过分析其受力情况，根据牛顿第二定律F=ma得到式(1)所示的原子的运动方程。

(1)

式中：mi为第i个原子的质量；ri为第i个原子的位置。

根据式(1)可以计算出原子的速度和加速度，指定积分步长后(通常是1～2 fs)，就能得到下一时刻原子新的坐标位置。依此重复计算下去，对比每一时刻原子的位置坐标，便可以得到整体的运动轨迹。

1.2 间隙密封的泄漏分析

液压缸作为液压传动系统的执行机构，在实际的生产工作中，由于间隙密封下的内泄漏问题无法从根本上避免，最终将导致整个传动系统的工作精度以及传递效率的降低。李琦[16]给出了泄漏量计算的详细推导及论证过程，在理想情况下，活塞与缸筒的偏心度ε为0。泄漏量计算公式为

(2)

式中：Δp为两端压差；μ为流体的动力学黏度；Q为总流量；d为活塞直径；h0为间隙大小；L为密封长度；v为活塞与缸筒的相对运动速度。

式(2)中等号右边第一项为压差泄漏量，第二项为剪切泄漏量，当运动方向与压差方向一致取“+”号，反之取“-”号。在工作过程中，由于活塞受压差所驱动，因此运动方向始终与压差方向相同，上式(2)中恒取“+”号。即：

(3)

在活塞运动的往复行程中，当内泄漏过大时，将造成液压缸的容积效率下降，达不到运行工况要求。虽然在反行程中流体介质逐渐朝密封副正行程相反一侧泄漏，进行一定的补偿抵消作用。但在正反2个行程中，其内泄漏均会对液压系统的精确控制产生不利影响，且传动误差将在活塞到达单向行程终点处达到最大值。

朱海燕等[17]指出活塞静止时其泄漏量与压差正相关，而活塞运动时，泄漏量将随着活塞运动速度而增大。考虑到由于密封副两端内泄漏问题引起的传动误差及效率传递误差将在活塞单向行程的终点处达到最大，因此文中针对单向行程中活塞在不同运动速度下的泄漏量进行了模拟，研究了表面润湿性对其所产生的影响，对于减少间隙密封设备在运动过程中的效率损失以及提高其工作可靠性有一定参考作用。

2 模拟过程及结果分析

2.1 模拟体系的建立

为模拟不同润湿性壁面对单向剪切泄漏的影响，考虑影响润湿特性的主要因素——固液界面相互作用[18]，文中以铝单质晶体的原子参数搭建了壁面结构，晶胞大小为边长0.404 95 nm的立方体，通过改变固液相互作用参数，得到不同润湿性的壁面。以SPC/E的水分子模型填充2个壁面的间隙，单个水分子的模型如图1所示。

图1 SPC/E水分子模型

为避免在晶面上产生“有序水分子单层”导致反常的疏水现象[19]，且FCC晶体各晶面表面能大小为FCC(1 1 0)>FCC(1 0 0)>FCC(1 1 1)，因此用(1 1 1)晶面作为表界面进行模型搭建。FCC晶体原子的空间排列及(1 1 1)晶面如图2所示。

图2 FCC晶体的原子空间排列及(1 1 1)晶面

2.2 势能模型的选取

水分子采用SPC/E刚球势能模型，其参数如表1所示。

表1 SPC/E水分子势能参数

SPC/E模型的势能函数包含两项，分别为范德华力和库仑力，总势能：

U(r)=UVDW+Ucoulomb

(4)

对于范德华非键结合势能，采用L-J势最常用的(12-6)势：

(5)

式中：r是2个粒子之间的距离；ε为势阱深度，等于势函数在极小值点rmin=21/6σ时的值，σ为当作用势为0时原子间的平衡距离，这2个参数是由粒子的种类所决定的[20]。

库仑势为

(6)

水分子与壁面之间的相互作用力也采用L-J势描述，对于不同的粒子之间的相互作用，参数可通过Lorentz-Berthelt混合规则来计算：

(7)

通过改变L-J势中的参数值可以调整固液相互作用力的大小，进而得到润湿性不同的壁面。

2.3 壁面润湿性的接触角表征

接触角是壁面润湿性的一个重要表征手段，因此文中针对上文所建立的不同壁面，模拟了水滴团簇在壁面上的润湿过程，采用COMPASS力场生成包含500个水分子的纳米水滴，壁面大小为10 nm×10 nm，在Z方向设置80 nm的真空层防止Z方向周期性边界条件的干扰，将水分子层放入，得到接触角的模拟体系如图3所示。

图3 接触角模拟体系(红色为氧原子；

对上述模拟体系用Forcite模块的Geometry optimization，选用Smart方法对结构进行优化，得到体系收敛时能量最低的状态，用Dynamics模块进行动力学计算。为了避免对水分子的转动及扩散产生影响，选用Nose-Hoova热浴法进行温控[21]，模拟温度为298 K，远程静电力(Electrostatic)采用PPPM求解器，范德华力(van der Waals)和氢键作用力(Hydrogen bond)采用Atom based方法，截断半径分别为1.25和0.3 nm，在NVT系综下进行50 ps的动力学驰豫，积分步长1 fs，共计50 000步，得到结果如图4所示。

图4 铝原子壁面接触角(红色虚线表示氢键)

由宽高法测量接触角[22]，计算公式见式(8)，接触角测量原理如图5所示。

(8)

图5 接触角测量方法示意

经测量得到接触角约为87°，与文献[23]中得出的水滴在纯铝上的接触角为82.7°基本一致。修改模拟体系中的固液相互作用参数，通常分子间作用力小于氢键力[24]，以固液作用力分别为5.0和8.8 kJ/mol进行接触角模拟，得到的模拟结果如图6所示，接触角分别为74°和41°。

图6 不同润湿性壁面的接触角

2.4 间隙密封单向剪切泄漏模拟

以上述不同润湿性的壁面参数进行剪切模拟，模拟体系如图7所示。

图7 剪切泄漏模拟体系

模拟体系大小为8.066 0 nm×4.946 2 nm×4.871 5 nm，水分子层中包含4 000个水分子，选用NVE[25]系综，在298 K的温度下进行模拟，剪切的行程为4.0 nm，分别在不同剪切速度以及不同润湿性壁面下进行模拟。模拟前先进行结构优化及动力学弛豫，体系达到稳定后进行剪切运动，水分子的运动如图8所示。

图8 水分子在剪切作用下的运动

剪切行程完成后分别统计了每组模拟条件下的最大剪应力、边界滑移长度以及模拟盒子外泄漏的水分子数，结果如图9—11所示。

图9所示为最大剪应力随壁面速度变化。可以看出，在剪切过程中，壁面所承受的最大剪应力大小受壁面的润湿性影响不大，而与剪切速度呈线性关系，随着剪切速度的提高，最大剪应力也相应增大。

图9 不同接触角时最大剪应力随壁面速度变化

图10所示为接触角分别为41°及87°时壁面表层水分子的滑移长度。可以看出，当剪切速度较高时，即使在润湿性良好的壁面上依然会产生细微的边界滑移现象，并随着壁面剪应力的提高不断增大；并且在润湿性较差的壁面上更易产生滑移现象，其滑移长度也更长，随着所受的剪应力继续增大，其滑移长度将趋向于稳定。

图10 接触角为41°及87°时边界层滑移长度随壁面速度变化

图11所示为接触角分别为41°及87°时润湿性壁面进行剪切模拟时的水分子单向泄漏量。可以看出，在接触角为87°时，活塞在不同运动速度下单向泄漏的水分子数量总体变化不大，处于波动状态，占总体水分子数的3.98%～4.88%；而当壁面较亲水时，泄漏量呈现出随着剪切速度增加而不断上升的趋势，2种壁面之间单向泄漏量差值在0.6 nm/ps的剪切速度下达到了水分子总数的2%。这意味着在高速剪切作用下，壁面润湿性对泄漏存在较大影响，且润湿性较差的壁面对防泄漏效果更好。

图11 接触角为41°及87°时泄漏量随壁面速度变化

实际生产中，由于受到发热量及制造水平的限制，剪切速度通常较低。为了进一步探究低速剪切下的单向行程泄漏情况，以更低的剪切速度进行模拟，发现在壁面的剪切速度为0.005 nm/ps时，其最大剪应力降低至0.885 GPa。此时最大剪应力不足以引起流体产生滑移，流体在几种固液边界上基本不产生滑移且泄漏量相当，由于剪切引起的单向行程泄漏量降低到了3%以下。

综上可知，在低速剪切时，亲水壁面上固液作用力较大，水分子易吸附于壁面之上，边界层基本不产生滑移，从而导致更大的泄漏量；随着剪切速度的提高，最大剪应力不断增大，在润湿性较差的壁面上，表层水分子受到内层水分子氢键力的作用，边界滑移长度与泄漏量逐渐趋向稳定，总体变化不大，呈现出波动状态，对于防泄漏效果更好。

3 结论

(1)分子动力学的模拟结果表明，在极高速的剪切作用下，流体在亲水性壁面的固液交界处仍然会产生细微的边界滑移现象。其滑移量的大小取决于壁面的润湿性程度以及剪切时的剪应力。

(2)剪切过程中的最大剪应力与壁面润湿性程度无关，与剪切运动速度呈正相关。在同样的剪应力作用下，润湿性较差的壁面更容易产生边界滑移，随着剪切速度的进一步增大，边界滑移的长度以及泄漏量将趋向于稳定值，而润湿性良好的壁面所能承受的极限剪应力较大，在不超过极限剪应力的范围内，流体的泄漏量将随着剪切速度的提升而增加。

(3)在实际生产中受限于制造水平及发热量的限制，剪切速度通常较小，难以产生极大的剪应力。模拟结果表明，当活塞运动速度为5 m/s时，最大剪应力大小为0.885 GPa，其剪应力不足以产生边界滑移，不同润湿性壁面下的水分子吸附于壁面之上，因此其单向行程泄漏量相当，此时对于提升液压设备工作的准确性与可靠性可以通过减小密封间隙或者降低活塞运动速度等实现。而在高速剪切运动时，壁面的润湿性对于泄漏量有着重要影响，除了上述手段之外，还可以通过在密封副上采用涂镀疏水涂层，降低其润湿性能，以此减小泄漏量，提升设备工作的精确度以及增加传递效率。

(4)文中以纯水作为流体介质进行了单向行程下由于剪切作用引起的泄漏研究，并得到了以上相关结论。但文中研究仍然存在不足之处，需要进一步进行研究，主要总结如下：在实际水压设备中大都采用资源充足的海水或者自来水，其中存在复杂的微量元素，海水中还存在大量矿质元素、钙镁化合物等，这些成分可能在一定程度上将对流体介质的理化性质产生影响，其成分以及含量对于泄漏量的影响行为有待进一步研究。