二維隨機蜂巢網格熔斷動力學過程和熔斷面標度性質的數值模擬*

李瑞濤 唐剛 夏輝 尋之朋 李嘉翔 朱磊

(中國礦業(yè)大學物理科學與技術學院, 徐州 221116)

石墨烯等材料具有典型的二維蜂巢結構,而隨機電阻絲模型則是研究非均勻材料斷裂十分有效的統(tǒng)計物理學模型.本文嘗試對二維蜂巢結構隨機電阻絲網絡熔斷的動力學過程及熔斷面性質進行數值模擬分析,以此來研究二維非均質蜂窩材料熔斷的動力學性質和熔斷面的動力學標度性質. 模擬研究表明, 二維隨機蜂窩網格的熔斷動力學過程和熔斷面具有明顯的標度性質, 得到的熔斷面整體和局域粗糙度指數分別為和 , 這兩者之間的明顯差異表明熔斷面具有奇異標度性. 通過對熔斷面極值高度的分析發(fā)現, 熔斷面高度的極值統(tǒng)計分布能很好地滿足Asym2sig型分布, 而不是最常見的三種極值統(tǒng)計分布. 本文的研究表明, 隨機電阻絲模型在模擬非均勻材料的電流熔斷過程和熔斷表面標度性的分析中同樣適用和有效.

1 引 言

二維蜂巢結構是覆蓋二維平面的最佳拓撲結構, 其構成是由一個個正六邊形單房、房口全部朝下或朝向一邊、背對背對稱排列組合而成. 同時二維蜂巢結構也是一種十分重要的材料結構形式, 這其中就包括石墨烯等重要材料的結構. 石墨烯是由碳原子以雜化軌道組成六角型呈蜂巢狀晶格的平面薄膜, 是一種可以只有一層原子厚度的二維材料. 由于其高導電率以及獨特的電子特性, 被認為是下一代電子材料中最有前途的候選材料之一, 并具有廣泛的應用前景[1]. 因而研究其導電性能和電流熔斷機理以及熔斷面的標度屬性具有重要的理論和實踐意義.

自從2004年二維石墨烯薄膜材料被發(fā)現以來, 人們對二維石墨烯材料的結構以及導電性能都進行了廣泛和深入的研究, 并取得了很多重要的理論和實驗成果[2]. 比如, Lu等[3]通過理論模型分析,發(fā)現外部電場對研究石墨烯系統(tǒng)的原子和電子結構都有重要影響, 指出控制石墨烯電子結構的重要性, 并且發(fā)現石墨烯上的吸附原子可以作為調節(jié)電子性質的有效工具等. 2013 年, Mour和 Marder[4]通過分子動力學模擬的方法對石墨烯的斷裂力學進行了深入的研究, 建立了斷裂的幾何模型, 得到了臨界裂紋長度和應力的近似表達式, 提出了改善石墨烯韌性的方法. 他們還發(fā)現, 裂紋的路徑和產生的邊緣結構依賴于初始裂紋的長度. 2016年,Ghorbanfekr-Kalashami等[5]通過反應力場(reaction force field)的方法對石墨烯的結構和力學性能進行了研究, 結果發(fā)現摻雜物的波紋改變了石墨烯斷裂表面的粗糙度, 并且與摻雜物的數量和局部排列也有關.

在材料斷裂方面, 非均勻材料的斷裂機理和斷裂規(guī)律、斷裂表面的標度性質等近年來一直都是活躍的實驗和理論研究領域[6]. 實驗發(fā)現, 在不同載荷下的幾種材料, 例如木材、蜂窩玻璃、混凝土和紙張[7-10]等, 其斷裂表面具有分形結構和標度性質, 并具有普適性質. 材料內部結構的無序性、非均勻性和斷裂過程的非線性決定了非均勻材料的斷裂機理、斷裂面形貌及其標度性質. 在非均勻材料斷裂的微觀機理和動力學過程的理論研究中, 通常是基于晶格模型的解析近似和數值模擬. 在對非均勻材料斷裂動力學過程的數值模擬方面, 隨機電阻絲模型則是最廣泛使用和十分有效的方法[6].

