1　引言

异种材料连接尤其是金属和玻璃的连接，在电子工业、汽车制造等领域展现出广阔的应用前景。传统的金属与玻璃连接工艺包括润湿封接^［1］，阳极键合^［2］，胶接^［3］，半固体连接^［4］等。与传统方法相比，异种材料超快激光焊接具有精度高、焊接速度快、热影响区小、连接质量高等优点，近年来受到广泛关注。到目前为止，大多数研究都是实验性的，重点对焊接参数、焊接质量评估和焊接接头的微观表征进行了大量的实验和理论分析。Utsumi等^［5］使用5 ns脉冲激光直接把玻璃和金属结合起来，发现焊接区域旁边出现大的颗粒，随着脉冲能量增加和辐照次数增加，颗粒也越来越多，连接强度越来越大。Ozeki等^［6］使用飞秒脉冲实现了玻璃和铜基板的直接微焊接，在二者交界面处可能存在玻璃与铜的混合区。Zhang等^［7］利用1 kHz飞秒激光脉冲实现了铝硅玻璃与金属的直接焊接，阐述了焊接过程：电子吸收能量，将能量传递给晶格，晶格达到热平衡，之后扩散、熔化、产生相变和等离子体形成，产生的反冲压力将铜离子反射到玻璃缝隙中，之后冷凝形成连接。Carter等^［8-9］利用皮秒激光实现金属和玻璃焊接，并对焊接截面进行表征，发现交界面处有Si和Al元素5 μm左右的过渡区域，揭示超快激光焊接存在扩散现象。Carter等^［10］还比较了6082铝合金与石英玻璃和BK7玻璃的皮秒激光微焊接，当功率提高到一定程度时，提高功率并不能明显提高剪切强度，且由于热膨胀系数差异，6082铝合金与BK7玻璃的焊接强度优于石英玻璃。Li等^［11］用激光对钛合金表面氧化处理，实现钛合金和高硼硅玻璃透射焊接，处理后的连接强度提高5倍，在交界处出现4 μm Si-Ti元素过渡层。张敏等^［12］讨论了激光封接玻璃与可伐合金的连接机理，发现连接区域主要为Fe-Si化合物，其他黏附物起到辅助连接作用。此外，张国栋等^［13］实现了硅基底与硅硼酸玻璃飞秒激光焊接，在两种材料的交界面处，存在大约10 μm的Si元素过渡区域。而在理论上，超快激光与金属-玻璃界面的相互作用机制尚未完全揭示。

超快激光与材料相互作用的过程，使用有限元等传统方法难以模拟，而分子动力学常用于研究超快激光烧蚀、表面改性等，这将是模拟微观尺度超快激光与异种材料相互作用的理想方法。吴寒等^［14］构建了Ar-Al相互作用势能函数并给出了飞秒激光在氩气环境中烧蚀铝靶分子动力学模拟结果，喷射物由中等大小的团簇和单原子混合物组成，前端是低密度的气态原子，后端喷射物则来自材料内部烧蚀后的分解和扩张等慢过程。Cheng等^［15］使用势能函数（Beest Kramer van Santen potentials，BKS）研究了玻璃超快激光烧蚀的分子动力学过程，模型考虑了电子的产生和复合，烧蚀过程在材料内部模拟得到了空泡，而不考虑自由电子效应，则在材料内部模拟未发现空泡。陈冰等^［16］采用结合双温方程的分子动力学方法，数值模拟了飞秒激光与CuZr非晶合金的相互作用，结果表明超快激光作用下CuZr非晶材料中原子加热速度比普通晶态金属慢得多，且烧蚀深度随着能量密度增大而增加。Karim等^［17］开发了一个耦合到分子动力学的双温模型（TTM-MD）来模拟在固体透明覆盖层空间约束条件下短脉冲激光与金属的相互作用，在约束下更容易产生局部高温高压条件，减弱应力波的产生，从而减少熔池深度，抑制金属表面破裂和喷发。Karim等^［18］通过进一步的实验研究，证明分子模拟能够揭示材料在焊接过程中微观结构的演变。孙诗壮等^［19］使用分子动力学方法计算了Mo、Si原子发生反射和再溅射的概率，以及原子的反射、再溅射角度和能量分布，发现当沉积原子传递给基底的能量降低时，发生反射的概率增加，发生再溅射的概率降低。吴雪峰等^［20］分析了飞秒激光仿真所采用的双温模型、分子动力学模型及复合模型的特点及其适用范围，为飞秒激光加工的理论研究提供依据。侯超剑等^［21］基于分子动力学仿真分析了纳米颗粒钎料在SiO₂基底上的激光熔融过程，发现为了获得可靠的纳米互连结点，激光辐照下温度对纳米颗粒与衬底之间吸附能的调控是影响纳米互连结点稳定性的主要因素。周榆等^［22］和牛泽伟等^［23］通过双温模型耦合分子动力学方法研究了飞秒激光作用于金属纳米结构的物理过程。这些研究显示了分子动力学方法在模拟飞秒激光与异种材料作用中的可行性。

