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

English Abstract

全文HTML

--> --> -->

二维蜂巢结构是覆盖二维平面的最佳拓扑结构, 其构成是由一个个正六边形单房、房口全部朝下或朝向一边、背对背对称排列组合而成. 同时二维蜂巢结构也是一种十分重要的材料结构形式, 这其中就包括石墨烯等重要材料的结构. 石墨烯是由碳原子以${\rm{s}}{{\rm{p}}^2}$杂化轨道组成六角型呈蜂巢状晶格的平面薄膜, 是一种可以只有一层原子厚度的二维材料. 由于其高导电率以及独特的电子特性, 被认为是下一代电子材料中最有前途的候选材料之一, 并具有广泛的应用前景^[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方程

${\nabla ^2}V = 0.$

在节点${x_{ij}}$处电流和电压取决于连接电阻丝的电导率并遵守基尔霍夫定律和欧姆定律. 当电阻丝中的电流达到熔断电流${i_{\rm{c}}}$时, 电阻丝烧断, 其中的电流变为零, 这样的断裂过程描述的是脆性材料. 为了更加准确地描述实际非均匀材料的实际断裂过程, 需要引入对断裂行为具有重要影响的无序分布. 二是假设电阻丝具有满足某种分布函数$P(x)$的断裂阈值. 按照强度变量${\alpha _t}$和$f({\alpha _t})$的均匀阈值分布($0 \leqslant t \leqslant 1$)给出:

${f_t}\left( {{\alpha _t}} \right) = 2 - {\alpha _t} , \; 0 \leqslant {\alpha _t} \leqslant 2,$

其中${\alpha _t} = \dfrac{{\log t}}{{\log L}}$, ${f_t}({\alpha _t}) = \dfrac{{\log {L^2}tP(t)}}{{\log L}}$. 因此, 在0—1之间的均匀概率分布, 表示阈值分布$P(t)$的两个控制参数${\varPhi _0 }$, ${\varPhi _\infty }$接近于0和无穷大, 分别赋予${\varPhi _0} = 1$, ${\varPhi _\infty } = \infty $, 其中定义${\varPhi _{{0 / \infty }}} = $$ \mathop {\lim }\limits_{t \to {0 / \infty }} \left( {\dfrac{{\log \left( {tP\left( t \right)} \right)}}{{\log \left( t \right)}}} \right)$. 对于幂律阈值分布, 标度不变阈值范围为

${f_t}\left( {{\alpha _t}} \right) = 2 - \beta {\alpha _t}, \; 0 \leqslant {\alpha _t} \leqslant {2 / \beta }, $

两个控制参数为${\varPhi _0} = \beta $, ${\varPhi _\infty } = \infty $. 基于两个控制变量${\varPhi _0}$和${\varPhi _\infty }$的值, 均匀阈值分布和幂律阈值分布($\beta \leqslant 2$)都属于相同的标度值(被损伤和局域化表征). 另外, 随机电阻丝模型还需要考虑时间离散化的极限情况, 即电流的弛豫时间远远小于外电势或者电流的变化时间, 对于任何形式的局域电流超过断裂阈值的情况均予以忽略. 由此模型的动力学过程可以转化为求解${\max _{ij}} = {{{i_{t,ij}}} / {{i_{c,ij}}}}$, 这样的情况常见于极值动力学, 通常最终会发生雪崩效应.
近年来, 在对随机电阻丝模型的研究中, 科学家们做了许多十分有价值的研究工作. Duxbury等^[14]发现在淬火随机介质中电击穿的尺寸效应, 即局域断裂理论, 证明了随机淬火介质中的有限缺陷部分可以定性地降低实际材料的电流击穿性能. Nukala等^[15]通过数值模拟的方法分析了强无序随机电阻丝模型的损伤成核化和局域化, 说明了断裂的过程和特点, 找到了损伤标度律并指出损伤在大尺寸上的不相关性. 两年后他们又分析了三维随机电阻丝模型断裂粗糙度和断裂面的标度特征, 指出损伤累积是以发散的方式进行的, 直至峰值载荷处, 然后发生局域化. 同时他们还发现. 整体表面粗糙度指数与峰值载荷后损伤轮廓的局域化长度是一致的. 对不同系统尺寸的数据进行分析, 发现断裂宽度分布可以很好地塌缩在一起. Toussaint和Hansen^[16]还从平均场理论对柱形电阻丝网络的断裂机制进行了研究, 通过分离和分析系统的相图, 找到了系统尺寸和损伤发生的特征尺寸之间的标度律. Jan ?ystein和Hansen^[17]在研究断裂面粗糙度的映射工作中指出, 在电阻丝模型中粗糙度指数是普遍的. 当晶格影响断裂生长时, 电阻丝模型的粗糙度指数会随着阈值分布而改变; 当影响消失时, 局部粗糙度指数趋于${\alpha _{{\rm{loc}}}} = 0.65 \pm 0.03$. 这些研究工作充分地证明了随机电阻丝模型在模拟非均匀材料断裂的微观机理和动力学过程中的适用性和有效性^[6].
在此前的模拟研究中, 因为结构简单、计算方便而且计算量比较小, 所以采用的网络结构大多是三角网格和菱形网络^[13], 尚未见对二维蜂巢结构直接进行模拟分析计算的工作. 二维蜂巢结构是最常见和最重要的网格结构之一, 也是石墨烯等重要材料所具有的晶格结构. 除了石墨烯以外, 如果用电阻丝断裂的阈值代表化学键的强度, 则随机电阻丝模型还可以用来研究其他的晶态材料, 对很多金属材料的强度的模拟分析研究是十分有帮助的. 因此, 对二维蜂巢结构的随机电阻丝网络的微观熔断机理和动力学过程进行模拟研究具有一定的理论和实践意义. 此外, 使用随机电阻丝模型对断裂微观机理和动力学过程分析模拟的比较多. 本文则尝试使用随机电阻丝模型来分析断裂面形貌及其标度性质.
本文对二维蜂巢结构随机电阻丝网络的熔断过程进行数值模拟分析, 目的是研究二维蜂巢结构随机电阻丝网络断裂的微观机理和动力学过程以及熔断面的形貌和标度性质. 通过对熔断面表面宽度和局域表面宽度的计算, 发现二维蜂巢结构随机电阻丝网络的熔断面呈现出标度性质并具有奇异标度性质, 得到的标度指数分别为$\alpha = 0.911 \pm $$0.005$和${\alpha _{{\rm{loc}}}} = 0.808 \pm 0.003$. 通过对二维蜂巢结构随机电阻丝网络熔断面极值高度的分析, 发现其极值高度的极大(小)值呈现一定的分布规律, 分布函数则能很好地满足Asym2sig型分布, 且同一系统尺寸下熔断面的相对极大和极小高度分布具有较好的对称性. 在模拟计算过程中, 使用节点分析法构建了系数矩阵, 并对系数矩阵进行Cholesky分解^[18], 并采用Sherman-Morrison-Woodbury算法快速对系数矩阵求逆. 通过这些加速算法和对其结构的优化大大地提高了计算效率, 使得本文的数值计算和分析工作能够顺利进行.

石墨烯中各个原子的化学性质是相同的, 但是几何环境并不完全相同, 所以在一个基元中有两种不同环境的原子, 构成石墨烯的复式晶格结构. 对于随机电阻丝模型, 很多文献选择的是三角形、菱形^[13]等. 本文则尝试使用随机电阻丝模型来研究二维蜂巢网格结构的熔断动力学过程和熔断面的动力学标度性质.
如图1所示, 本文考虑基底尺寸为$L$的网络, 初始时刻所有电阻丝是完好无损的, 完整系统的总键数目为${N_{{\rm{el}}}} = (2 \times L - 1) \times L$, 且具有相同的电导率, 但每根电阻丝的阈值$t$不同, 符合[0, 1]之间的均匀阈值概率分布, 以此来模拟晶格的非均匀性. 电阻丝的熔断是不可逆的, 一旦流经电阻丝网络的电流值超过阈值$t$时, 电阻丝网络就会被烧断. 在电阻丝网络上下端施加恒定电压$V$, 为方便计算设置电压差为$\Delta V = 1$; 在水平方向上则采用周期性边界条件. 在模拟计算中, 基于基尔霍夫定律来确定每根电阻丝上的电流值. 系统中每根电阻丝$j$, 流经的电流值${i_j}$以及熔断阈值${t_j}$可以被确定, 具有最大比值${{{i_j}} / {{t_j}}}$的电阻丝会被首先烧断, 之后电阻丝网络的电流值再进行重新分配, 构建新的系数矩阵来确定下一根电阻丝的熔断. 重复这个过程直至电阻丝网络完全断开. 其中 ${{{i_j}} / {{t_j}}}$定义为电阻丝断裂的相对强度, 自变量则为基底尺寸$L$. 本文所研究的系统尺寸范围为L = 16, 32, 64, 128, 256, 384, 512, 768, 为使得模拟结果准确有效, 对各个尺寸的统计次数分别为${N_{{\rm{config}}}} = {10^7}$, ${10^6}$, $3 \times {10^5}$, $2 \times {10^4}$, ${10^4}$, $5 \times {10^3}$, $2.5 \times {10^3}$, $1.5 \times {10^3}$. 对于随机电阻丝模型的模拟计算过程, 这里以$2 \times 2$的正方格子为例进行一些简单的说明.
图 1 石墨烯蜂巢结构随机电阻丝网络通电熔断示意图
Figure1. Schematic diagram of random fuse model electric fuse in graphene honeycomb structure.

如图2所示, 共有6个不同的电流值${I_1}$—${I_6}$, 还有两个节点${A_1}$和${A_2}$. 对于电流方向, 定义从下到上, 从左到右为正. 分别对节点和每一个封闭的正方形框列基尔霍夫电流和电压方程:
图 2 2 × 2的正方格子电流流向示意图
Figure2. 2 × 2 square lattice current flow diagram.

${A_1}$: ${I_1} + {I_4} = {I_3} + {I_5}$,
${A_2}$: ${I_2} + {I_3} = {I_4} + {I_6}$;
框1: ${I_1} - {I_2} + {I_3}{\rm{ = }}0$,
框2: ${I_5} - {I_6} - {I_3}{\rm{ = }}0$.
然后根据周期性边界条件得到一个基尔霍夫电压方程:

${I_3}{\rm{ + }}{I_4}{\rm{ = }}0.$

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

${I_1}{\rm{ + }}{I_5}{\rm{ = 1}} .$

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

$\left\{ {\begin{array}{*{20}{l}} {{I_1} + {I_4} = {I_3} + {I_5},} \\ {{I_2} + {I_3} = {I_4} + {I_6},} \\ {{I_1} - {I_2} + {I_3} = 0,} \\ {{I_5} - {I_6} - {I_3} = 0,} \\ {{I_3} + {I_4} = 0,} \\ {{I_1} + {I_5} = 1.} \end{array}} \right.$

求解得到${I_1}$—${I_6}$的值, 然后通过和阈值比较决定需要断裂的电阻丝. 比如${I_1}$断裂了, 方程组需要进行修改, 修改的方法是在基尔霍夫电流方程中, ${I_1} = 0$; 在基尔霍夫电压方程中, 消去${I_1}$, 得到

$\left\{ {\begin{array}{*{20}{l}} {{I_4} = {I_3} + {I_5},} \\ {{I_2} + {I_3} = {I_4} + {I_6},} \\ {{I_5} - {I_6} - {I_3} = 0,} \\ {{I_3} + {I_4} = 0,} \\ {{I_2} - {I_3} + {I_5} = 1.} \end{array}} \right.$

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

${{C}} = \left( {\begin{array}{*{20}{c}} 1&0&{ - 1}&1&{ - 1}&0 \\ 0&1&1&{ - 1}&0&{ - 1} \end{array}} \right).$

根据节点分析法有${ {CI}} = 0$; ${ U} = {{ C}^{\rm{T}}}{{ U}_n}$, 其中${ I}$是支路电流, ${ U}$是支路电压, ${{ U}_n}$是节点电压. 经过代换得到: ${{ Y}_n}{{ U}_n} = { C}{{ I}_{\rm{s}}} - { {AY}}{{ U}_{\rm{s}}}$, ${{ Y}_n} = { {CY}}{{ C}^{\rm{T}}}$, ${{ I}_{\rm{s}}}$是电流源向量, ${{{U}}_{\rm{s}}}$是电压源向量, ${ Y}$是支路导纳矩阵. 图1中, 因为下边界电压为1, 上边界电压为0, 没有电流源, 所以${{ I}_{\rm{s}}} = 0$, ${{{U}}_{\rm{s}}} = [1{\rm{ \;\; 1 \;\; 0 \;\; 0 \;\; 0 \;\; 0}}]$, ${ Y} = $${\rm{diag}}[1{\rm{ \;\; 1 \;\; 1 \;\; 1 \;\; 1 \;\; 1}}]$. 解方程可以得到${{ U}_n}$, 进而得到各条电阻丝上的电流.
当电阻丝网断裂的时候, 只要改变关联矩阵中的数据就可以计算. 例如I₁断裂，则关联矩阵变为

${{C}} = \left( {\begin{array}{*{20}{c}} 0&0&{ - 1}&1&{ - 1}&0 \\ 0&1&1&{ - 1}&0&{ - 1} \end{array}} \right).$

在列完方程后, 即${{A}}x = {{B}}$. 在解这个线性方程组的时候, 最常见的方法是${{{A}}^{ - 1}}{{A}}x = {{{A}}^{ - 1}}{{B}}$, 这个方法是最简单的方法, 但是常规求逆这个过程计算过程比较慢.
Sherman-Morrison-Woodbury算法是一种快速求逆的方法. 对于这种方法的做法是设

$\begin{split}{{A}} =\; &\left( {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&{{a_{13}}} \\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}} \\ {{a_{31}}}&{{a_{32}}}&{{a_{33}}} \end{array}} \right) , \; \\{{B}} =\;&{{{A}}^{ - 1}} =\left( {\begin{array}{*{20}{c}} {{b_{11}}}&{{b_{12}}}&{{b_{13}}} \\ {{b_{21}}}&{{b_{22}}}&{{b_{23}}} \\ {{b_{31}}}&{{b_{32}}}&{{b_{33}}} \end{array}} \right), \\ {{ v}^{\rm{T}}} =\; & [1\;{\rm{ 0 }}\; - {\rm{1}}],\end{split}$

如果
${{C}} = \left( {\begin{array}{*{20}{c}} {{a_{11}} - {k_{ij}}}&{{a_{12}}}&{{a_{13}} + {k_{ij}}} \\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}} \\ {{a_{31}} + {k_{ij}}}&{{a_{32}}}&{{a_{33}} - {k_{ij}}} \end{array}} \right) = {{A}} - {k_{ij}}{ v}{{ v}^{\rm{T}}}$,
那么${{{C}}^{ - 1}} = \left[ {{{{A}}^{ - 1}} + {k_{ij}}\dfrac{{{ u}{{ u}^{\rm{T}}}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm T}}{ u}} \right)}}} \right]$,
其中${ u} = {{{A}}^{ - 1}}{ v}$, 那么${ x} = {{{C}}^{ - 1}}{ D}$,