隨機電阻絲模型是Arcangelis等[11]在1985年引入的. 在隨機電阻絲模型中, 是將材料看成由電阻絲組成的網格, 在網格兩端加上電壓, 用電阻絲中的電流強度來表示材料的應力, 用網格中各個電阻絲的斷裂電流閾值的隨機分布來模擬材料的非均勻性, 用電阻絲的熔斷過程來模擬實際材料的斷裂過程. 研究表明, 隨機電阻絲模型可以較準確地模擬實際材料的斷裂過程[12], 并能夠得到斷裂過程中的基本特征, 而且這樣的模型相對簡單和容易處理[13].

在隨機電阻絲模型中有兩個基本假設. 一是假設模型中電阻絲具有不可逆的熔斷性質. 要求電阻絲網絡模型滿足連續(xù)Laplace方程

變閾值范圍為

近年來, 在對隨機電阻絲模型的研究中, 科學家們做了許多十分有價值的研究工作. Duxbury等[14]發(fā)現在淬火隨機介質中電擊穿的尺寸效應,即局域斷裂理論, 證明了隨機淬火介質中的有限缺陷部分可以定性地降低實際材料的電流擊穿性能.Nukala等[15]通過數值模擬的方法分析了強無序隨機電阻絲模型的損傷成核化和局域化, 說明了斷裂的過程和特點, 找到了損傷標度律并指出損傷在大尺寸上的不相關性. 兩年后他們又分析了三維隨機電阻絲模型斷裂粗糙度和斷裂面的標度特征, 指出損傷累積是以發(fā)散的方式進行的, 直至峰值載荷處, 然后發(fā)生局域化. 同時他們還發(fā)現. 整體表面粗糙度指數與峰值載荷后損傷輪廓的局域化長度是一致的. 對不同系統(tǒng)尺寸的數據進行分析, 發(fā)現斷裂寬度分布可以很好地塌縮在一起. Toussaint和Hansen[16]還從平均場理論對柱形電阻絲網絡的斷裂機制進行了研究, 通過分離和分析系統(tǒng)的相圖, 找到了系統(tǒng)尺寸和損傷發(fā)生的特征尺寸之間的標度律. Jan ?ystein 和 Hansen[17]在研究斷裂面粗糙度的映射工作中指出, 在電阻絲模型中粗糙度指數是普遍的. 當晶格影響斷裂生長時, 電阻絲模型的粗糙度指數會隨著閾值分布而改變; 當影響消失時, 局部粗糙度指數趨于. 這些研究工作充分地證明了隨機電阻絲模型在模擬非均勻材料斷裂的微觀機理和動力學過程中的適用性和有效性[6].

在此前的模擬研究中, 因為結構簡單、計算方便而且計算量比較小, 所以采用的網絡結構大多是三角網格和菱形網絡[13], 尚未見對二維蜂巢結構直接進行模擬分析計算的工作. 二維蜂巢結構是最常見和最重要的網格結構之一, 也是石墨烯等重要材料所具有的晶格結構. 除了石墨烯以外, 如果用電阻絲斷裂的閾值代表化學鍵的強度, 則隨機電阻絲模型還可以用來研究其他的晶態(tài)材料, 對很多金屬材料的強度的模擬分析研究是十分有幫助的. 因此, 對二維蜂巢結構的隨機電阻絲網絡的微觀熔斷機理和動力學過程進行模擬研究具有一定的理論和實踐意義. 此外, 使用隨機電阻絲模型對斷裂微觀機理和動力學過程分析模擬的比較多. 本文則嘗試使用隨機電阻絲模型來分析斷裂面形貌及其標度性質.

