二维随机蜂巢网格熔断动力学过程和熔断面标度性质的数值模拟*

李瑞涛 唐刚 夏辉 寻之朋 李嘉翔 朱磊

(中国矿业大学物理科学与技术学院, 徐州 221116)

石墨烯等材料具有典型的二维蜂巢结构,而随机电阻丝模型则是研究非均匀材料断裂十分有效的统计物理学模型.本文尝试对二维蜂巢结构随机电阻丝网络熔断的动力学过程及熔断面性质进行数值模拟分析,以此来研究二维非均质蜂窝材料熔断的动力学性质和熔断面的动力学标度性质. 模拟研究表明, 二维随机蜂窝网格的熔断动力学过程和熔断面具有明显的标度性质, 得到的熔断面整体和局域粗糙度指数分别为和 , 这两者之间的明显差异表明熔断面具有奇异标度性. 通过对熔断面极值高度的分析发现, 熔断面高度的极值统计分布能很好地满足Asym2sig型分布, 而不是最常见的三种极值统计分布. 本文的研究表明, 随机电阻丝模型在模拟非均匀材料的电流熔断过程和熔断表面标度性的分析中同样适用和有效.

1 引 言

二维蜂巢结构是覆盖二维平面的最佳拓扑结构, 其构成是由一个个正六边形单房、房口全部朝下或朝向一边、背对背对称排列组合而成. 同时二维蜂巢结构也是一种十分重要的材料结构形式, 这其中就包括石墨烯等重要材料的结构. 石墨烯是由碳原子以杂化轨道组成六角型呈蜂巢状晶格的平面薄膜, 是一种可以只有一层原子厚度的二维材料. 由于其高导电率以及独特的电子特性, 被认为是下一代电子材料中最有前途的候选材料之一, 并具有广泛的应用前景[1]. 因而研究其导电性能和电流熔断机理以及熔断面的标度属性具有重要的理论和实践意义.

自从2004年二维石墨烯薄膜材料被发现以来, 人们对二维石墨烯材料的结构以及导电性能都进行了广泛和深入的研究, 并取得了很多重要的理论和实验成果[2]. 比如, Lu等[3]通过理论模型分析,发现外部电场对研究石墨烯系统的原子和电子结构都有重要影响, 指出控制石墨烯电子结构的重要性, 并且发现石墨烯上的吸附原子可以作为调节电子性质的有效工具等. 2013 年, Mour和 Marder[4]通过分子动力学模拟的方法对石墨烯的断裂力学进行了深入的研究, 建立了断裂的几何模型, 得到了临界裂纹长度和应力的近似表达式, 提出了改善石墨烯韧性的方法. 他们还发现, 裂纹的路径和产生的边缘结构依赖于初始裂纹的长度. 2016年,Ghorbanfekr-Kalashami等[5]通过反应力场(reaction force field)的方法对石墨烯的结构和力学性能进行了研究, 结果发现掺杂物的波纹改变了石墨烯断裂表面的粗糙度, 并且与掺杂物的数量和局部排列也有关.

在材料断裂方面, 非均匀材料的断裂机理和断裂规律、断裂表面的标度性质等近年来一直都是活跃的实验和理论研究领域[6]. 实验发现, 在不同载荷下的几种材料, 例如木材、蜂窝玻璃、混凝土和纸张[7-10]等, 其断裂表面具有分形结构和标度性质, 并具有普适性质. 材料内部结构的无序性、非均匀性和断裂过程的非线性决定了非均匀材料的断裂机理、断裂面形貌及其标度性质. 在非均匀材料断裂的微观机理和动力学过程的理论研究中, 通常是基于晶格模型的解析近似和数值模拟. 在对非均匀材料断裂动力学过程的数值模拟方面, 随机电阻丝模型则是最广泛使用和十分有效的方法[6].