${ x} = \left[ {{{{A}}^{ - 1}} + {k_{ij}}\dfrac{{{ u}{{ u}^{\rm{T}}}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm{T}}}{ u}{ u}} \right)}}} \right] *{ D}.$

在这个算法中有个技巧, 因为${ u}{{ u}^{\rm{T}}}{ D}$这个三个矩阵相乘, 通常情况下是不能交换位置的, 但是由于${ u}$和${{ u}^{\rm{T}}}$是互为转置, 所以${ u}{{ u}^{\rm{T}}}{ D} = {{ u}^{\rm{T}}}{ D}{ u}$, 这种计算方法计算速度快, 所以

${ x} = \left[ {{{{A}}^{ - 1}}{ D} + {k_{ij}}\dfrac{{{{ u}^{\rm{T}}}{ D}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm{T}}}{ u}} \right)}}{ u}} \right].$

当${{ v}^{\rm{T}}} = \left[ {0\; {\rm{ 0\; 1}}} \right]$的时候, 同样可以计算,

$\begin{split}{{C}} &= \left( {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&{{a_{13}}} \\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}} \\ {{a_{31}}}&{{a_{32}}}&{{a_{33}} - {k_{ij}}} \end{array}} \right) = {{A}} - {k_{ij}}{ v}{{ v}^{\rm{T}}}, \\ {{{C}}^{ - 1}} &= \left[ {{{{A}}^{ - 1}} + {k_{ij}}\frac{{{ u}{{ u}^{\rm{T}}}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm{T}}}{ u}} \right)}}} \right],\end{split}$