本文對二維蜂巢結構隨機電阻絲網絡的熔斷過程進行數值模擬分析, 目的是研究二維蜂巢結構隨機電阻絲網絡斷裂的微觀機理和動力學過程以及熔斷面的形貌和標度性質. 通過對熔斷面表面寬度和局域表面寬度的計算, 發(fā)現二維蜂巢結構隨機電阻絲網絡的熔斷面呈現出標度性質并具有奇異標度性質, 得到的標度指數分別為和. 通過對二維蜂巢結構隨機電阻絲網絡熔斷面極值高度的分析, 發(fā)現其極值高度的極大(小)值呈現一定的分布規(guī)律, 分布函數則能很好地滿足Asym2sig型分布, 且同一系統(tǒng)尺寸下熔斷面的相對極大和極小高度分布具有較好的對稱性. 在模擬計算過程中, 使用節(jié)點分析法構建了系數矩陣, 并對系數矩陣進行Cholesky分解[18], 并采用 Sherman-Morrison-Woodbury 算法快速對系數矩陣求逆. 通過這些加速算法和對其結構的優(yōu)化大大地提高了計算效率, 使得本文的數值計算和分析工作能夠順利進行.

2 二維蜂巢隨機電阻絲網格及其熔斷過程的數值模擬計算方法

石墨烯中各個原子的化學性質是相同的, 但是幾何環(huán)境并不完全相同, 所以在一個基元中有兩種不同環(huán)境的原子, 構成石墨烯的復式晶格結構. 對于隨機電阻絲模型, 很多文獻選擇的是三角形、菱形[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.

最后通過初始電壓再得到一個基爾霍夫電壓方程:

聯(lián)立以上6個方程得到一個方程組:

然后繼續(xù)解方程, 斷裂, 直到所有電阻絲上的電流都為0. 這樣, 就可以得到電阻絲網的斷裂順序了.而對于解隨機電阻絲模型的節(jié)點分析法則是, 對于同樣的網絡和電流方向, 先給出節(jié)點和電流的關聯(lián)矩陣:

當電阻絲網斷裂的時候, 只要改變關聯(lián)矩陣中的數據就可以計算. 例如I1斷裂，則關聯(lián)矩陣變?yōu)?/p>

Sherman-Morrison-Woodbury算法是一種快速求逆的方法. 對于這種方法的做法是設

模擬過程中, 使用節(jié)點分析法構建系數矩陣,該系數矩陣為稀疏矩陣. 然后再對稀疏矩陣進行Cholesky分解. Cholesky分解是一種求解大型線性方程組的一種常見方法, 通過將對稱矩陣分解成一個上三角矩陣和下三角矩陣從而加速線性方程組求解速度的方法, 即. 對于上下兩個三角矩陣采用Sherman-Morrison-Woodbury算法快速求逆, 通過一些加速算法和對其結構的優(yōu)化可以大幅度加快分解速度和運算速度, 從而提高計算機的模擬效率. 在本文的模擬分析計算中, 通過使用以上方法大大地提高了計算的效率, 使得本文的模擬計算工作得以順利地進行.

3 二維蜂巢隨機電阻絲網絡熔斷面粗糙度及極值高度的計算分析

3.1 二維蜂巢隨機電阻絲網絡熔斷面粗糙度的計算分析

材料斷裂現象是一種十分復雜的隨機過程, 但其斷裂面通常都能夠形成自仿射的分形結構. 其粗糙度通常用表面寬度進行描述, 其定義為

圖3 整體表面寬度 隨系統(tǒng)尺寸 的對數-對數曲線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 二維蜂巢隨機電阻絲網絡熔斷面極值高度的分析

除熔斷面粗糙度之外, 本文還分析了熔斷面相對高度的極值分布行為. 在自然界中, 極值事件對非平衡系統(tǒng)有著非常重要的作用[21]. 很多復雜的物理系統(tǒng)大多受到極值的影響[22]. 理論研究揭示出極值事件在描述非平衡系統(tǒng)方面具有重要意義. 對于非平衡系統(tǒng), 極值統(tǒng)計早期用來分析生長表面的極值高度分布情況, 并且取得了很好的研究成果. 材料斷裂過程是一種非平衡動力學過程, 研究其極值分布具有重要的意義[23,24]. 截止目前, 文獻中很少有通過極值統(tǒng)計對其進行研究. 在表面界面生長領域, 研究發(fā)現, 對于基底尺寸為的有限系統(tǒng), 相對高度的最大值在飽和區(qū)域的分布滿足[25,26]