随机电阻丝模型是Arcangelis等[11]在1985年引入的. 在随机电阻丝模型中, 是将材料看成由电阻丝组成的网格, 在网格两端加上电压, 用电阻丝中的电流强度来表示材料的应力, 用网格中各个电阻丝的断裂电流阈值的随机分布来模拟材料的非均匀性, 用电阻丝的熔断过程来模拟实际材料的断裂过程. 研究表明, 随机电阻丝模型可以较准确地模拟实际材料的断裂过程[12], 并能够得到断裂过程中的基本特征, 而且这样的模型相对简单和容易处理[13].

在随机电阻丝模型中有两个基本假设. 一是假设模型中电阻丝具有不可逆的熔断性质. 要求电阻丝网络模型满足连续Laplace方程

变阈值范围为

近年来, 在对随机电阻丝模型的研究中, 科学家们做了许多十分有价值的研究工作. Duxbury等[14]发现在淬火随机介质中电击穿的尺寸效应,即局域断裂理论, 证明了随机淬火介质中的有限缺陷部分可以定性地降低实际材料的电流击穿性能.Nukala等[15]通过数值模拟的方法分析了强无序随机电阻丝模型的损伤成核化和局域化, 说明了断裂的过程和特点, 找到了损伤标度律并指出损伤在大尺寸上的不相关性. 两年后他们又分析了三维随机电阻丝模型断裂粗糙度和断裂面的标度特征, 指出损伤累积是以发散的方式进行的, 直至峰值载荷处, 然后发生局域化. 同时他们还发现. 整体表面粗糙度指数与峰值载荷后损伤轮廓的局域化长度是一致的. 对不同系统尺寸的数据进行分析, 发现断裂宽度分布可以很好地塌缩在一起. Toussaint和Hansen[16]还从平均场理论对柱形电阻丝网络的断裂机制进行了研究, 通过分离和分析系统的相图, 找到了系统尺寸和损伤发生的特征尺寸之间的标度律. Jan Øystein 和 Hansen[17]在研究断裂面粗糙度的映射工作中指出, 在电阻丝模型中粗糙度指数是普遍的. 当晶格影响断裂生长时, 电阻丝模型的粗糙度指数会随着阈值分布而改变; 当影响消失时, 局部粗糙度指数趋于. 这些研究工作充分地证明了随机电阻丝模型在模拟非均匀材料断裂的微观机理和动力学过程中的适用性和有效性[6].

在此前的模拟研究中, 因为结构简单、计算方便而且计算量比较小, 所以采用的网络结构大多是三角网格和菱形网络[13], 尚未见对二维蜂巢结构直接进行模拟分析计算的工作. 二维蜂巢结构是最常见和最重要的网格结构之一, 也是石墨烯等重要材料所具有的晶格结构. 除了石墨烯以外, 如果用电阻丝断裂的阈值代表化学键的强度, 则随机电阻丝模型还可以用来研究其他的晶态材料, 对很多金属材料的强度的模拟分析研究是十分有帮助的. 因此, 对二维蜂巢结构的随机电阻丝网络的微观熔断机理和动力学过程进行模拟研究具有一定的理论和实践意义. 此外, 使用随机电阻丝模型对断裂微观机理和动力学过程分析模拟的比较多. 本文则尝试使用随机电阻丝模型来分析断裂面形貌及其标度性质.

本文对二维蜂巢结构随机电阻丝网络的熔断过程进行数值模拟分析, 目的是研究二维蜂巢结构随机电阻丝网络断裂的微观机理和动力学过程以及熔断面的形貌和标度性质. 通过对熔断面表面宽度和局域表面宽度的计算, 发现二维蜂巢结构随机电阻丝网络的熔断面呈现出标度性质并具有奇异标度性质, 得到的标度指数分别为和. 通过对二维蜂巢结构随机电阻丝网络熔断面极值高度的分析, 发现其极值高度的极大(小)值呈现一定的分布规律, 分布函数则能很好地满足Asym2sig型分布, 且同一系统尺寸下熔断面的相对极大和极小高度分布具有较好的对称性. 在模拟计算过程中, 使用节点分析法构建了系数矩阵, 并对系数矩阵进行Cholesky分解[18], 并采用 Sherman-Morrison-Woodbury 算法快速对系数矩阵求逆. 通过这些加速算法和对其结构的优化大大地提高了计算效率, 使得本文的数值计算和分析工作能够顺利进行.