其中${ u} = {{{A}}^{ - 1}}{ v}$, 那么$x = {{{C}}^{ - 1}}{ D}$,

$\begin{split}{ x} =& \left[ {{{{A}}^{ - 1}} + {k_{ij}}\frac{{{ u}{{ u}^{\rm{T}}}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm{T}}}{ u}} \right)}}} \right]*{ D}, \\ { x} =& \left[ {{{{A}}^{ - 1}}{ D} + {k_{ij}}\frac{{{{ u}^{\rm{T}}}{ D}}}{{\left( {1 - {k_{ij}}{{ v}^{\rm T}}{ u}} \right)}}{ u}} \right].\end{split}$

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

2

-->

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

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

$W\left( {L,t} \right) = \sqrt {\frac{1}{L}\left\langle {\sum\nolimits_x {{{\left[ {h\left( {x,t} \right) - \bar h\left( {x,t} \right)} \right]}^2}} } \right\rangle } ,$

其中$L$表示基底的系统尺寸, $h\left( {x,t} \right)$表示$t$时刻在基底格点$x$处的断裂高度, $\bar h\left( {x,t} \right)$表示断裂高度的平均值. 通常情况下, 在有限尺寸系统中, 表面宽度$W\left( {L,t} \right)$满足Family-Vicsek^[19]提出的动力学标度律