20 世紀 30 年代, Fisher和 Tippett[27]在對獨立同分布的極大(小)值漸進分布進行理論研究時提出了三種極值分布: Ⅰ型的Gumbel分布, Ⅱ型的 Frechet分布和Ⅲ型的 Weibull分布. 廣義Gumbel分布函數作為描述不同分布下樣本容量中極大(小)的分布, 有著非常重要的作用, 如文獻[28—30]指出, 關聯(lián)物理系統(tǒng)中可用Gumbel分布函數來描述全局漲落. 但有研究結果表明Gumbel分布函數不能很好地描述非平衡飽和表面的極值統(tǒng)計分布行為. Oliveira等[31]在研究Kardar-Parisi-Zhang (KPZ)和 Villain-Lai-Das Sarma(VLDS) 普適類的表面界面生長模型時發(fā)現, 具有不對稱局域高度分布的表面界面生長模型的極大和極小值分布不同. Wen等[32]在研究 1+1維Wolf-Villain模型飽和表面的生長高度極值統(tǒng)計分布時發(fā)現相對極小值分布不滿足廣義Gumbel分布. 而研究發(fā)現, 在很多的極值統(tǒng)計中, 也常常滿足的是Asym2sig型函數分布. 如Cui等[33]研究HSSF-CW關于示蹤物的停滯時間密度分布和Brar[34]研究光致發(fā)光譜的強度的極值分布, 在我們前面工作中[35,36], 也發(fā)現是滿足Asym2sig型分布. 本文研究發(fā)現二維蜂巢隨機電阻絲網絡熔斷面極值高度具有較好的統(tǒng)計行為, 極值高度的分布較好地符合Asym2sig型函數分布, 其表達式為

圖5 不同系統(tǒng)尺寸下石墨烯蜂巢結構隨機電阻絲網絡熔斷面相對極大高度分布Fig. 5. Relative maximum height distribution of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

圖6 不同系統(tǒng)尺寸下石墨烯蜂巢結構隨機電阻絲網絡熔斷面相對極小高度分布Fig. 6. Relative minimum height distribution of the fracture surface of random fuse model with graphene honeycomb structure under different system sizes.

圖5 為不同系統(tǒng)尺寸下石墨烯蜂巢結構的隨機電阻絲網絡熔斷面相對極大高度的概率分布圖.橫坐標定義為,表示相對極大高度的統(tǒng)計平均值,表示相對極大高度的標準差; 縱坐標定義為相對極大高度與峰值處極值高度的對應統(tǒng)計次數的比值,范圍在 [0, 1]之間. 其中, 離散的點為數值模擬結果, 實線為Asym2sig函數擬合曲線, 相關參數的擬合值如表2所列, 影響峰值寬度參數的在三種基底尺寸下擬合所得數值在誤差范圍內可認為是相等的, 說明不同的基底尺寸不影響熔斷面極值高度的分布. 圖6為不同系統(tǒng)尺寸下石墨烯蜂巢結構的隨機電阻絲網絡熔斷面相對極小高度的概率分布圖.

表2 系 統(tǒng) 尺 寸 為 L = 384, 512, 768 時Asym2sig函數擬合的參數Table 2. Parameters of Asym2sig function fitting when the system size is L = 384, 512, 768.

為了進一步說明熔斷面相對極值滿足的標度規(guī)律, 本文對縱坐標做半對數處理后發(fā)現: 在不同系統(tǒng)尺寸下石墨烯蜂巢結構的隨機電阻絲網絡熔斷面的相對極大(小)高度分布依然呈現出較好的標度規(guī)律, 如圖7和圖8所示.

很顯然, 圖7和圖8表明在不同系統(tǒng)尺寸下石墨烯蜂巢結構的隨機電阻絲網絡熔斷面的相對極大(小)高度分布滿足一定的標度規(guī)律, 很好地服從Asym2sig峰值分布函數.

圖7 不同系統(tǒng)尺寸下石墨烯蜂巢結構隨機電阻絲網絡熔斷面的相對極大高度的半對數分布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 不同系統(tǒng)尺寸下石墨烯蜂巢結構隨機電阻絲網絡熔斷面的相對極小高度的半對數分布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所示. 結果顯示, 同一系統(tǒng)尺寸下熔斷面的相對極大(小)高度分布能夠很好地重合在一起, 表明熔斷表面極大(小)高度分布具有對稱性.