由于石英玻璃中原子间主要的相互作用力为共价键，其精确计算需要较高的计算成本，此时简化的相互作用势能函数就具有重要意义。本文初步实现飞秒激光作用石英玻璃和铝界面的模拟研究，构建了SiO₂粒子间相互作用Lennard-Jones（LJ）势能函数、铝和SiO₂相互作用LJ势能函数和飞秒激光作用铝和石英玻璃界面的分子动力学模型。所有模拟研究均使用分子动力学方法，该方法专注于系统中每个粒子的运动，揭示了超快激光作用异种材料界面的微观机制，提供了飞秒激光焊接异质材料的一般物理图像。

2　分子动力学模型的构建

使用由桑迪亚国家实验室开发的LAMMPS开源代码包^［24］，在广州国家超级计算机中心（NSCC-GZ）天河二号超级计算机上进行模拟求解，仿真域被划分为小子域并分配给每个处理器，执行并行分子动力学模拟。

2.1　SiO₂粒子之间的相互作用势

石英玻璃由具有复杂原子间相互作用和拓扑结构的无定形SiO₂组成，其中包含了氧原子-氧原子、氧原子-硅原子、硅原子-硅原子的相互作用势等，增加粒子模拟的复杂性，导致计算成本增加。因此对SiO₂粒子之间的相互作用势能函数进行简化，将SiO₂整体简化为单一粒子，建立SiO₂粒子之间的相互作用势函数。2个SiO₂粒子之间的相互作用势指定为LJ势^［25］，可表示为

式中：r_ij为2个SiO₂粒子之间的距离；σ为势能（PE）为零时2个分子之间的距离；ε为最小势能的绝对值。为了建立SiO₂粒子相互作用势函数，需指定LJ势中提到的分子质量m、σ和ε。这些参数可根据材料的部分宏观性质（密度ρ=2.2 g/cm²、熔化温度T_m=1938 K和弹性常数C₁₁=76.6 GPa）^［26］拟合得到。

通常，一种宏观性质与多个微观参数有关。因此，在拟合过程中应同时调整相关微观参数，以获得最佳理论值。此外，为简化模型，假定SiO₂粒子的晶体结构为面心立方晶格（fcc），具有LJ原子间势的fcc晶体的晶格常数为2^2/3σ^［27］。而m和σ之间存在精确的关系，即m=ρσ³。未确定参数由3个减少到2个，只需确定σ和ε，即可确定SiO₂粒子之间的LJ势。通过模拟SiO₂粒子的微观加热过程和弛豫后对晶格施加应变的过程，将微观参数与石英玻璃熔化温度T_m和弹性常数C₁₁建立联系，确定SiO₂粒子LJ势函数中的m、ε和σ参数值。

2.2　铝原子与SiO₂粒子之间的相互作用势