2 二维蜂巢随机电阻丝网格及其熔断过程的数值模拟计算方法

石墨烯中各个原子的化学性质是相同的, 但是几何环境并不完全相同, 所以在一个基元中有两种不同环境的原子, 构成石墨烯的复式晶格结构. 对于随机电阻丝模型, 很多文献选择的是三角形、菱形[13]等. 本文则尝试使用随机电阻丝模型来研究二维蜂巢网格结构的熔断动力学过程和熔断面的动力学标度性质.

图1 石墨烯蜂巢结构随机电阻丝网络通电熔断示意图Fig. 1. Schematic diagram of random fuse model electric fuse in graphene honeycomb structure.

然后根据周期性边界条件得到一个基尔霍夫电压方程:

图2 2 × 2 的正方格子电流流向示意图Fig. 2. 2 × 2 square lattice current flow diagram.

最后通过初始电压再得到一个基尔霍夫电压方程:

联立以上6个方程得到一个方程组:

然后继续解方程, 断裂, 直到所有电阻丝上的电流都为0. 这样, 就可以得到电阻丝网的断裂顺序了.而对于解随机电阻丝模型的节点分析法则是, 对于同样的网络和电流方向, 先给出节点和电流的关联矩阵:

当电阻丝网断裂的时候, 只要改变关联矩阵中的数据就可以计算. 例如I1断裂，则关联矩阵变为

Sherman-Morrison-Woodbury算法是一种快速求逆的方法. 对于这种方法的做法是设

模拟过程中, 使用节点分析法构建系数矩阵,该系数矩阵为稀疏矩阵. 然后再对稀疏矩阵进行Cholesky分解. Cholesky分解是一种求解大型线性方程组的一种常见方法, 通过将对称矩阵分解成一个上三角矩阵和下三角矩阵从而加速线性方程组求解速度的方法, 即. 对于上下两个三角矩阵采用Sherman-Morrison-Woodbury算法快速求逆, 通过一些加速算法和对其结构的优化可以大幅度加快分解速度和运算速度, 从而提高计算机的模拟效率. 在本文的模拟分析计算中, 通过使用以上方法大大地提高了计算的效率, 使得本文的模拟计算工作得以顺利地进行.

3 二维蜂巢随机电阻丝网络熔断面粗糙度及极值高度的计算分析

3.1 二维蜂巢随机电阻丝网络熔断面粗糙度的计算分析

材料断裂现象是一种十分复杂的随机过程, 但其断裂面通常都能够形成自仿射的分形结构. 其粗糙度通常用表面宽度进行描述, 其定义为

图3 整体表面宽度 随系统尺寸 的对数-对数曲线Fig. 3. The log-logarithmic curve of the global surface width W with the system size L.

图4 局域表面宽度 随局域尺寸 l 的对数-对数曲线Fig. 4. The Log-logarithmic curve of local surface width w with local size l.

从表1可以看出, 二维菱形结构随机电阻丝模型断裂面的整体与局域粗糙度指数分别为和, 二维三角形结构随机电阻丝模型断裂面的整体与局域粗糙度指数分别为和0.003. 这两种模型结构是各向同性的, 整体与局域粗糙度的结果说明断裂面不存在奇异标度行为. 本文所研究的模型是各向异性的, 说明模型熔断表面存在着奇异标度行为.

表1 二维菱形、三角形及石墨烯蜂巢结构电阻丝网络熔断面整体与局域的粗糙度指数Table 1. Roughness index of the global and local of the burnout surface of two-dimensional diamond,triangle and graphene honeycomb structures.