圖9 系統(tǒng)尺寸 L = 384 的熔斷面的相對極大 (小)高度分布 Fig. 9. Relatively maximum (minimum) height distribution of fracture surface with system size L = 384.

圖10 系統(tǒng)尺寸 L = 512 的熔斷面的相對極大 (小)高度分布Fig. 10. Relatively maximum (minimum) height distribution of fracture surface with system size L = 512.

4 結 論

本文構建了六邊形蜂巢結構的隨機電阻絲網絡模型, 通過對電阻絲網絡施加電壓使其斷裂, 發(fā)現該模型斷裂的機理以及斷裂面的一些標度性質,從理論上豐富了隨機電阻絲模型的應用領域, 同時具有一定的實踐意義, 因為蜂窩結構廣泛應用于材料力學、電學等性能的研究. 如王曉芳等[37]在摩擦材料的制備和性能研究中引入蜂巢結構, 采用結構仿生學原理, 結果發(fā)現, 蜂巢結構的引入大大改善了摩擦試樣的物理性能、力學性能和摩擦磨損性能. 吳海華等[38]提出了一種制備填充型導電復合材料的方法, 采用蜂窩多孔石墨骨架, 獲得新型導電復合材料, 結果表明, 蜂窩數量為18個時, 導電復合材料的電導率和抗彎強度都有明顯的提升, 應用蜂窩結構后材料的力學和電學性能得到了提高.McGregor等[39]基于增材制造 (additive manufacturing, AM), 利用連續(xù)液體界面生產的方法研究六邊形晶格結構的力學性能, 提出選擇六邊形網絡為研究對象是因為其廣泛的應用性和完善的理論基礎[40]. 考慮六邊形網格的機械零部件性能接近于預期, 他們研究發(fā)現復制材料結構的斷裂模式依賴于基底材料的幾何形狀或者是材料的各向異性, 結果表明聚合物AM對具有晶格結構的力學零件具有很大的應用空間. 本文工作中發(fā)現具有各向異性結構的材料通過通電使其斷裂也能夠找到其斷裂方式, 同時斷裂面的性質能為蜂窩結構材料的力學性能等研究提供借鑒.

二維蜂巢結構是十分重要的晶格結構, 石墨烯等材料就具有這種二維蜂巢結構. 近年來, 隨機電阻絲模型在非均勻材料斷裂的數值模擬研究中被廣泛應用并取得了許多十分有價值的研究成果[6].通過查閱文獻發(fā)現, 單層石墨烯的實驗研究大多是關于摻雜對石墨烯表面粗糙度的影響以及外部電場對石墨烯結構的影響, 少有通過實驗研究手段直接給出粗糙度研究的, 因此我們的工作目前還無法與相關實驗進行直接的對比分析. 本文工作的意義在于: 電阻絲模型能夠很好地應用于二維蜂窩結構熔斷面標度性質的研究分析, 并得出熔斷面具有標度性質和奇異標度性的結論.

在本文的模擬計算過程中, 通過使用節(jié)點分析法構建了系數矩陣, 并對系數矩陣進行Cholesky分解, 然后采用 Sherman-Morrison-Woodbury算法快速對系數矩陣求逆等技術, 大大優(yōu)化了計算流程和計算效率, 使得本文的數值模擬計算和分析工作得以順利進行. 通過對粗糙度的計算, 發(fā)現熔斷面呈現動力學標度性質并具有奇異標度性[41]; 通過對熔斷面極值高度的分析, 發(fā)現其極值高度能很好地符合Asym2sig峰值分布函數, 至于這種分布與我們模擬過程的微觀機理有什么內在聯(lián)系和其特殊意義, 在我們的工作范圍內, 目前還很難給出進一步明確的解釋. 本文工作表明, 隨機電阻絲模型不僅適用于非均勻材料斷裂動力學過程的模擬, 而且也同樣適用于斷裂面動力學標度性質的分析.