$W\left( {L,t} \right) = {L^\alpha }f\left( {{t / {{L^z}}}} \right),$

其中标度函数$f\left( v \right)$具有渐近行为: 当$v < < 1$时, $f\left( v \right) \propto {v^\beta }$; 当$v \gg 1$时, $f\left( v \right)$趋于常数. (5)式中$\alpha $和$z$分别为粗糙度指数和动力学指数, 用来确定粗糙断裂表面所属的普适类^[20]. 当空间关联长度趋向$L$时, 网络完全断开, 粗糙度指数表现出$W\left( L \right) \propto $${L^\alpha }$的标度关系. 图3给出了石墨烯蜂巢结构的随机电阻丝模型熔断面整体粗糙度$W$随系统尺寸$L$变化的双对数关系, 得到的粗糙度指数满足幂律关系$W\left( L \right) \propto {L^\alpha }$, 其粗糙度指数为$\alpha = 0.911 \pm 0.005$.
图 3 整体表面宽度$W$随系统尺寸$L$的对数-对数曲线
Figure3. The log-logarithmic curve of the global surface width W with the system size L.

本文对每个系统尺寸熔断表面随机取限度为$l \left( {l \ll L} \right)$的局域窗口, 计算其局域粗糙度$w\left( {l,t} \right)$,

