单晶硅表面纳米结构润湿行为的分子动力学模拟

王浩杰，曹自洋，张洋精，朱译文

（苏州科技大学 机械工程学院，江苏 苏州 215009）

表面润湿性是一种液体在一种固体表面铺展的能力或倾向性，表面润湿性是固体表面的重要特性。接触角一般用来表征表面的润湿性。表面润湿性是由固体表面的表面自由能和表面粗糙度共同影响决定的，其中固体表面的几何构型可以通过人为加工而进行改变[1-2]，因此表面微结构对于润湿性的影响特别值得关注。

近些年来，表面微结构对于润湿性的影响已经进行了大量的研究，Wenzel和Cassie等人提出的Wenzel理论[3]和Cassie理论[4]是现在表面润湿性主要的理论，润湿性的接触状态分别对应Wenzel接触状态和Cassie接触状态两种润湿状态。Yoshimitsu等人通过在硅片的表面上构制规则的微槽阵列和方柱阵列，并且通过控制硅片表面微结构的凸台宽度，来控制硅片表面粗糙度因子的大小，以此来比较微小水滴在两种微结构上表面接触角的实验测量值和理论计算值[5]。然而，接触角测量实验目前以微米级别为主，而对于纳米级别的水滴以及纳米级别的微结构，测量实验非常难精准。分子动力学模拟方法是一种计算机模拟固体材料表面微观性能的方法，分子动力学模拟方法使得可以观察测量纳米级别的接触角以及对于固体表面纳米结构得以控制。近几年借助分子动力学，纳米级别表面微结构对于润湿性的影响的研究成果相继出现。Chen等人利用分子动力学模拟的方法，研究出了纳米水滴在纳米阵列表面上存在的四种润湿模式，分别为Wenzel接触模式、亚稳定的Cassie接触模式、稳定的Cassie接触模式及混合润湿接触模式[6]。颜笑等通过分子动力学模拟方法研究了纳米水滴在纳米表面沟槽结构上的接触状态，研究了不同固-液相互作用下，纳米水滴的表面接触角和润湿模式随沟槽的宽度和沟槽的深度的变化规律[7]。上述研究中，微结构以光栅、方柱等基础常规结构为主，对于其他特殊微结构，如二级结构、凹槽结构、梯形结构未曾考虑且没有进行实验验证。

本文拟采用MD模拟技术，利用lammps软件，研究纳米液滴在不同形貌的单晶硅纳米表面结构对表面润湿性的影响。

1 单晶硅纳米微结构表面的构建与接触角理论模型

1.1 接触角理论模型

材料表面自由能平衡状态下存在着液体表面张力和固体表面张力，固液表面的接触角由两个表面张力决定，此接触角模型满足了Young's方程[8]：

其中，σsv、σsl、σlv分别为固相-固相、固相-液相和气相-液相之间的界面张力；θY是水滴在光滑固体表面上，固、液、气在三相平衡的时候的接触角，即本征接触角。图1所示为结构示意图：

图1 液滴张力在理论光滑表面上结构示意图

不过在真实情况下，固体表面上可能会出现缺陷、灰尘等，所以不存在理论上光滑的固体表面。Wenzel接触角理论模型[3]和Cassie-Baxter接触角理论模型[4]用来描述水滴在粗糙固体表面上的润湿状态。Wenzel理论模型和Cassie-Baxter理论模型接触角示意图如图2所示。

图2 Wenzel模型与Cassie-Baxter模型接触角示意图

Wenzel理论模型的方程为

Cassie-Baxter理论模型的方程为

式中，r是粗糙度因子，即实际固液接触面积与液滴水平投影面积之比。在Cassie-Baxter理论模型中f是相面积分数，是液滴和固体表面接触面积所占交叉界面的比值。

1.2 单晶硅纳米微结构表面的构建

为了更加充分的研究单晶硅不同表面纳米结构的润湿行为，本文在单晶硅表面构建不同的粗糙壁面。充分考虑未来的可加工性，分别构建了：光栅形、方柱形、凹槽形、梯形、二级结构，对应结构如图3所示[9]。

图3 光栅形、方柱形、凹槽形、梯形、二级结构形貌与尺寸示意图

其中，a为凸台宽度，b为凸台间的间距，h为凸台宽度，c为梯形结构下底长，l为梯形侧边长；在二级结构中，底部凸台宽度、凸台间距、凸台高度分别为a1、b1、h1，二级凸台宽度、凸台间距、凸台高度分别为a2、b2、h2。各个结构的单位为单晶硅的晶格常数，即一个晶胞的边长。

通过计算光栅形、方柱形、凹槽形、梯形、二级结构粗糙壁面的粗糙度因子r与结构参数的关系，如式（4）-（8）所示。