选择LJ势作为铝原子和SiO₂粒子之间的相互作用势函数。LJ势函数中参数σ为铝和SiO₂粒子平衡时间距为两种材料各自的粒子间平衡距离的算术平均值^［17］。LJ势函数中参数ε根据铝液滴和石英玻璃之间润湿过程中的黏附功拟合求解。实验上石英玻璃和铝在1073 K下的黏附功为0.844 J/m^2［28-29］，建立微观铝与SiO₂基材之间分子动力学模型模拟润湿过程，模型中使用的SiO₂粒子之间的势能函数为式（1）得到的LJ势函数，铝原子之间的势函数使用改进的嵌入原子势函数，该模型在激光烧蚀模拟中与实验结果较一致^［30］。模拟中使铝原子和SiO₂粒子分别在1073 K时弛豫，铝原子移动到石英玻璃基底表面，系统稳定在1073 K。铝液滴在SiO₂表面延伸并逐渐达到稳定状态，黏附功可以根据润湿面积和润湿过程前后的总能量差来计算。

2.3　飞秒激光作用玻璃-铝界面皮秒时间尺度的分子动力学过程

超快激光作用物质，在时间尺度大体可分为以下几个过程：在1~1000 fs的时间内，物质内产生自由电子、等离子体，随后电子与晶格进行能量交换，在1~100 ps时间内，出现相变、库伦爆炸、超临界流体等现象，随后在1~100 ns发生熔化流动、组织形变、表面形变、冷凝等过程^［31］。模拟频率为1 kHz，脉冲宽度为35 fs，能量密度为1.0 J/cm²的飞秒激光作用于石英玻璃和铝的交界面。在单脉冲作用后，下一个脉冲到来之前，1 ms时间内已完成了上述过程，因此模拟计算只考虑单个脉冲的作用，未考虑多脉冲累加效应。通过飞秒激光作用在石英玻璃和铝交界面的分子动力学过程模拟，小规模地定性揭示焊接物理现象，计算模型的示意图如图1所示。使用经典的分子动力学方法（MD）结合铝的双温度模型（TTM）来描述激光与铝的相互作用。石英玻璃部分只使用MD方法，因为脉冲能量吸收主要发生在铝区域。在传统的TTM中，有两个热传导方程分别描述了电子和晶格的热传导过程。而在TTM-MD方法中，晶格方程被分子动力学方程所取代，即晶格中的热传导是通过原子或分子运动实现的^［32］。

图 1. 石英玻璃-铝焊接区域分子动力学模拟示意图

Fig. 1. Schematic of the fused silica-Al molecular dynamics welding simulation

下载图片 查看所有图片

脉宽为35 fs、能量密度为1.0 J/cm²的飞秒激光脉冲沿x轴从左到右穿过石英玻璃区域，照射在铝和石英玻璃的界面上。模拟区域外设定为真空环境，y、z方向设定为周期性边界条件。铝区域的长度为120 nm，石英玻璃区域长度为80 nm，在y、z方向上宽度和高度均设置为10 nm。由于模拟区域y-z截面（10 nm×10 nm）取自实际焊接区域中的一块小区域，施加在模拟区域的激光可以认为是均匀分布的。此外，小尺度模型模拟会引起应力波在模拟边界的强烈反射，应力波大约需要20 ps传输到系统的左右边界。为避免应力波反射对两种材料界面的影响，仿真时间尺度设为30 ps。

3　模拟结果与分析讨论

3.1　SiO₂粒子之间相互作用势参数求解

模拟中调整SiO₂粒子LJ势函数中的σ和ε参数，计算不同参数时石英玻璃的熔点和弹性常数，使模拟的宏观性质与石英玻璃的实验结果足够接近，以此来确定LJ势函数中的参数值。不同的σ和ε参数值下模拟石英玻璃熔化温度T_m和弹性常数C₁₁值，如表1所示。