$w\left( {l,t} \right) = \frac{1}{{\sqrt l }}{\left\langle {\sum\nolimits_x {{{\left[ {h\left( {x,t} \right) - {{\bar h}_l}\left( t \right)} \right]}^2}} } \right\rangle ^{{1 / 2}}},$

方法同(4)式, ${\bar h_l} = \left( {{1 / l}} \right)\sum\nolimits_x {h\left( {x,t} \right)} $表示在$t$时刻尺寸为$l$的局域窗口范围内的平均高度, 如图4. 结果发现${\alpha _{{\rm{loc}}}} = 0.808 \pm 0.003$, 与整体粗糙度相比存在一定的差异. 从计算所得到的整体粗糙度和局域粗糙度指数的不同可以看出, 石墨烯蜂巢结构电阻丝网络熔断面存在奇异标度行为, 这与二维菱形结构及三角形结构^[13]的结果是不同的. 表1比较了三种结构电阻丝网络熔断面的整体及局域的标度行为.
图 4 局域表面宽度$w$随局域尺寸$l$的对数-对数曲线
Figure4. The Log-logarithmic curve of local surface width w with local size l.

模型	$\alpha $	${\alpha _{{\rm{loc}}}}$
菱形	0.752 ± 0.008	0.758 ± 0.012
三角形	0.772 ± 0.013	0.776 ± 0.003
石墨烯蜂巢结构	0.911 ± 0.005	0.808 ± 0.003

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