3.2 二维蜂巢随机电阻丝网络熔断面极值高度的分析

除熔断面粗糙度之外, 本文还分析了熔断面相对高度的极值分布行为. 在自然界中, 极值事件对非平衡系统有着非常重要的作用[21]. 很多复杂的物理系统大多受到极值的影响[22]. 理论研究揭示出极值事件在描述非平衡系统方面具有重要意义. 对于非平衡系统, 极值统计早期用来分析生长表面的极值高度分布情况, 并且取得了很好的研究成果. 材料断裂过程是一种非平衡动力学过程, 研究其极值分布具有重要的意义[23,24]. 截止目前, 文献中很少有通过极值统计对其进行研究. 在表面界面生长领域, 研究发现, 对于基底尺寸为的有限系统, 相对高度的最大值在饱和区域的分布满足[25,26]

20 世纪 30 年代, Fisher和 Tippett[27]在对独立同分布的极大(小)值渐进分布进行理论研究时提出了三种极值分布: Ⅰ型的Gumbel分布, Ⅱ型的 Frechet分布和Ⅲ型的 Weibull分布. 广义Gumbel分布函数作为描述不同分布下样本容量中极大(小)的分布, 有着非常重要的作用, 如文献[28—30]指出, 关联物理系统中可用Gumbel分布函数来描述全局涨落. 但有研究结果表明Gumbel分布函数不能很好地描述非平衡饱和表面的极值统计分布行为. Oliveira等[31]在研究Kardar-Parisi-Zhang (KPZ)和 Villain-Lai-Das Sarma(VLDS) 普适类的表面界面生长模型时发现, 具有不对称局域高度分布的表面界面生长模型的极大和极小值分布不同. Wen等[32]在研究 1+1维Wolf-Villain模型饱和表面的生长高度极值统计分布时发现相对极小值分布不满足广义Gumbel分布. 而研究发现, 在很多的极值统计中, 也常常满足的是Asym2sig型函数分布. 如Cui等[33]研究HSSF-CW关于示踪物的停滞时间密度分布和Brar[34]研究光致发光谱的强度的极值分布, 在我们前面工作中[35,36], 也发现是满足Asym2sig型分布. 本文研究发现二维蜂巢随机电阻丝网络熔断面极值高度具有较好的统计行为, 极值高度的分布较好地符合Asym2sig型函数分布, 其表达式为

图5 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面相对极大高度分布Fig. 5. Relative maximum height distribution of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

图6 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面相对极小高度分布Fig. 6. Relative minimum height distribution of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

图5 为不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面相对极大高度的概率分布图.横坐标定义为,表示相对极大高度的统计平均值,表示相对极大高度的标准差; 纵坐标定义为相对极大高度与峰值处极值高度的对应统计次数的比值,范围在 [0, 1]之间. 其中, 离散的点为数值模拟结果, 实线为Asym2sig函数拟合曲线, 相关参数的拟合值如表2所列, 影响峰值宽度参数的在三种基底尺寸下拟合所得数值在误差范围内可认为是相等的, 说明不同的基底尺寸不影响熔断面极值高度的分布. 图6为不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面相对极小高度的概率分布图.

表2 系 统 尺 寸 为 L = 384, 512, 768 时Asym2sig函数拟合的参数Table 2. Parameters of Asym2sig function fitting when the system size is L = 384, 512, 768.

为了进一步说明熔断面相对极值满足的标度规律, 本文对纵坐标做半对数处理后发现: 在不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面的相对极大(小)高度分布依然呈现出较好的标度规律, 如图7和图8所示.

很显然, 图7和图8表明在不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面的相对极大(小)高度分布满足一定的标度规律, 很好地服从Asym2sig峰值分布函数.

图7 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面的相对极大高度的半对数分布Fig. 7. Semi-logarithmic distribution of the relative maximum height of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

图8 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面的相对极小高度的半对数分布Fig. 8. Semi-logarithmic distribution of the relative minimum height of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