由表1可知，保持ε在0.200 eV，随着σ的增加，C₁₁呈下降趋势，逐渐从87.094 GPa减小到63.457 GPa，而熔化温度保持在1920±10 K。结果显示材料系数σ值越大，柔软度更高。保持σ在0.280 nm下，随着ε的增加，T_m和C₁₁均呈线性上升。根据上述参数组得到的结果，最终选择了参数σ=0.285 nm，ε=0.200 eV，此时模拟的宏观性质为T_m=1925.0 K，C₁₁=75.413 GPa，与实验结果T_m=1938.0 K，C₁₁=76.6 GPa比较接近。

3.2　铝原子与SiO₂粒子之间的相互作用势参数

通过模拟铝液滴和石英玻璃基底之间的润湿过程，确定铝原子与SiO₂粒子之间LJ势函数中的σ和ε值，浸润过程模拟图像如图2所示。

图 2. 铝液滴向石英玻璃移动后不同时间的润湿模拟图像（由VMD^［33］绘制）。（a）2 ps；（b）10 ps；（c）20 ps；（d）30 ps；（e）100 ps；（f）200 ps

Fig. 2. Snapshot of wetting simulation at different time after the aluminum is moved towards the fused silica (Drawn by VMD^[33]). (a) 2 ps; (b) 10 ps; (c) 20 ps; (d) 30 ps; (e) 100 ps; (f) 200 ps

下载图片 查看所有图片

由图2（a）~图2（c）可发现铝液滴在接触石英玻璃基底后开始在基底上扩展。而由图2（d）~图2（f）可发现铝液滴在SiO₂表面扩散到一定程度后，铝液滴保持稳定的形状，系统达到稳定状态。取SiO₂粒子与铝原子相互作用LJ势参数σ为0.2695 nm。为了计算宏观系统黏附功，需要计算润湿过程中释放的能量和最终润湿面积。为此研究了润湿过程中ε不同值下系统势能随时间变化，并将其绘制在图3（a），势能差随着ε的变化如图3（b）所示。

图 3. 润湿过程系统势能变化。（a）ε值在0.07~0.12 eV下系统势能变化；（b）不同ε值下系统势能差变化

Fig. 3. PE change during the wetting process. (a) PE change for different ε values during 0.07-0.12 eV; (b) ΔPE following the change of ε

下载图片 查看所有图片

由图3（a）可知，不同ε值下，系统势能随时间逐渐降低，然后在特定值上下波动，系统趋于稳定。随着ε值的增加，润湿过程中系统势能逐渐降低，即势能差增加。由图3（b）可知，随着ε值的增加，势能差呈线性增长。由于润湿区域的形状不规则，润湿面积采用近似方法计算。铝-石英玻璃界面被分成小网格，润湿面积可以通过计算铝原子覆盖的网格来计算。得到不同ε下润湿面积随时间变化如图4（a）所示，模拟得到的平均润湿面积随ε变化如图4（b）所示。

图 4. 润湿过程的润湿面积变化。（a）不同ε下润湿模拟过程的润湿面积变化；（b）不同ε取值对应的润湿面积

Fig. 4. Wetting area variation during the wetting process. (a) Wetting area variation during wetting simulations for different ε; (b) wetting area following the change of ε

下载图片 查看所有图片

图4（a）显示的是整个润湿过程开始后不久，润湿面积随着铝液滴与石英玻璃玻璃接触面积的扩大而快速增加，随后整个系统趋于稳定，润湿面积保持在某一特定值附近，与图2（d）、图2（e）和图2（f）的润湿图像保持一致。图4（b）浸润面积为系统稳定后润湿面积随时间的平均值，从中可以看出，润湿面积随着ε值的增加呈线性增加。最终，通过系统势能差和润湿面积比值计算不同ε值下的系统黏附功，如图5所示。图5显示系统的黏附功随着ε值呈线性增长，对二者的关系进行线性拟合以定量描述，根据实验上二者的黏附功为0.844 J/m²，确定ε为0.0970 eV。

图 5. 黏附功的拟合结果

Fig. 5. Fitting result of the work of adhesion

