1. 项目概述从原子视角窥探玻璃的“筋骨”玻璃这种我们日常生活中再熟悉不过的材料其本质是一种非晶态固体。与金属或晶体不同玻璃内部的原子排列没有长程有序的周期性呈现出一种“冻结的液体”状态。这种无序性赋予了玻璃透明、各向同性等独特性质但也使其在力学性能上尤其是断裂行为充满了复杂性和不确定性。为什么有些玻璃一碰就碎而有些却能承受巨大的弯曲这个问题的答案深埋在原子尺度的微观世界里。传统的材料研究多依赖于宏观实验但面对非晶材料这种“黑箱”宏观数据往往难以揭示其内在的物理机制。分子动力学模拟就像一台超级显微镜允许我们追踪每一个原子的运动轨迹在计算机中“制备”出具有特定微观结构的材料模型并对其施加虚拟的载荷从而直观地观察从弹性变形到最终断裂的全过程。本次分享的项目正是运用这一利器聚焦于一个经典的模型体系——二维二氧化硅玻璃。选择二维体系并非为了简化而是为了剥离三维空间的复杂性更清晰地聚焦于非晶结构本身对力学行为的影响。我们构建了三种具有代表性的微观结构模型纳米晶玻璃、准晶玻璃和连续随机网络玻璃旨在系统性地回答一个核心问题原子尺度上的无序结构差异如何一步步决定材料的宏观断裂韧性对于材料科学、物理或相关领域的从业者无论是初涉模拟的研究生还是希望深化对非晶材料理解的技术人员这篇文章将带你深入一场从原子建模到宏观性能解读的完整旅程。我们将不仅复现模拟的关键步骤更会拆解每一个结果背后的物理图像分享参数设置背后的考量以及如何从海量的原子轨迹数据中提炼出有价值的科学见解。2. 核心模型构建与结构表征解析分子动力学模拟的第一步也是决定模拟成败与物理意义的关键一步就是构建能够真实反映材料特性的原子模型。对于二维二氧化硅其基本结构单元是硅氧四面体在二维平面上的投影即一个硅原子与三个氧原子键合形成的三角形单元这些单元通过共用氧原子连接成网络。2.1 三种非晶结构模型的生成策略我们构建了三种截然不同的二维二氧化硅玻璃模型它们代表了非晶结构中从“有序”到“无序”的谱系。2.1.1 纳米晶玻璃模型顾名思义NCG模型由许多微小的晶粒镶嵌在非晶基质中构成。我们的生成方法是首先在模拟盒子中随机撒播多个晶核点每个晶核以二维β-鳞石英的晶体结构一种二氧化硅的晶体相为模板进行生长。通过控制晶核的数量和生长速率我们可以获得不同晶粒尺寸和体积分数的纳米晶结构。当晶粒生长到相互接触或达到预设尺寸后剩余的空间则通过快速淬火即高温熔体快速冷却的方式形成非晶基质。这个过程模拟了某些通过控制结晶工艺制备的玻璃陶瓷材料。注意晶粒与基质的界面是关键。在模拟中需要确保界面处的原子通过能量最小化进行了充分弛豫以避免引入不切实际的高能界面从而影响断裂起始位置的判断。2.1.2 准晶玻璃模型Para模型比NCG更加无序但保留了中程有序的特征。我们采用了一种“拓扑约束”生长算法。首先生成一个完全随机的初始原子构型然后通过引入一组约束条件如键长、键角的目标分布函数这些函数来源于对实验或高质量CRN结构的统计进行迭代优化。在优化过程中系统被允许形成一些局域的、类晶的短程有序团簇但这些团簇之间没有明确的晶界呈现出一种“近程有序、中程弥散”的状态。这类似于某些通过特殊工艺制备的、结构高度均匀的非晶合金。2.1.3 连续随机网络模型CRN模型是描述理想玻璃态最经典的模型。它假设网络是连续、无限且完全无序的但每个原子都满足确定的配位数硅为4氧为2。我们采用Wooten-Winer-Weaire算法或其变体来生成。该算法从一个晶态种子开始通过一系列随机的键切换操作来引入无序同时严格保持网络的连通性和配位数约束直到系统的能量和径向分布函数收敛到典型玻璃态的特征。最终得到的结构没有任何周期性是三种模型中最无序的。2.2 关键结构表征量的计算与物理意义模型建好后如何定量描述它们的结构差异这就需要借助一系列结构表征函数。2.2.1 径向分布函数探测短程有序的“尺子”RDF是材料科学中最基础的结构表征工具。对于我们的体系我们计算了Si-Si、Si-O和O-O的RDF。计算原理是以某个原子为中心统计在距离r到rdr的球壳内发现其他原子的平均数量然后对体系中的所有原子和所有构型进行平均并归一化。第一峰位置对应最近邻原子间的平均距离。例如Si-O RDF的第一峰位置直接给出了Si-O键的平均键长这是判断势函数合理性和结构弛豫程度的重要指标。第一峰面积积分后可以得到配位数。Si-O第一峰面积应接近4Si的配位数O-Si第一峰面积应接近2O的配位数。峰的宽度和衰减峰越尖锐说明键长分布越集中有序度越高。从第二、第三峰开始对于晶体RDF会出现清晰的峰对于CRN玻璃这些峰会迅速展宽并湮灭这正是长程无序的体现。将模拟CRN的RDF与实验测得的RDF如图1e,f进行对比是验证模型可靠性的黄金标准。2.2.2 环分布统计揭示中程有序的“拓扑地图”在二维网络中环由Si和O原子交替连接形成的闭合多边形是描述拓扑结构的核心。我们统计了网络中所有最小环的边数通常是4-10元环。对于完美的β-鳞石英晶体只存在6元环。而在非晶网络中会出现4、5、7、8元环等多种拓扑缺陷。环分布图如图1d所示NCG模型由于含有晶粒其环分布会在6元环处出现一个尖锐的峰同时在其他位置有非晶基质贡献的宽分布。Para模型的分布较宽但可能在6元环附近仍有凸起。CRN模型的分布最宽且平缓是典型的无特征分布。环分布的差异直接影响了网络的刚性、自由体积分布进而影响力学性能。2.2.3 角度分布函数与结构相似性参数ADF统计了诸如Si-O-Si和O-Si-O等键角的分布。对于二氧化硅Si-O-Si键角分布通常在120°-180°之间中心约在145°。ADF的宽度反映了网络键角扭曲的程度。结构相似性参数s则是一个更全局的量度用于量化局部原子环境与某个参考结构如理想的晶体环境的差异其概率分布函数图1g可以直观显示体系内结构异质性的程度。通过这些细致的结构表征我们不仅在视觉上图1a-c更在定量上严格区分了三种模型为后续关联结构与性能奠定了坚实的基础。模拟得到的CRN的RDF和ADF与实验数据高度吻合图1e,f,h,i这让我们对模拟的可靠性充满信心。3. 分子动力学模拟流程与关键参数设置有了可靠的原子模型下一步就是让它们“动”起来接受力学测试。分子动力学模拟的本质是数值求解牛顿运动方程其流程和参数设置直接决定了模拟的物理真实性和计算效率。3.1 势函数选择描述原子间相互作用的“宪法”对于二氧化硅体系原子间的相互作用通常由经验势函数来描述。我们选择了基于BKS模型的修正势函数。BKS势是一种经典的二体势包含了长程库仑相互作用和短程Buckingham势的形式。E_ij q_i q_j / r_ij A_ij exp(-B_ij r_ij) - C_ij / r_ij^6其中第一项是库仑项q_i为原子电荷第二项和第三项是短程排斥和吸引项。我们采用Wolf方法对长程库仑力进行截断和修正以在保证精度的前提下大幅提升计算速度。实操心得势函数的选择是分子动力学模拟的基石。对于二氧化硅除了BKS还有TTAM、CHIK等势函数。选择时需权衡精度与速度。BKS势在描述玻璃态结构如RDF、密度方面表现良好且计算效率较高适合进行大规模断裂模拟。务必在正式模拟前用小体系测试所选势函数能否复现实验已知的键长、键角、弹性模量等基本性质。3.2 系统弛豫与平衡态获取在施加载荷前必须让初始构建的模型达到平衡态即体系能量和压力等宏观量在平均值附近小幅波动。能量最小化首先对初始模型进行共轭梯度法能量最小化消除原子间的过度重叠和不合理的高应力。NPT系综弛豫在恒温恒压系综下运行一段时间通常数十皮秒让体系密度弛豫到平衡值。温度通过Nosé-Hoover热浴控制压力通过Parrinello-Rahman控压方法控制。对于二维体系需注意控制面内压力而垂直于平面方向第三维的压力设为0。NVT系综弛豫在平衡体积下切换至恒温恒容系综继续弛豫确保温度均匀体系处于稳定的玻璃态。我们通常将温度设定在远低于玻璃化转变温度Tg的数值如300K。关键参数弛豫时间需要足够长。我们通过监测体系总能量、温度、压力以及RDF是否不再有漂移来判断平衡。对于含有数千原子的体系NPTNVT弛豫总计需要至少100-200皮秒。3.3 单轴拉伸模拟的实现细节我们采用位移控制的方式进行单轴拉伸模拟这是研究断裂行为的常用方法。加载方向沿模拟盒子的x方向armchair方向施加应变。同时允许y方向自由收缩泊松效应而z方向二维平面外保持固定厚度。应变率这是一个至关重要的参数。分子动力学模拟受限于计算资源使用的应变率通常为10^8 - 10^9 /s远高于实验速率10^-3 - 10^0 /s。如此高的应变率会抑制原子的热激活运动可能导致材料表现得比实际更脆。为了部分抵消这一效应我们采用了相对“慢”的应变率在MD尺度下例如5×10^8 /s并通过与不同应变率的对比测试评估应变率效应的影响。模拟步骤在每个模拟时间步通常为1飞秒将模拟盒子在x方向的长度增加一个微小量ΔL ε˙ * L_x0 * Δt其中ε˙是应变率L_x0是初始长度Δt是时间步长。然后重新标定所有原子的x坐标仿射变形接着在NVT系综下运行若干步让原子弛豫到新的构型并计算维里应力。应力计算体系在变形过程中的应力张量通过维里定理计算包含了动能项和原子间相互作用力的贡献。我们主要关注沿拉伸方向的轴向应力σ_xx。关键技巧为了捕捉断裂起始的精确位置我们在整个拉伸过程中以很高的频率如每100步输出原子的位置、速度和应力数据。这些海量的轨迹文件是后续所有分析的源头。4. 断裂起始机制从最弱环节到空洞形核材料不会无缘无故断裂断裂总是从最薄弱的环节开始。通过分子动力学模拟我们可以像侦探一样在断裂发生前就锁定这些“嫌疑位置”并观察“犯罪现场”是如何形成的。4.1 维里系数预测断裂起始点的微观应力探针宏观应力是体系整体的平均值而断裂往往由局部极高的应力集中引发。维里系数V_θ是我们定义的一个局部应力度量。对于每个原子i我们计算其在一个局部区域例如其最近邻壳层内的应力张量然后通过坐标变换得到沿拉伸方向θ的正应力分量V_θ(i)。这个值直观反映了该原子所在局部区域被拉伸的紧张程度。如图2a所示在未变形的平衡结构中我们就可以绘制出V_θ的空间分布图。那些呈现高正值红色区域的原子团簇意味着即使在零应变下该处的局部网络也处于一种“预拉伸”的亚稳态。它们就是潜在的断裂起始点。我们的模拟结果惊人地显示最终发生空洞形核的位置图2a中圆圈标注与这些高V_θ区域高度重合。4.2 空洞形核的临界条件与结构关联我们统计了所有模拟中空洞形核时的应变值ε_nuc。分析发现空洞形核应变与平衡结构中的两个关键结构参数存在强烈的相关性与最拉伸状态的关系对于每个样本我们找到其平衡结构中V_θ的最小值即最拉伸的状态。如图2b所示V_θ_min与ε_nuc之间存在近似指数衰减的关系。V_θ_min越小即初始就越拉伸材料越早发生空洞形核韧性越差。这为预测材料的断裂起始韧性提供了一个基于初始结构的微观判据。与Si-O-Si键长的关系我们测量了所有Si-O-Si桥氧键的长度即两个硅原子之间的距离。统计发现空洞往往在那些具有较长Si-O-Si键的附近形核图2d。较长的键通常意味着键能较弱或者该处的网络拓扑结构较为松散自由体积较大在拉伸下更容易被拉开。注意事项空洞形核是一个统计过程。即使在同一模型中由于原子级结构的随机涨落每次模拟的形核位置和应变也可能不同。因此我们需要对每种结构进行大量如50-100次的独立模拟并统计其累积概率分布函数CDF如图2c所示。CRN的CDF最宽说明其形核应变的波动性最大这与它结构的均匀性没有明显的优先形核点有关而NCG的CDF较窄因为晶界作为明确的缺陷总是优先形核。4.3 断裂事件的原子尺度快照通过分析断裂瞬间的原子构型图3d, g, j我们可以清晰地看到空洞形核的微观过程在NCG中空洞几乎无一例外地在晶界处形核。晶界处原子排列不规则键合较弱是天然的弱点。在Para中空洞在那些中程有序被破坏的局部松散区域形核。在CRN中空洞形核位置看起来更随机但追溯其平衡结构总能找到对应的V_θ高点或长键位置。形核后空洞并非立即导致灾难性断裂。它首先会稳定生长直到达到一个临界尺寸。5. 裂纹扩展动力学与韧性差异分析空洞形核只是断裂的序幕随后的裂纹扩展阶段才真正决定了材料的断裂韧性和断裂形貌。这是三种结构模型表现出最大差异的环节。5.1 应力-应变响应与宏观韧性指标图3a展示了三种玻璃和完美晶体在单轴拉伸下的应力-应变曲线。这是最宏观的性能对比晶体表现出典型的脆性断裂特征——线弹性变形直至应力峰值随后应力骤降断裂应变小。NCG玻璃曲线与晶体类似但峰值应力强度和断裂应变略有降低。断裂是突然的一旦裂纹从晶界萌生便迅速扩展贯穿样品。Para玻璃曲线出现了一定的非线性塑性特征断裂应变明显增加。应力在达到峰值后并非瞬间跌落而是有一个缓慢下降的过程。CRN玻璃表现出最高的断裂应变和最显著的“延性”。应力-应变曲线在峰值后有一个很长的下降尾巴意味着裂纹扩展过程是缓慢、受阻的。峰值应力代表材料的强度而应力-应变曲线下的面积则代表了材料断裂前吸收的能量是衡量断裂韧性的直观指标。显然CRN Para NCG ≈ Crystal。5.2 裂纹扩展的两种模式与Griffith理论修正我们追踪了裂纹尺寸归一化为样品尺寸随应变的变化图4a,b。发现裂纹扩展存在两种模式直接快速扩展主要发生在NCG和晶体中。一旦主裂纹形成便迅速扩展几乎不遇到阻力。受阻缓慢扩展伴随空洞化主要发生在CRN和部分Para结构中。主裂纹前端会引发大量次级空洞的形核图4c中void数量增加。这些空洞或与主裂纹汇合coalescence或通过剪切带与主裂纹桥接bridging如图4e所示。这个过程消耗了大量额外能量显著提高了韧性。经典的Griffith理论认为当裂纹尺寸达到临界值c_cr 4γ / (πEε^2)时裂纹将失稳扩展。其中γ是表面能E是杨氏模量ε是外加应变。我们从模拟中提取了γ和E计算了c_cr图4d中虚线。结果显示对于NCG观测到的扩展裂纹尺寸与c_cr吻合较好。但对于CRN许多裂纹在远小于c_cr时就开始缓慢扩展而大量形核的次级空洞尺寸则小于c_cr且不扩展非传播空洞。这说明在非晶材料中由于结构无序导致的应力场异质性图4f裂纹尖端并非处于理想的K场应力强度因子场局部波动很大。Griffith理论基于连续介质力学和均匀材料假设在原子尺度的无序结构面前需要修正。5.3 断裂路径粗糙度与“无序捕获”效应我们定量分析了断裂后表面的粗糙度图3b。粗糙度定义为断裂轮廓线相对于平均高度的均方根值。NCG断裂面最光滑因为裂纹倾向于沿着相对平直的晶界扩展。CRN断裂面最粗糙。裂纹在扩展过程中不断被局域的无序结构所偏折、分叉甚至暂时“被困住”disorder-trapping effect需要更大的驱动力才能继续前进如图4g所示。这种“之”字形的扩展路径极大地增加了新生表面积从而消耗了更多能量表面能γ本质上就是创造新表面所需的能量。Para粗糙度介于两者之间。“空洞引导”的裂纹路径是CRN高韧性的另一个关键机制图4g插图。裂纹并非总是直线前进而是倾向于朝着预先存在的或新形核的次级空洞方向扩展因为连接两个空洞所需的能量可能低于直接撕裂完整的网络。这种路径选择进一步增加了裂纹扩展的曲折度。6. 模拟实践中的常见问题与深度排查指南分子动力学模拟非晶材料断裂是一个复杂的过程从建模到分析处处是坑。以下是我在大量实践中总结出的典型问题与解决方案。6.1 模型构建与弛豫阶段问题问题现象可能原因排查方法与解决方案能量最小化不收敛或收敛后能量极高1. 初始原子构型存在严重重叠如原子间距小于0.5 Å。2. 势函数参数文件未正确载入或单位制不匹配。3. 边界条件设置错误。1. 检查初始模型生成脚本确保随机撒点或结构生长算法有最小距离限制。2. 使用可视化软件如OVITO、VMD检查初始构型手动删除或调整过于接近的原子。3. 仔细核对势函数文件路径、格式以及模拟程序中设定的质量、距离、能量单位是否与势函数参数一致。NPT弛豫后密度严重偏离实验值1. 控压算法参数如驰豫时间设置不当。2. 势函数本身对平衡体积的预测不准。3. 温度设置过高接近或超过Tg。1. 调整Parrinello-Rahman控压算法的质量参数或等效的阻尼参数使其与体系规模匹配。可以先在更小的体系上调试。2. 这是根本性问题。需要更换或修正势函数。可查阅文献看所用势函数对密度的预测能力如何。3. 对于玻璃弛豫温度应远低于Tg对于二氧化硅Tg约1200K-1500K模拟中常设在300K-500K。RDF与实验数据对不上特别是第二峰1. 淬火速率过快或过慢。2. 势函数对中程有序的描述能力不足。3. 体系尺寸太小统计性不足。1. 调整从熔体到玻璃的冷却速率。通常需要尝试不同的淬火速率如10K/ps, 1K/ps观察其对RDF第二峰形状的影响。2. 尝试不同的势函数。BKS势对短程有序描述好但对中程有序RDF第二峰的描述可能不如一些多体势。3. 增大模型尺寸。对于CRN至少需要数千个原子才能获得有统计意义的RDF。6.2 拉伸模拟与数据分析阶段问题问题现象可能原因排查方法与解决方案应力-应变曲线噪声过大1. 应变率过高系统远离准静态平衡。2. 温度控制算法引起热波动干扰应力信号。3. 应力输出频率太低未进行时间平均。1. 在计算资源允许的情况下尽可能降低应变率。同时确保在每次应变增量后给予足够的时间步数进行NVT弛豫通常50-100步。2. 使用Nosé-Hoover热浴时适当增加其弛豫时间常数以减少温度振荡。也可以尝试Berendsen热浴虽不严格但更平稳。3. 计算应力时不要只用单个时间步的瞬时值。应在每个应变状态下取一段弛豫时间内的应力平均值作为该应变下的应力。断裂行为不可重复同种模型每次模拟断裂应变差异巨大1. 体系尺寸太小个别原子的热涨落或局部缺陷对整体行为影响过大。2. 初始速度分布不同随机种子不同。3. 对于CRN等均匀结构断裂本身具有统计性。1.这是最重要的原因。必须进行有限尺寸效应分析。逐步增大体系尺寸如从1000原子到10000原子观察断裂应变的分布是否收敛。发表结论所基于的体系尺寸必须足够大。2. 这正是非晶材料断裂的内在特性。解决方案不是追求单次模拟的“正确”而是进行系综平均。对同一种结构用不同的随机种子生成多个独立样本至少20-30个进行拉伸模拟然后统计分析其断裂应变、形核位置的概率分布。3. 明确区分是“误差”还是“涨落”。通过大量样本的统计可以获得断裂应变的平均值和标准差这才是描述材料性能的可靠指标。无法清晰识别空洞形核的起始时刻和位置1. 轨迹文件输出频率不够高错过了形核的瞬间。2. 空洞识别算法阈值设置不合理。1. 在接近预估的断裂应变区间大幅提高轨迹输出频率如每10步输出一次。2. 空洞识别通常基于原子局部密度或配位数。例如可以计算每个原子周围一定截断半径内的邻居数。如果邻居数低于某个阈值如理想配位数的一半则认为该原子处于空洞内部。需要尝试不同的截断半径和阈值并与可视化结果对照找到能稳健识别初始微小空洞的参数。裂纹扩展路径分析困难无法自动提取粗糙度1. 断裂后原子位置混乱难以定义清晰的断裂面。2. 自定义分析脚本编写复杂。1. 一种有效方法是在断裂完成后对体系进行一个短时间的低温弛豫如50K下弛豫10ps让原子振动减小断裂面原子会相对聚集得更清晰。2. 利用强大的后处理工具。OVITO软件及其Python接口是神器。可以编写脚本自动识别断裂面上的原子例如通过坐标和局部原子密度然后对这些原子的坐标进行多项式拟合得到平均断面再计算各原子到该面的距离的均方根值RMS作为粗糙度。社区有很多现成的脚本可以参考和修改。6.3 性能优化与计算资源管理大规模分子动力学模拟极其消耗计算资源。一次包含数万原子、数百万时间步的模拟在普通工作站上可能需要数周时间。并行计算务必使用支持并行如MPI、OpenMP的分子动力学软件如LAMMPS。根据体系特点选择最优的并行分区策略。邻居列表更新合理设置邻居列表的截断半径和更新频率。对于断裂模拟原子相对运动剧烈需要比平衡模拟更频繁地更新邻居列表。增量式应变加载与其在每一步都重新标定所有原子坐标不如采用“变形-弛豫”的循环。每次施加一个有限应变增量如0.25%然后在NVT下弛豫足够步数再施加下一个增量。这更接近准静态过程且允许在弛豫阶段使用更大的时间步长。数据输出策略全轨迹输出所有原子的位置和速度会产生TB级的数据。只输出你分析所必需的数据。例如对于应力-应变曲线只需输出系统整体的应力和应变对于局部分析可以只在高频输出特定原子组的信息或者每隔很多步才输出一次全轨迹用于可视化。通过这套系统的建模、模拟、分析和排错流程我们才能从原子运动的噪声中提炼出决定非晶材料断裂行为的清晰物理图像。这项研究不仅加深了对玻璃这一古老材料本质的理解其揭示的“结构异质性调控韧性”的普适原理也为设计更高性能的金属玻璃、高分子玻璃等非晶材料提供了宝贵的微观设计思路。
分子动力学模拟揭秘：非晶材料断裂韧性的原子尺度起源
1. 项目概述从原子视角窥探玻璃的“筋骨”玻璃这种我们日常生活中再熟悉不过的材料其本质是一种非晶态固体。与金属或晶体不同玻璃内部的原子排列没有长程有序的周期性呈现出一种“冻结的液体”状态。这种无序性赋予了玻璃透明、各向同性等独特性质但也使其在力学性能上尤其是断裂行为充满了复杂性和不确定性。为什么有些玻璃一碰就碎而有些却能承受巨大的弯曲这个问题的答案深埋在原子尺度的微观世界里。传统的材料研究多依赖于宏观实验但面对非晶材料这种“黑箱”宏观数据往往难以揭示其内在的物理机制。分子动力学模拟就像一台超级显微镜允许我们追踪每一个原子的运动轨迹在计算机中“制备”出具有特定微观结构的材料模型并对其施加虚拟的载荷从而直观地观察从弹性变形到最终断裂的全过程。本次分享的项目正是运用这一利器聚焦于一个经典的模型体系——二维二氧化硅玻璃。选择二维体系并非为了简化而是为了剥离三维空间的复杂性更清晰地聚焦于非晶结构本身对力学行为的影响。我们构建了三种具有代表性的微观结构模型纳米晶玻璃、准晶玻璃和连续随机网络玻璃旨在系统性地回答一个核心问题原子尺度上的无序结构差异如何一步步决定材料的宏观断裂韧性对于材料科学、物理或相关领域的从业者无论是初涉模拟的研究生还是希望深化对非晶材料理解的技术人员这篇文章将带你深入一场从原子建模到宏观性能解读的完整旅程。我们将不仅复现模拟的关键步骤更会拆解每一个结果背后的物理图像分享参数设置背后的考量以及如何从海量的原子轨迹数据中提炼出有价值的科学见解。2. 核心模型构建与结构表征解析分子动力学模拟的第一步也是决定模拟成败与物理意义的关键一步就是构建能够真实反映材料特性的原子模型。对于二维二氧化硅其基本结构单元是硅氧四面体在二维平面上的投影即一个硅原子与三个氧原子键合形成的三角形单元这些单元通过共用氧原子连接成网络。2.1 三种非晶结构模型的生成策略我们构建了三种截然不同的二维二氧化硅玻璃模型它们代表了非晶结构中从“有序”到“无序”的谱系。2.1.1 纳米晶玻璃模型顾名思义NCG模型由许多微小的晶粒镶嵌在非晶基质中构成。我们的生成方法是首先在模拟盒子中随机撒播多个晶核点每个晶核以二维β-鳞石英的晶体结构一种二氧化硅的晶体相为模板进行生长。通过控制晶核的数量和生长速率我们可以获得不同晶粒尺寸和体积分数的纳米晶结构。当晶粒生长到相互接触或达到预设尺寸后剩余的空间则通过快速淬火即高温熔体快速冷却的方式形成非晶基质。这个过程模拟了某些通过控制结晶工艺制备的玻璃陶瓷材料。注意晶粒与基质的界面是关键。在模拟中需要确保界面处的原子通过能量最小化进行了充分弛豫以避免引入不切实际的高能界面从而影响断裂起始位置的判断。2.1.2 准晶玻璃模型Para模型比NCG更加无序但保留了中程有序的特征。我们采用了一种“拓扑约束”生长算法。首先生成一个完全随机的初始原子构型然后通过引入一组约束条件如键长、键角的目标分布函数这些函数来源于对实验或高质量CRN结构的统计进行迭代优化。在优化过程中系统被允许形成一些局域的、类晶的短程有序团簇但这些团簇之间没有明确的晶界呈现出一种“近程有序、中程弥散”的状态。这类似于某些通过特殊工艺制备的、结构高度均匀的非晶合金。2.1.3 连续随机网络模型CRN模型是描述理想玻璃态最经典的模型。它假设网络是连续、无限且完全无序的但每个原子都满足确定的配位数硅为4氧为2。我们采用Wooten-Winer-Weaire算法或其变体来生成。该算法从一个晶态种子开始通过一系列随机的键切换操作来引入无序同时严格保持网络的连通性和配位数约束直到系统的能量和径向分布函数收敛到典型玻璃态的特征。最终得到的结构没有任何周期性是三种模型中最无序的。2.2 关键结构表征量的计算与物理意义模型建好后如何定量描述它们的结构差异这就需要借助一系列结构表征函数。2.2.1 径向分布函数探测短程有序的“尺子”RDF是材料科学中最基础的结构表征工具。对于我们的体系我们计算了Si-Si、Si-O和O-O的RDF。计算原理是以某个原子为中心统计在距离r到rdr的球壳内发现其他原子的平均数量然后对体系中的所有原子和所有构型进行平均并归一化。第一峰位置对应最近邻原子间的平均距离。例如Si-O RDF的第一峰位置直接给出了Si-O键的平均键长这是判断势函数合理性和结构弛豫程度的重要指标。第一峰面积积分后可以得到配位数。Si-O第一峰面积应接近4Si的配位数O-Si第一峰面积应接近2O的配位数。峰的宽度和衰减峰越尖锐说明键长分布越集中有序度越高。从第二、第三峰开始对于晶体RDF会出现清晰的峰对于CRN玻璃这些峰会迅速展宽并湮灭这正是长程无序的体现。将模拟CRN的RDF与实验测得的RDF如图1e,f进行对比是验证模型可靠性的黄金标准。2.2.2 环分布统计揭示中程有序的“拓扑地图”在二维网络中环由Si和O原子交替连接形成的闭合多边形是描述拓扑结构的核心。我们统计了网络中所有最小环的边数通常是4-10元环。对于完美的β-鳞石英晶体只存在6元环。而在非晶网络中会出现4、5、7、8元环等多种拓扑缺陷。环分布图如图1d所示NCG模型由于含有晶粒其环分布会在6元环处出现一个尖锐的峰同时在其他位置有非晶基质贡献的宽分布。Para模型的分布较宽但可能在6元环附近仍有凸起。CRN模型的分布最宽且平缓是典型的无特征分布。环分布的差异直接影响了网络的刚性、自由体积分布进而影响力学性能。2.2.3 角度分布函数与结构相似性参数ADF统计了诸如Si-O-Si和O-Si-O等键角的分布。对于二氧化硅Si-O-Si键角分布通常在120°-180°之间中心约在145°。ADF的宽度反映了网络键角扭曲的程度。结构相似性参数s则是一个更全局的量度用于量化局部原子环境与某个参考结构如理想的晶体环境的差异其概率分布函数图1g可以直观显示体系内结构异质性的程度。通过这些细致的结构表征我们不仅在视觉上图1a-c更在定量上严格区分了三种模型为后续关联结构与性能奠定了坚实的基础。模拟得到的CRN的RDF和ADF与实验数据高度吻合图1e,f,h,i这让我们对模拟的可靠性充满信心。3. 分子动力学模拟流程与关键参数设置有了可靠的原子模型下一步就是让它们“动”起来接受力学测试。分子动力学模拟的本质是数值求解牛顿运动方程其流程和参数设置直接决定了模拟的物理真实性和计算效率。3.1 势函数选择描述原子间相互作用的“宪法”对于二氧化硅体系原子间的相互作用通常由经验势函数来描述。我们选择了基于BKS模型的修正势函数。BKS势是一种经典的二体势包含了长程库仑相互作用和短程Buckingham势的形式。E_ij q_i q_j / r_ij A_ij exp(-B_ij r_ij) - C_ij / r_ij^6其中第一项是库仑项q_i为原子电荷第二项和第三项是短程排斥和吸引项。我们采用Wolf方法对长程库仑力进行截断和修正以在保证精度的前提下大幅提升计算速度。实操心得势函数的选择是分子动力学模拟的基石。对于二氧化硅除了BKS还有TTAM、CHIK等势函数。选择时需权衡精度与速度。BKS势在描述玻璃态结构如RDF、密度方面表现良好且计算效率较高适合进行大规模断裂模拟。务必在正式模拟前用小体系测试所选势函数能否复现实验已知的键长、键角、弹性模量等基本性质。3.2 系统弛豫与平衡态获取在施加载荷前必须让初始构建的模型达到平衡态即体系能量和压力等宏观量在平均值附近小幅波动。能量最小化首先对初始模型进行共轭梯度法能量最小化消除原子间的过度重叠和不合理的高应力。NPT系综弛豫在恒温恒压系综下运行一段时间通常数十皮秒让体系密度弛豫到平衡值。温度通过Nosé-Hoover热浴控制压力通过Parrinello-Rahman控压方法控制。对于二维体系需注意控制面内压力而垂直于平面方向第三维的压力设为0。NVT系综弛豫在平衡体积下切换至恒温恒容系综继续弛豫确保温度均匀体系处于稳定的玻璃态。我们通常将温度设定在远低于玻璃化转变温度Tg的数值如300K。关键参数弛豫时间需要足够长。我们通过监测体系总能量、温度、压力以及RDF是否不再有漂移来判断平衡。对于含有数千原子的体系NPTNVT弛豫总计需要至少100-200皮秒。3.3 单轴拉伸模拟的实现细节我们采用位移控制的方式进行单轴拉伸模拟这是研究断裂行为的常用方法。加载方向沿模拟盒子的x方向armchair方向施加应变。同时允许y方向自由收缩泊松效应而z方向二维平面外保持固定厚度。应变率这是一个至关重要的参数。分子动力学模拟受限于计算资源使用的应变率通常为10^8 - 10^9 /s远高于实验速率10^-3 - 10^0 /s。如此高的应变率会抑制原子的热激活运动可能导致材料表现得比实际更脆。为了部分抵消这一效应我们采用了相对“慢”的应变率在MD尺度下例如5×10^8 /s并通过与不同应变率的对比测试评估应变率效应的影响。模拟步骤在每个模拟时间步通常为1飞秒将模拟盒子在x方向的长度增加一个微小量ΔL ε˙ * L_x0 * Δt其中ε˙是应变率L_x0是初始长度Δt是时间步长。然后重新标定所有原子的x坐标仿射变形接着在NVT系综下运行若干步让原子弛豫到新的构型并计算维里应力。应力计算体系在变形过程中的应力张量通过维里定理计算包含了动能项和原子间相互作用力的贡献。我们主要关注沿拉伸方向的轴向应力σ_xx。关键技巧为了捕捉断裂起始的精确位置我们在整个拉伸过程中以很高的频率如每100步输出原子的位置、速度和应力数据。这些海量的轨迹文件是后续所有分析的源头。4. 断裂起始机制从最弱环节到空洞形核材料不会无缘无故断裂断裂总是从最薄弱的环节开始。通过分子动力学模拟我们可以像侦探一样在断裂发生前就锁定这些“嫌疑位置”并观察“犯罪现场”是如何形成的。4.1 维里系数预测断裂起始点的微观应力探针宏观应力是体系整体的平均值而断裂往往由局部极高的应力集中引发。维里系数V_θ是我们定义的一个局部应力度量。对于每个原子i我们计算其在一个局部区域例如其最近邻壳层内的应力张量然后通过坐标变换得到沿拉伸方向θ的正应力分量V_θ(i)。这个值直观反映了该原子所在局部区域被拉伸的紧张程度。如图2a所示在未变形的平衡结构中我们就可以绘制出V_θ的空间分布图。那些呈现高正值红色区域的原子团簇意味着即使在零应变下该处的局部网络也处于一种“预拉伸”的亚稳态。它们就是潜在的断裂起始点。我们的模拟结果惊人地显示最终发生空洞形核的位置图2a中圆圈标注与这些高V_θ区域高度重合。4.2 空洞形核的临界条件与结构关联我们统计了所有模拟中空洞形核时的应变值ε_nuc。分析发现空洞形核应变与平衡结构中的两个关键结构参数存在强烈的相关性与最拉伸状态的关系对于每个样本我们找到其平衡结构中V_θ的最小值即最拉伸的状态。如图2b所示V_θ_min与ε_nuc之间存在近似指数衰减的关系。V_θ_min越小即初始就越拉伸材料越早发生空洞形核韧性越差。这为预测材料的断裂起始韧性提供了一个基于初始结构的微观判据。与Si-O-Si键长的关系我们测量了所有Si-O-Si桥氧键的长度即两个硅原子之间的距离。统计发现空洞往往在那些具有较长Si-O-Si键的附近形核图2d。较长的键通常意味着键能较弱或者该处的网络拓扑结构较为松散自由体积较大在拉伸下更容易被拉开。注意事项空洞形核是一个统计过程。即使在同一模型中由于原子级结构的随机涨落每次模拟的形核位置和应变也可能不同。因此我们需要对每种结构进行大量如50-100次的独立模拟并统计其累积概率分布函数CDF如图2c所示。CRN的CDF最宽说明其形核应变的波动性最大这与它结构的均匀性没有明显的优先形核点有关而NCG的CDF较窄因为晶界作为明确的缺陷总是优先形核。4.3 断裂事件的原子尺度快照通过分析断裂瞬间的原子构型图3d, g, j我们可以清晰地看到空洞形核的微观过程在NCG中空洞几乎无一例外地在晶界处形核。晶界处原子排列不规则键合较弱是天然的弱点。在Para中空洞在那些中程有序被破坏的局部松散区域形核。在CRN中空洞形核位置看起来更随机但追溯其平衡结构总能找到对应的V_θ高点或长键位置。形核后空洞并非立即导致灾难性断裂。它首先会稳定生长直到达到一个临界尺寸。5. 裂纹扩展动力学与韧性差异分析空洞形核只是断裂的序幕随后的裂纹扩展阶段才真正决定了材料的断裂韧性和断裂形貌。这是三种结构模型表现出最大差异的环节。5.1 应力-应变响应与宏观韧性指标图3a展示了三种玻璃和完美晶体在单轴拉伸下的应力-应变曲线。这是最宏观的性能对比晶体表现出典型的脆性断裂特征——线弹性变形直至应力峰值随后应力骤降断裂应变小。NCG玻璃曲线与晶体类似但峰值应力强度和断裂应变略有降低。断裂是突然的一旦裂纹从晶界萌生便迅速扩展贯穿样品。Para玻璃曲线出现了一定的非线性塑性特征断裂应变明显增加。应力在达到峰值后并非瞬间跌落而是有一个缓慢下降的过程。CRN玻璃表现出最高的断裂应变和最显著的“延性”。应力-应变曲线在峰值后有一个很长的下降尾巴意味着裂纹扩展过程是缓慢、受阻的。峰值应力代表材料的强度而应力-应变曲线下的面积则代表了材料断裂前吸收的能量是衡量断裂韧性的直观指标。显然CRN Para NCG ≈ Crystal。5.2 裂纹扩展的两种模式与Griffith理论修正我们追踪了裂纹尺寸归一化为样品尺寸随应变的变化图4a,b。发现裂纹扩展存在两种模式直接快速扩展主要发生在NCG和晶体中。一旦主裂纹形成便迅速扩展几乎不遇到阻力。受阻缓慢扩展伴随空洞化主要发生在CRN和部分Para结构中。主裂纹前端会引发大量次级空洞的形核图4c中void数量增加。这些空洞或与主裂纹汇合coalescence或通过剪切带与主裂纹桥接bridging如图4e所示。这个过程消耗了大量额外能量显著提高了韧性。经典的Griffith理论认为当裂纹尺寸达到临界值c_cr 4γ / (πEε^2)时裂纹将失稳扩展。其中γ是表面能E是杨氏模量ε是外加应变。我们从模拟中提取了γ和E计算了c_cr图4d中虚线。结果显示对于NCG观测到的扩展裂纹尺寸与c_cr吻合较好。但对于CRN许多裂纹在远小于c_cr时就开始缓慢扩展而大量形核的次级空洞尺寸则小于c_cr且不扩展非传播空洞。这说明在非晶材料中由于结构无序导致的应力场异质性图4f裂纹尖端并非处于理想的K场应力强度因子场局部波动很大。Griffith理论基于连续介质力学和均匀材料假设在原子尺度的无序结构面前需要修正。5.3 断裂路径粗糙度与“无序捕获”效应我们定量分析了断裂后表面的粗糙度图3b。粗糙度定义为断裂轮廓线相对于平均高度的均方根值。NCG断裂面最光滑因为裂纹倾向于沿着相对平直的晶界扩展。CRN断裂面最粗糙。裂纹在扩展过程中不断被局域的无序结构所偏折、分叉甚至暂时“被困住”disorder-trapping effect需要更大的驱动力才能继续前进如图4g所示。这种“之”字形的扩展路径极大地增加了新生表面积从而消耗了更多能量表面能γ本质上就是创造新表面所需的能量。Para粗糙度介于两者之间。“空洞引导”的裂纹路径是CRN高韧性的另一个关键机制图4g插图。裂纹并非总是直线前进而是倾向于朝着预先存在的或新形核的次级空洞方向扩展因为连接两个空洞所需的能量可能低于直接撕裂完整的网络。这种路径选择进一步增加了裂纹扩展的曲折度。6. 模拟实践中的常见问题与深度排查指南分子动力学模拟非晶材料断裂是一个复杂的过程从建模到分析处处是坑。以下是我在大量实践中总结出的典型问题与解决方案。6.1 模型构建与弛豫阶段问题问题现象可能原因排查方法与解决方案能量最小化不收敛或收敛后能量极高1. 初始原子构型存在严重重叠如原子间距小于0.5 Å。2. 势函数参数文件未正确载入或单位制不匹配。3. 边界条件设置错误。1. 检查初始模型生成脚本确保随机撒点或结构生长算法有最小距离限制。2. 使用可视化软件如OVITO、VMD检查初始构型手动删除或调整过于接近的原子。3. 仔细核对势函数文件路径、格式以及模拟程序中设定的质量、距离、能量单位是否与势函数参数一致。NPT弛豫后密度严重偏离实验值1. 控压算法参数如驰豫时间设置不当。2. 势函数本身对平衡体积的预测不准。3. 温度设置过高接近或超过Tg。1. 调整Parrinello-Rahman控压算法的质量参数或等效的阻尼参数使其与体系规模匹配。可以先在更小的体系上调试。2. 这是根本性问题。需要更换或修正势函数。可查阅文献看所用势函数对密度的预测能力如何。3. 对于玻璃弛豫温度应远低于Tg对于二氧化硅Tg约1200K-1500K模拟中常设在300K-500K。RDF与实验数据对不上特别是第二峰1. 淬火速率过快或过慢。2. 势函数对中程有序的描述能力不足。3. 体系尺寸太小统计性不足。1. 调整从熔体到玻璃的冷却速率。通常需要尝试不同的淬火速率如10K/ps, 1K/ps观察其对RDF第二峰形状的影响。2. 尝试不同的势函数。BKS势对短程有序描述好但对中程有序RDF第二峰的描述可能不如一些多体势。3. 增大模型尺寸。对于CRN至少需要数千个原子才能获得有统计意义的RDF。6.2 拉伸模拟与数据分析阶段问题问题现象可能原因排查方法与解决方案应力-应变曲线噪声过大1. 应变率过高系统远离准静态平衡。2. 温度控制算法引起热波动干扰应力信号。3. 应力输出频率太低未进行时间平均。1. 在计算资源允许的情况下尽可能降低应变率。同时确保在每次应变增量后给予足够的时间步数进行NVT弛豫通常50-100步。2. 使用Nosé-Hoover热浴时适当增加其弛豫时间常数以减少温度振荡。也可以尝试Berendsen热浴虽不严格但更平稳。3. 计算应力时不要只用单个时间步的瞬时值。应在每个应变状态下取一段弛豫时间内的应力平均值作为该应变下的应力。断裂行为不可重复同种模型每次模拟断裂应变差异巨大1. 体系尺寸太小个别原子的热涨落或局部缺陷对整体行为影响过大。2. 初始速度分布不同随机种子不同。3. 对于CRN等均匀结构断裂本身具有统计性。1.这是最重要的原因。必须进行有限尺寸效应分析。逐步增大体系尺寸如从1000原子到10000原子观察断裂应变的分布是否收敛。发表结论所基于的体系尺寸必须足够大。2. 这正是非晶材料断裂的内在特性。解决方案不是追求单次模拟的“正确”而是进行系综平均。对同一种结构用不同的随机种子生成多个独立样本至少20-30个进行拉伸模拟然后统计分析其断裂应变、形核位置的概率分布。3. 明确区分是“误差”还是“涨落”。通过大量样本的统计可以获得断裂应变的平均值和标准差这才是描述材料性能的可靠指标。无法清晰识别空洞形核的起始时刻和位置1. 轨迹文件输出频率不够高错过了形核的瞬间。2. 空洞识别算法阈值设置不合理。1. 在接近预估的断裂应变区间大幅提高轨迹输出频率如每10步输出一次。2. 空洞识别通常基于原子局部密度或配位数。例如可以计算每个原子周围一定截断半径内的邻居数。如果邻居数低于某个阈值如理想配位数的一半则认为该原子处于空洞内部。需要尝试不同的截断半径和阈值并与可视化结果对照找到能稳健识别初始微小空洞的参数。裂纹扩展路径分析困难无法自动提取粗糙度1. 断裂后原子位置混乱难以定义清晰的断裂面。2. 自定义分析脚本编写复杂。1. 一种有效方法是在断裂完成后对体系进行一个短时间的低温弛豫如50K下弛豫10ps让原子振动减小断裂面原子会相对聚集得更清晰。2. 利用强大的后处理工具。OVITO软件及其Python接口是神器。可以编写脚本自动识别断裂面上的原子例如通过坐标和局部原子密度然后对这些原子的坐标进行多项式拟合得到平均断面再计算各原子到该面的距离的均方根值RMS作为粗糙度。社区有很多现成的脚本可以参考和修改。6.3 性能优化与计算资源管理大规模分子动力学模拟极其消耗计算资源。一次包含数万原子、数百万时间步的模拟在普通工作站上可能需要数周时间。并行计算务必使用支持并行如MPI、OpenMP的分子动力学软件如LAMMPS。根据体系特点选择最优的并行分区策略。邻居列表更新合理设置邻居列表的截断半径和更新频率。对于断裂模拟原子相对运动剧烈需要比平衡模拟更频繁地更新邻居列表。增量式应变加载与其在每一步都重新标定所有原子坐标不如采用“变形-弛豫”的循环。每次施加一个有限应变增量如0.25%然后在NVT下弛豫足够步数再施加下一个增量。这更接近准静态过程且允许在弛豫阶段使用更大的时间步长。数据输出策略全轨迹输出所有原子的位置和速度会产生TB级的数据。只输出你分析所必需的数据。例如对于应力-应变曲线只需输出系统整体的应力和应变对于局部分析可以只在高频输出特定原子组的信息或者每隔很多步才输出一次全轨迹用于可视化。通过这套系统的建模、模拟、分析和排错流程我们才能从原子运动的噪声中提炼出决定非晶材料断裂行为的清晰物理图像。这项研究不仅加深了对玻璃这一古老材料本质的理解其揭示的“结构异质性调控韧性”的普适原理也为设计更高性能的金属玻璃、高分子玻璃等非晶材料提供了宝贵的微观设计思路。