从表1可以看出, 二维菱形结构随机电阻丝模型断裂面的整体与局域粗糙度指数分别为$\alpha = 0.752 \pm 0.008$和${\alpha _{{\rm{loc}}}} = 0.758 \pm 0.012$, 二维三角形结构随机电阻丝模型断裂面的整体与局域粗糙度指数分别为$\alpha = 0.772 \pm 0.013$和${\alpha _{{\rm{loc}}}} = 0.776 \pm $0.003. 这两种模型结构是各向同性的, 整体与局域粗糙度的结果说明断裂面不存在奇异标度行为. 本文所研究的模型是各向异性的, 说明模型熔断表面存在着奇异标度行为.
2 -->

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

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

$P\left( {{h_{\max }},L} \right) \propto {L^{ - \alpha }}f\left( {{{{h_{\max }}} / {{L^\alpha }}}} \right),$

其中$f\left( x \right)$为标度函数, $\alpha $为整体粗糙度指数; 高度的平均值为$\left\langle h \right\rangle $; ${h_{\max }}, \; {h_{\min }}$为生长高度的最大(小)值, $R{h_{\max }} = {h_{\max }} - \left\langle h \right\rangle $, $R{h_{\min }} = \left| {{h_{\min }} - \left\langle h \right\rangle } \right|$为相对极大(小)高度. 通常, 不同系统尺寸下的相对高度极大(小)值${h_{\rm m}}$的概率分布满足以下关系:

$P\left( {{h_{\rm m}}} \right) \propto f\left( {\frac{{{h_{\rm m}} - \left\langle {{h_{\rm m}}} \right\rangle }}{\delta }} \right),$

其中$\delta $表示相对高度极大(小)值${h_{\rm m}}$的标准差, $P\left( {{h_{\rm m}}} \right){\rm{d}}{h_{\rm m}}$表示区间$\left[ {{h_{\rm m}},{h_{\rm m}} + {\rm{d}}{h_{\rm m}}} \right]$内分布的概率.
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型函数分布, 其表达式为

$\begin{split}y =\;& {y_0} + A{\left[ {1 + {{\rm{e}}^{{{ - \left( {x - {x_{\rm{c}}} + {{{\omega _1}} / 2}} \right)} / {{\omega _2}}}}}} \right]^{ - 1}}\\&\times \left\{ {1 - {{\left[ {1 + {{\rm{e}}^{{{ - \left( {x - {x_{\rm{c}}} + {{{\omega _1}} / 2}} \right)} / {{\omega _3}}}}}} \right]}^{ - 1}}} \right\},\end{split}$

其中${y_0}$为偏移量; $A$为振幅; ${x_{\rm{c}}}$为水平方向的中心点; ${\omega _1}, \; {\omega _2}, \; {\omega _3}$为宽度参量, 并且${\omega _2}, \; {\omega _3}$决定峰值位置. 其如图5和图6所示, 并且符合(7)式和(8)式所表示的标度关系.
图 5 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面相对极大高度分布
Figure5. Relative maximum height distribution of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

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

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

	384 max&min	512 max&min	768 max&min
y0	–0.001 ± 0.012	–0.002 ± 0.011	–0.001 ± 0.010
	–0.004 ± 0.008	–0.004 ± 0.007	–0.008 ± 0.007
xc	–0.57 ± 0.02	–0.58 ± 0.02	–0.59 ± 0.02
	–0.70 ± 0.02	–0.68 ± 0.02	–0.72 ± 0.02
A	1.18 ± 0.12	1.15 ± 0.08	1.28 ± 0.11
	1.38 ± 0.13	1.38 ± 0.11	1.49 ± 0.15
${\omega _1}$	0.85 ± 0.08	0.90 ± 0.06	0.80 ± 0.07
	0.65 ± 0.07	0.69 ± 0.06	0.58 ± 0.08
${\omega _2}$	0.13 ± 0.02	0.11 ± 0.01	0.11 ± 0.01
	0.11 ± 0.01	0.10 ± 0.01	0.12 ± 0.01
${\omega _3}$	0.25 ± 0.03	0.27 ± 0.03	0.32 ± 0.03
	0.28 ± 0.02	0.33 ± 0.02	0.30 ± 0.02

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