下载图片 查看所有图片

3.3　焊接模拟的结果

基于上述建立的SiO₂-SiO₂、铝-SiO₂相互作用势函数，模拟石英玻璃和铝的焊接微观过程，其中的铝原子之间的势能函数使用改进的嵌入原子势函数。向铝-石英玻璃界面施加强度为1.0 J/cm²的35 fs激光脉冲，得到前30 ps的粒子演化，如图6所示。

图 6. 分子动力学模拟铝-石英玻璃焊接过程图像

Fig. 6. Snapshots of the aluminum-silica glass welding process simulated by molecular dynamics

下载图片 查看所有图片

由图6可知，在1 ps时，铝区域完全无序，这意味着最初的铝内部结构被破坏，原子在随机剧烈热运动。接下来从1 ps到30 ps，铝原子产生的应力使得铝原子向石英玻璃内部扩散，两种材料的界面向石英玻璃一侧移动。在30 ps时，明显可以看出铝原子进入到石英玻璃区域，而SiO₂粒子也进入到铝中，界面附近的两种粒子由于剧烈的热运动而呈现出扩散趋势。此外，界面附近的SiO₂粒子变得无序，而距离界面约10 nm的SiO₂粒子依旧保持有序状态。整个系统空间上温度和应力随时间变化如图7（a）、图7（b）所示。

图 7. 系统的温度、应力分布。（a）不同演化时间下的系统温度随空间变化；（b）不同演化时间下系统y-z平面法向应力随空间变化

Fig. 7. Temperature and stress distribution of the system. (a) Temperature profile of the system at different time; (b) normal stress profile in the y-z plane of the system at different time

下载图片 查看所有图片

图7（a）显示激光脉冲作用在铝-石英玻璃界面（x=0 nm）后铝一侧区域温度迅速上升，温度在长约为100 nm的范围内保持稳定，而石英玻璃一侧区域由于热传导系数较低，温度上升相对较慢。铝一侧的温度，在1~10 ps迅速从1000 K左右增长到7500 K左右，而在10~30 ps仅从7500 K增长到10000 K左右，温度梯度变化随着时间逐渐降低。在石英玻璃一侧，在1~30 ps，温度稳定增长到2000 K左右，温度梯度变化先增加后减小。随着铝原子的热扩散，更高能量的铝原子进入到石英玻璃一侧，整个温度交界处随着铝原子的扩散，向石英玻璃一侧移动，温度主要集中在Al一侧100 nm的深度内。在文献［17］中也报道了相似的温度云图仿真结果，激光作用transparent overlayer-Ag的双层结构时，在皮秒时间尺度内，Ag温度上升到9000 K左右，温度主要集中在transparent overlayer-Ag交界面下100 nm深度内。

图7（b）显示的是y-z平面的法向应力σ_x随时间在不同空间上的演化。铝一侧的应力在1~10 ps，从4 GPa迅速增长到20 GPa，随后在10~30 ps，迅速降低到5 GPa左右。在激光作用的一瞬间，应力梯度非常大，应力在10 ps内迅速变化，之后应力从x=0 nm开始向Al内部传递并衰减，在时空上呈现传递特性。与文献［34］报道的飞秒激光作用Al表面产生的应力波时空传递特性相似，应力波峰值最高达15 GPa，随着向Al内部的传递而衰减。而在石英玻璃一侧交界面处，1~10 ps，应力也快速增长到10 GPa，随后10~30 ps下降到4 GPa左右。沿着玻璃内部，应力从交界面处向石英玻璃内部传递。铝-石英玻璃交界面处存在温度和应力梯度表明两种颗粒正发生剧烈的碰撞和混合。文献［17］报道了相近仿真结果，在激光作用于transparent overlayer-Ag交界面后，Ag内部应力最高可达30 GPa，交界两侧温度和应力存在梯度变化。