通过计算光栅形、方柱形、凹槽形、梯形、二级结构粗糙壁面的其相面积分数f与结构参数的关系，如式（9）-（15）所示。

经过计算，可以得到五种不同结构的r，f的值如表1所列。

表1 不同结构r、f的计算值

2 分子动力学模拟模型和方法

Choong等人研究发现，单晶硅100晶面相较于单晶硅111和110晶面，在相同加工条件下，加工出的结构更加显著，能够更好地进行研究[10]。因此，本文的研究对象选择了单晶硅100晶面，在模拟过程中，选取了相同的单晶硅100晶面进行研究，在单晶硅表面构建不同的结构。模拟体系的初始构型如图4所示，模拟盒子尺寸为86.9 nm×86.9 nm×148.5 nm，2 048个水分子随机分布成纳米水滴，水分子的初速度由随机数发生器确定。对单晶硅晶面进行剪切来构造单晶硅100晶面，对剪切后的单晶硅晶胞进行超胞命令，构成了十六个原子层厚度，此外单晶硅晶胞在XY方向上进行扩张延展，形成了面积为7 551 nm2的平面。

图4 模拟体系的初始构型

本次分子动力学模拟仿真中，系统的初始温度设定为接近于普通室温的298 K，选取了NVT正则系综进行模拟，即在恒定的温度、恒定的粒子数、恒定的环境和一定的体积下进行分子动力学仿真。采用Nose-Hoover温度耦合算法维持体系温度以及压力的恒定。本研究在模拟空间的X、Y方向上，选择并且施加周期性边界条件，在模拟空间的Z方向上，对于边界施加镜面反射边界条件，在仿真模拟过程中设置的步长为2.0 fs，仿真时间为1 000 ps。保证系统达到平衡。模拟数据采用LAMMPS软件来进行计算，接触角确定方法参考文献[11]。

本文模拟的水滴采用能够充分表征溶液体系的水分子结构的SPC模型，SPC的势能参数如表2[10]所列。

表2 SPC模型的参考值

Lennard-Jones双体势能用于水分子中的氧原子之间，考虑到短程排斥和远程吸引力。LJ作用势是一种短程作用关系，一般用来描述液体和气体之间的相互作用以及范德华力占主导地位的材料。Tersoff多体势函数描述了单晶硅中Si原子的选择。 Si原子和水的氧原子之间的势能也用LJ势函数来描述。总势能由两部分组成，短程LJ势能和长程静电势能。如式（14）[11]所示。

水分子中O原子与壁原子的相互作用仍采用LJ的12-6势能形式，如式（15）所示，相应的参数可以利用Lorentz-Berthelot混合规则确定，如式（16）[12]所示。

3 结果与讨论

3.1 分子动力学模拟结果

通过MD模拟，图5为系统稳定后水滴在平整表面的情况，此时水滴近似为球形，分子模拟的接触角为72.02。在实际测量中，被抛光后的单晶硅表面液滴的测量值为67.3。

图5 平整表面模拟结果

通过MD模拟得到的纳米液滴在五种形貌粗糙壁面上的接触状态及接触角，如图6和表3所示。

图6 不同形貌粗糙壁面的接触角状态

表3 不同形貌粗糙壁面的接触角

3.2 接触角预测模型的分析与讨论

通过式（2）和式（3），根据构建的两种模型进行数值分析，可以得到两种模型下的接触角数值，实际表观接触角与两种模型的接触角数值如图7所示。

图7 不同纳米微结构接触角的关系图

可以看到在分子层面上，具有纳米微结构的单晶硅表面疏水性优于平整表面，具有二级结构的单晶硅表面疏水性能优于其他结构，凹槽结构表面疏水性能最差；具有纳米微结构的单晶硅表面的润湿行为更加趋向于Cassie-Baxter模型。

4 结论

本文采用分子动力学模拟技术，探讨了单晶硅表面不同纳米结构粗糙表面润湿行为，通过表面形貌分析和接触角预测模型的分析与讨论，得出以下结论：

（1）在纳米尺度下，构造了纳米微结构的单晶硅，其表面疏水性能增强，从平整表面的亲水性转化为疏水性。

（2）通过对单晶硅不同纳米微结构表面润湿性的分子动力学模拟，单晶硅的二级结构表观接触角测量值最大，凹槽结构表面表观接触角测量值最小。

（3）单晶硅表面不同纳米结构粗糙表面润湿行为趋向于Cassie-Baxter空气柱模型。不过，由于忽略了动能、湿度等外在因素对于液滴接触角的影响，导致模拟数据与理论数据存在部分误差。