本文还对同一尺寸下熔断面极值高度的极大值和极小值进行了比较, 如图9和图10所示. 结果显示, 同一系统尺寸下熔断面的相对极大(小)高度分布能够很好地重合在一起, 表明熔断表面极大(小)高度分布具有对称性.

图9 系统尺寸 L = 384 的熔断面的相对极大 (小)高度分布 Fig. 9. Relatively maximum (minimum) height distribution of fracture surface with system size L = 384.

图10 系统尺寸 L = 512 的熔断面的相对极大 (小)高度分布Fig. 10. Relatively maximum (minimum) height distribution of fracture surface with system size L = 512.

4 结 论

本文构建了六边形蜂巢结构的随机电阻丝网络模型, 通过对电阻丝网络施加电压使其断裂, 发现该模型断裂的机理以及断裂面的一些标度性质,从理论上丰富了随机电阻丝模型的应用领域, 同时具有一定的实践意义, 因为蜂窝结构广泛应用于材料力学、电学等性能的研究. 如王晓芳等[37]在摩擦材料的制备和性能研究中引入蜂巢结构, 采用结构仿生学原理, 结果发现, 蜂巢结构的引入大大改善了摩擦试样的物理性能、力学性能和摩擦磨损性能. 吴海华等[38]提出了一种制备填充型导电复合材料的方法, 采用蜂窝多孔石墨骨架, 获得新型导电复合材料, 结果表明, 蜂窝数量为18个时, 导电复合材料的电导率和抗弯强度都有明显的提升, 应用蜂窝结构后材料的力学和电学性能得到了提高.McGregor等[39]基于增材制造 (additive manufacturing, AM), 利用连续液体界面生产的方法研究六边形晶格结构的力学性能, 提出选择六边形网络为研究对象是因为其广泛的应用性和完善的理论基础[40]. 考虑六边形网格的机械零部件性能接近于预期, 他们研究发现复制材料结构的断裂模式依赖于基底材料的几何形状或者是材料的各向异性, 结果表明聚合物AM对具有晶格结构的力学零件具有很大的应用空间. 本文工作中发现具有各向异性结构的材料通过通电使其断裂也能够找到其断裂方式, 同时断裂面的性质能为蜂窝结构材料的力学性能等研究提供借鉴.

二维蜂巢结构是十分重要的晶格结构, 石墨烯等材料就具有这种二维蜂巢结构. 近年来, 随机电阻丝模型在非均匀材料断裂的数值模拟研究中被广泛应用并取得了许多十分有价值的研究成果[6].通过查阅文献发现, 单层石墨烯的实验研究大多是关于掺杂对石墨烯表面粗糙度的影响以及外部电场对石墨烯结构的影响, 少有通过实验研究手段直接给出粗糙度研究的, 因此我们的工作目前还无法与相关实验进行直接的对比分析. 本文工作的意义在于: 电阻丝模型能够很好地应用于二维蜂窝结构熔断面标度性质的研究分析, 并得出熔断面具有标度性质和奇异标度性的结论.

在本文的模拟计算过程中, 通过使用节点分析法构建了系数矩阵, 并对系数矩阵进行Cholesky分解, 然后采用 Sherman-Morrison-Woodbury算法快速对系数矩阵求逆等技术, 大大优化了计算流程和计算效率, 使得本文的数值模拟计算和分析工作得以顺利进行. 通过对粗糙度的计算, 发现熔断面呈现动力学标度性质并具有奇异标度性[41]; 通过对熔断面极值高度的分析, 发现其极值高度能很好地符合Asym2sig峰值分布函数, 至于这种分布与我们模拟过程的微观机理有什么内在联系和其特殊意义, 在我们的工作范围内, 目前还很难给出进一步明确的解释. 本文工作表明, 随机电阻丝模型不仅适用于非均匀材料断裂动力学过程的模拟, 而且也同样适用于断裂面动力学标度性质的分析.