为了进一步说明熔断面相对极值满足的标度规律, 本文对纵坐标做半对数处理后发现: 在不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面的相对极大(小)高度分布依然呈现出较好的标度规律, 如图7和图8所示.
图 7 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面的相对极大高度的半对数分布
Figure7. 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 不同系统尺寸下石墨烯蜂巢结构随机电阻丝网络熔断面的相对极小高度的半对数分布
Figure8. Semi-logarithmic distribution of the relative minimum height of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

很显然, 图7和图8表明在不同系统尺寸下石墨烯蜂巢结构的随机电阻丝网络熔断面的相对极大(小)高度分布满足一定的标度规律, 很好地服从Asym2sig峰值分布函数.
本文还对同一尺寸下熔断面极值高度的极大值和极小值进行了比较, 如图9和图10所示. 结果显示, 同一系统尺寸下熔断面的相对极大(小)高度分布能够很好地重合在一起, 表明熔断表面极大(小)高度分布具有对称性.
图 9 系统尺寸L = 384的熔断面的相对极大(小)高度分布
Figure9. Relatively maximum (minimum) height distribution of fracture surface with system size L = 384.

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

本文构建了六边形蜂巢结构的随机电阻丝网络模型, 通过对电阻丝网络施加电压使其断裂, 发现该模型断裂的机理以及断裂面的一些标度性质, 从理论上丰富了随机电阻丝模型的应用领域, 同时具有一定的实践意义, 因为蜂窝结构广泛应用于材料力学、电学等性能的研究. 如王晓芳等^[37]在摩擦材料的制备和性能研究中引入蜂巢结构, 采用结构仿生学原理, 结果发现, 蜂巢结构的引入大大改善了摩擦试样的物理性能、力学性能和摩擦磨损性能. 吴海华等^[38]提出了一种制备填充型导电复合材料的方法, 采用蜂窝多孔石墨骨架, 获得新型导电复合材料, 结果表明, 蜂窝数量为18个时, 导电复合材料的电导率和抗弯强度都有明显的提升, 应用蜂窝结构后材料的力学和电学性能得到了提高. McGregor等^[39]基于增材制造(additive manufacturing, AM), 利用连续液体界面生产的方法研究六边形晶格结构的力学性能, 提出选择六边形网络为研究对象是因为其广泛的应用性和完善的理论基础^[40]. 考虑六边形网格的机械零部件性能接近于预期, 他们研究发现复制材料结构的断裂模式依赖于基底材料的几何形状或者是材料的各向异性, 结果表明聚合物AM对具有晶格结构的力学零件具有很大的应用空间. 本文工作中发现具有各向异性结构的材料通过通电使其断裂也能够找到其断裂方式, 同时断裂面的性质能为蜂窝结构材料的力学性能等研究提供借鉴.
二维蜂巢结构是十分重要的晶格结构, 石墨烯等材料就具有这种二维蜂巢结构. 近年来, 随机电阻丝模型在非均匀材料断裂的数值模拟研究中被广泛应用并取得了许多十分有价值的研究成果^[6]. 通过查阅文献发现, 单层石墨烯的实验研究大多是关于掺杂对石墨烯表面粗糙度的影响以及外部电场对石墨烯结构的影响, 少有通过实验研究手段直接给出粗糙度研究的, 因此我们的工作目前还无法与相关实验进行直接的对比分析. 本文工作的意义在于: 电阻丝模型能够很好地应用于二维蜂窝结构熔断面标度性质的研究分析, 并得出熔断面具有标度性质和奇异标度性的结论.
在本文的模拟计算过程中, 通过使用节点分析法构建了系数矩阵, 并对系数矩阵进行Cholesky分解, 然后采用Sherman-Morrison-Woodbury算法快速对系数矩阵求逆等技术, 大大优化了计算流程和计算效率, 使得本文的数值模拟计算和分析工作得以顺利进行. 通过对粗糙度的计算, 发现熔断面呈现动力学标度性质并具有奇异标度性^[41]; 通过对熔断面极值高度的分析, 发现其极值高度能很好地符合Asym2sig峰值分布函数, 至于这种分布与我们模拟过程的微观机理有什么内在联系和其特殊意义, 在我们的工作范围内, 目前还很难给出进一步明确的解释. 本文工作表明, 随机电阻丝模型不仅适用于非均匀材料断裂动力学过程的模拟, 而且也同样适用于断裂面动力学标度性质的分析.