为了研究铝-SiO₂粒子的混合区域变化，将沿x轴方向同时包含铝原子和SiO₂粒子的区域定义为混合区域，来分析焊接过程中焊接区域的微观分子运动，混合区域的长度和中心位置如图8所示。由图8（a）可知，随着时间演化，混合区域一直在扩大，混合速率逐渐变大。在30 ps末，混合区域的长度增长到10 nm。另一方面，随着铝原子不断地向石英玻璃一侧碰撞，混合区域中心向石英玻璃区域移动，在30 ps时混合区域中心到了-20 nm，说明混合区域在玻璃-铝交界处整体向着石英玻璃一侧移动，并逐渐扩大，最终形成大面积混合区，宏观实验上可观察到SiO₂-Al的混合区^［9］，在其他异质材料焊接中也观察到类似的混合区^{［11，13］}。为研究混合区原子的状态，混合区域中铝和SiO₂粒子的温度和应力时空分布图如图9所示。

图 8. 混合区域几何参数随时间变化。（a）混合区域长度随时间变化；（b）混合区域中心位置随时间变化

Fig. 8. Geometric parameters of the mixed region over time. (a) Length of the mixed region over time; (b) central position of the mixed region over time

下载图片 查看所有图片

图 9. 混合区热力学性质云图。（a）铝原子温度云图；（b）SiO₂温度云图；（c）铝原子应力云图；（d）SiO₂应力云图

Fig. 9. Thermodynamic property contours in the mixed region. (a) Temperature contour plot of aluminum atoms; (b) temperature contour plot of fused silica; (c) stress contour plot of aluminum atoms; (d) stress contour plot of fused silica

下载图片 查看所有图片

很明显由图9（a）可知，混合区域的温度相当高，达到数千开尔文。激光脉冲作用后20 ps时铝原子的温度达到8000 K，在铝一侧的空间高温区域宽度随着时间逐渐扩展到10 nm。由图9（b）可知，混合区域中SiO₂粒子的温度在20 ps时达到了6000 K，在30 ps时达到10000 K，且整个高温区域宽度大约10 nm，其结果与图8结果保持一致。图9（c）、图9（d）显示混合区中铝原子部分的应力比SiO₂粒子低，但均达到了3 GPa以上，由图中可以看出混合区中高应力的SiO₂粒子都集中在较小的区域内，且偏向石英玻璃一侧，容易在玻璃内部积累应力，产生裂纹，使两种材料连接强度降低。

总之，经过飞秒激光脉冲后，整个模拟系统达到了高能量的非平衡状态。此外，铝中的温度和应力迅速增加到异常高的水平。铝中的能量通过粒子之间的相互碰撞传递给SiO₂粒子。而在界面附近的粒子由于率先吸收激光脉冲能量达到高能状态，导致两材料界面局部产生高温和强应力，从而使得两种粒子在界面附近扩散混合。

4　结论

构建了一个分子动力学模型，用于模拟飞秒激光作用在铝和石英玻璃界面的分子动力学演变过程。首先，基于玻璃的宏观性质，即单向拉伸模量和熔点，建立了石英玻璃的LJ相互作用势函数，确定其中σ=0.285 nm，ε=0.200 eV。其次，根据两种材料之间的黏附功，建立了铝和石英玻璃的LJ相互作用势。LJ势能函数的参数由σ=0.2695 nm，ε=0.0970 eV确定。最后，使用上面计算的参数模拟了小尺寸铝和石英玻璃的超快激光焊接微观过程。根据模拟结果，飞秒激光脉冲发射后，局部超小焊接区瞬时温度高达10000 K，高温区域大约10 nm。在两种材料内部，应力高达20 GPa，在两种材料混合区域，应力均达到了3 GPa以上。且在30 ps时，更多的铝原子进入到石英玻璃一侧，混合区域中心向着石英玻璃一侧移动，逐渐扩大，最终到了-20 nm位置。而后，两种物质的混合区由于相互侵入和高温高压的持续碰撞而不断扩大。这种现象揭示了两种材料结合的扩散机制过程，为超快激光异质材料焊接提供了理论基础。